Commit 43d81252 authored by frey_m's avatar frey_m
Browse files

Merge branch 'master' of gitlab.psi.ch:OPAL/src into 554-unused-classes-in-src

parents c8f39cb9 4848122f
//
// Class CavityAutophaser
//
// This class determines the phase of an RF cavity for which the reference particle
// is accelerated to the highest energy.
//
// Copyright (c) 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Algorithms/CavityAutophaser.h"
#include "Algorithms/Vektor.h"
#include "AbsBeamline/RFCavity.h"
......@@ -19,8 +40,7 @@ CavityAutophaser::CavityAutophaser(const PartData &ref,
itsCavity_m(cavity)
{
double zbegin = 0.0, zend = 0.0;
PartBunchBase<double, 3> *fakeBunch = NULL;
cavity->initialise(fakeBunch, zbegin, zend);
cavity->getDimensions(zbegin, zend);
initialR_m = Vector_t(0, 0, zbegin);
}
......@@ -138,10 +158,10 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
itsCavity_m->getName() + "_AP.dat"
});
std::ofstream out(fname);
track(initialR_m, initialP_m, t + tErr, dt, newPhase, &out);
track(t + tErr, dt, newPhase, &out);
out.close();
} else {
track(initialR_m, initialP_m, t + tErr, dt, newPhase, NULL);
track(t + tErr, dt, newPhase, NULL);
}
INFOMSG(level1 << std::fixed << std::setprecision(4)
......@@ -209,7 +229,7 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
if (element->getAutophaseVeto()) {
double basePhase = std::fmod(element->getFrequencym() * t, Physics::two_pi);
double phase = std::fmod(originalPhase - basePhase + Physics::two_pi, Physics::two_pi);
double E = track(initialR_m, initialP_m, t, dt, phase);
double E = track(t, dt, phase);
std::pair<double, double> status(originalPhase, E);//-basePhase, E);
return status;
}
......@@ -220,7 +240,7 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
const int numRefinements = Options::autoPhase;
int j = -1;
double E = track(initialR_m, initialP_m, t, dt, phi);
double E = track(t, dt, phi);
double Emax = E;
do {
......@@ -228,7 +248,7 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
Emax = E;
initialPhase = phi;
phi -= dphi;
E = track(initialR_m, initialP_m, t, dt, phi);
E = track(t, dt, phi);
} while(E > Emax);
if(j == 0) {
......@@ -240,20 +260,20 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
Emax = E;
initialPhase = phi;
phi += dphi;
E = track(initialR_m, initialP_m, t, dt, phi);
E = track(t, dt, phi);
} while(E > Emax);
}
for(int rl = 0; rl < numRefinements; ++ rl) {
dphi /= 2.;
phi = initialPhase - dphi;
E = track(initialR_m, initialP_m, t, dt, phi);
E = track(t, dt, phi);
if(E > Emax) {
initialPhase = phi;
Emax = E;
} else {
phi = initialPhase + dphi;
E = track(initialR_m, initialP_m, t, dt, phi);
E = track(t, dt, phi);
if(E > Emax) {
initialPhase = phi;
Emax = E;
......@@ -262,15 +282,13 @@ std::pair<double, double> CavityAutophaser::optimizeCavityPhase(double initialPh
}
Phimax = std::fmod(initialPhase + Physics::two_pi, Physics::two_pi);
E = track(initialR_m, initialP_m, t, dt, Phimax + originalPhase);
E = track(t, dt, Phimax + originalPhase);
std::pair<double, double> status(Phimax, E);
return status;
}
double CavityAutophaser::track(Vector_t /*R*/,
Vector_t /*P*/,
double t,
double CavityAutophaser::track(double t,
const double dt,
const double phase,
std::ofstream *out) const {
......
//
// Class CavityAutophaser
//
// This class determines the phase of an RF cavity for which the reference particle
// is accelerated to the highest energy.
//
// Copyright (c) 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef CAVITYAUTOPHASER
#define CAVITYAUTOPHASER
......@@ -22,9 +43,7 @@ private:
double t,
double dt);
double track(Vector_t R,
Vector_t P,
double t,
double track(double t,
const double dt,
const double phase,
std::ofstream *out = NULL) const;
......
//
// Class OrbitThreader
//
// This class determines the design path by tracking the reference particle through
// the 3D lattice.
//
// Copyright (c) 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/CavityAutophaser.h"
......@@ -31,7 +53,6 @@ OrbitThreader::OrbitThreader(const PartData &ref,
double maxDiffZBunch,
double t,
double dt,
size_t maxIntegSteps,
double zstop,
OpalBeamline &bl) :
r_m(r),
......@@ -39,7 +60,6 @@ OrbitThreader::OrbitThreader(const PartData &ref,
pathLength_m(s),
time_m(t),
dt_m(dt),
maxIntegSteps_m(maxIntegSteps),
zstop_m(zstop),
itsOpalBeamline_m(bl),
errorFlag_m(0),
......@@ -82,16 +102,16 @@ OrbitThreader::OrbitThreader(const PartData &ref,
}
distTrackBack_m = std::min(pathLength_m, std::max(0.0, maxDiffZBunch));
computeBoundingBox();
}
void OrbitThreader::execute() {
double initialPathLength = pathLength_m;
double maxDistance = computeMaximalImplicitDrift();
auto allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
std::set<std::string> visitedElements;
trackBack(2 * maxDistance);
trackBack();
Vector_t nextR = r_m / (Physics::c * dt_m);
integrator_m.push(nextR, p_m, dt_m);
......@@ -109,7 +129,8 @@ void OrbitThreader::execute() {
double initialS = pathLength_m;
Vector_t initialR = r_m;
Vector_t initialP = p_m;
integrate(elementSet, maxIntegSteps_m, 2 * maxDistance);
double maxDistance = computeDriftLengthToBoundingBox(elementSet, r_m, p_m);
integrate(elementSet, maxDistance);
registerElement(elementSet, initialS, initialR, initialP);
......@@ -146,7 +167,7 @@ void OrbitThreader::execute() {
processElementRegister();
}
void OrbitThreader::integrate(const IndexMap::value_t &activeSet, size_t /*maxSteps*/, double maxDrift) {
void OrbitThreader::integrate(const IndexMap::value_t &activeSet, double maxDrift) {
static size_t step = 0;
CoordinateSystemTrafo labToBeamline = itsOpalBeamline_m.getCSTrafoLab2Local();
const double oldPathLength = pathLength_m;
......@@ -281,7 +302,7 @@ double OrbitThreader::getMaxDesignEnergy(const IndexMap::value_t &elementSet) co
return designEnergy;
}
void OrbitThreader::trackBack(double maxDrift) {
void OrbitThreader::trackBack() {
dt_m *= -1;
double initialPathLength = pathLength_m;
......@@ -289,11 +310,12 @@ void OrbitThreader::trackBack(double maxDrift) {
integrator_m.push(nextR, p_m, dt_m);
nextR *= Physics::c * dt_m;
maxDrift = std::min(maxDrift, distTrackBack_m);
while (std::abs(initialPathLength - pathLength_m) < distTrackBack_m) {
auto elementSet = itsOpalBeamline_m.getElements(nextR);
integrate(elementSet, 1000, maxDrift);
double maxDrift = computeDriftLengthToBoundingBox(elementSet, r_m, -p_m);
maxDrift = std::min(maxDrift, distTrackBack_m);
integrate(elementSet, maxDrift);
nextR = r_m / (Physics::c * dt_m);
integrator_m.push(nextR, p_m, dt_m);
......@@ -388,75 +410,30 @@ void OrbitThreader::setDesignEnergy(FieldList &allElements, const std::set<std::
}
}
double OrbitThreader::computeMaximalImplicitDrift() {
void OrbitThreader::computeBoundingBox() {
FieldList allElements = itsOpalBeamline_m.getElementByType(ElementBase::ANY);
double maxDrift = 0.0;
MarkerRep start("#S");
CoordinateSystemTrafo toEdge(r_m, getQuaternion(p_m, Vector_t(0, 0, 1)));
start.setElementLength(0.0);
start.setCSTrafoGlobal2Local(toEdge);
std::shared_ptr<Component> startPtr(static_cast<Marker*>(start.clone()));
allElements.push_front(ClassicField(startPtr, 0.0, 0.0));
FieldList::iterator it = allElements.begin();
const FieldList::iterator end = allElements.end();
globalBoundingBox_m.lowerLeftCorner = Vector_t(std::numeric_limits<double>::max());
globalBoundingBox_m.upperRightCorner = Vector_t(-std::numeric_limits<double>::max());
for (; it != end; ++ it) {
auto element1 = it->getElement();
const auto &toEdge = element1->getCSTrafoGlobal2Local();
auto toEnd = element1->getEdgeToEnd() * toEdge;
Vector_t end1 = toEnd.transformFrom(Vector_t(0, 0, 0));
Vector_t directionEnd = toEnd.rotateFrom(Vector_t(0, 0, 1));
if (element1->getType() == ElementBase::RBEND ||
element1->getType() == ElementBase::SBEND ||
element1->getType() == ElementBase::RBEND3D ) {
auto toBegin = element1->getEdgeToBegin() * toEdge;
BendBase *bend = static_cast<BendBase*>(element1.get());
double angleRelativeToFace = bend->getEntranceAngle() - bend->getBendAngle();
directionEnd = toBegin.rotateFrom(Vector_t(sin(angleRelativeToFace), 0, cos(angleRelativeToFace)));
}
double minDistanceLocal = std::numeric_limits<double>::max();
FieldList::iterator it2 = allElements.begin();
for (; it2 != end; ++ it2) {
if (it == it2) continue;
auto element2 = it2->getElement();
const auto &toEdge = element2->getCSTrafoGlobal2Local();
auto toBegin = element2->getEdgeToBegin() * toEdge;
auto toEnd = element2->getEdgeToEnd() * toEdge;
Vector_t begin2 = toBegin.transformFrom(Vector_t(0, 0, 0));
Vector_t end2 = toEnd.transformFrom(Vector_t(0, 0, 0));
Vector_t directionBegin = toBegin.rotateFrom(Vector_t(0, 0, 1));
if (element2->getType() == ElementBase::RBEND ||
element2->getType() == ElementBase::SBEND ||
element2->getType() == ElementBase::RBEND3D ) {
BendBase *bend = static_cast<BendBase*>(element2.get());
double E1 = bend->getEntranceAngle();
directionBegin = toBegin.rotateFrom(Vector_t(sin(E1), 0, cos(E1)));
}
double distance = euclidean_norm(begin2 - end1);
double directionProjection = dot(directionEnd, directionBegin);
bool overlapping = dot(begin2 - end1, directionBegin) < 0.0? true: false;
if (!overlapping &&
directionProjection > 0.999 &&
minDistanceLocal > distance) {
minDistanceLocal = distance;
}
}
if (maxDrift < minDistanceLocal &&
minDistanceLocal != std::numeric_limits<double>::max()) {
maxDrift = minDistanceLocal;
if (it->getElement()->getType() == ElementBase::MARKER) {
continue;
}
ElementBase::BoundingBox other = it->getBoundingBoxInLabCoords();
globalBoundingBox_m.getCombinedBoundingBox(other);
}
}
maxDrift = std::min(maxIntegSteps_m * std::abs(dt_m) * Physics::c, maxDrift);
double OrbitThreader::computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
const Vector_t & position,
const Vector_t & direction) const {
if (elements.empty() ||
(elements.size() == 1 && (*elements.begin())->getType() == ElementBase::ElementType::DRIFT)) {
boost::optional<Vector_t> intersectionPoint = globalBoundingBox_m.getPointOfIntersection(position, direction);
return intersectionPoint ? euclidean_norm(intersectionPoint.get() - position): 10.0;
}
return maxDrift;
}
return std::numeric_limits<double>::max();
}
\ No newline at end of file
//
// Class OrbitThreader
//
// This class determines the design path by tracking the reference particle through
// the 3D lattice.
//
// Copyright (c) 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020 Christof Metzger-Kraus
//
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#ifndef OPAL_ORBITTHREADER_H
#define OPAL_ORBITTHREADER_H
......@@ -20,7 +41,6 @@ public:
double maxDiffZBunch,
double t,
double dT,
size_t maxIntegSteps,
double zstop,
OpalBeamline &bl);
......@@ -48,8 +68,6 @@ private:
/// the time step
double dt_m;
/// the number of time steps to track
const size_t maxIntegSteps_m;
/// final position in path length
const double zstop_m;
......@@ -64,6 +82,8 @@ private:
std::ofstream logger_m;
size_t loggingFrequency_m;
ElementBase::BoundingBox globalBoundingBox_m;
struct elementPosition {
double startField_m;
double endField_m;
......@@ -78,8 +98,8 @@ private:
std::multimap<std::string, elementPosition> elementRegistry_m;
void trackBack(double maxDrift);
void integrate(const IndexMap::value_t &activeSet, size_t maxSteps, double maxDrift = 10.0);
void trackBack();
void integrate(const IndexMap::value_t &activeSet, double maxDrift = 10.0);
bool containsCavity(const IndexMap::value_t &activeSet);
void autophaseCavities(const IndexMap::value_t &activeSet, const std::set<std::string> &visitedElements);
double getMaxDesignEnergy(const IndexMap::value_t &elementSet) const;
......@@ -87,7 +107,11 @@ private:
void registerElement(const IndexMap::value_t &elementSet, double, const Vector_t &r, const Vector_t &p);
void processElementRegister();
void setDesignEnergy(FieldList &allElements, const std::set<std::string> &visitedElements);
double computeMaximalImplicitDrift();
void computeBoundingBox();
// double computeMaximalImplicitDrift();
double computeDriftLengthToBoundingBox(const std::set<std::shared_ptr<Component>> & elements,
const Vector_t & position,
const Vector_t & direction) const;
};
inline
......
......@@ -4,7 +4,9 @@
// The visitor class for tracking particles with time as independent
// variable.
//
// Copyright (c) 200x - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 200x - 2014, Christof Kraus, Paul Scherrer Institut, Villigen PSI, Switzerland
// 2015 - 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020, Christof Metzger-Kraus
// All rights reserved
//
// This file is part of OPAL.
......@@ -75,8 +77,7 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
fieldEvaluationTimer_m(IpplTimings::getTimer("External field eval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
particleMatterStatus_m(false),
totalParticlesInSimulation_m(0)
particleMatterStatus_m(false)
{}
ParallelTTracker::ParallelTTracker(const Beamline &beamline,
......@@ -108,8 +109,7 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
fieldEvaluationTimer_m(IpplTimings::getTimer("External field eval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
particleMatterStatus_m(false),
totalParticlesInSimulation_m(0)
particleMatterStatus_m(false)
{
for (unsigned int i = 0; i < zstop.size(); ++ i) {
stepSizes_m.push_back(dt[i], zstop[i], maxSteps[i]);
......@@ -193,7 +193,6 @@ void ParallelTTracker::execute() {
prepareSections();
double minTimeStep = stepSizes_m.getMinTimeStep();
unsigned long long totalNumSteps = stepSizes_m.getNumStepsFinestResolution();
itsOpalBeamline_m.activateElements();
......@@ -248,7 +247,6 @@ void ParallelTTracker::execute() {
-rmin(2),
itsBunch_m->getT(),
(back_track? -minTimeStep: minTimeStep),
totalNumSteps,
stepSizes_m.getFinalZStop() + 2 * rmax(2),
itsOpalBeamline_m);
......@@ -257,7 +255,6 @@ void ParallelTTracker::execute() {
saveCavityPhases();
numParticlesInSimulation_m = itsBunch_m->getTotalNum();
totalParticlesInSimulation_m = itsBunch_m->getTotalNum();
setTime();
......@@ -588,7 +585,6 @@ void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
ne = itsBunch_m->destroyT();
}
numParticlesInSimulation_m = itsBunch_m->getTotalNum();
totalParticlesInSimulation_m -= ne;
deletedParticles_m = true;
}
......@@ -599,7 +595,7 @@ void ParallelTTracker::computeExternalFields(OrbitThreader &oth) {
if (ne > 0) {
msg << level1 << "* Deleted " << ne << " particles, "
<< "remaining " << totalParticlesInSimulation_m << " particles" << endl;
<< "remaining " << numParticlesInSimulation_m << " particles" << endl;
}
}
......@@ -1390,4 +1386,4 @@ void ParallelTTracker::evenlyDistributeParticles() {
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End:
// End:
\ No newline at end of file
......@@ -4,7 +4,9 @@
// The visitor class for tracking particles with time as independent
// variable.
//
// Copyright (c) 200x - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// Copyright (c) 200x - 2014, Christof Kraus, Paul Scherrer Institut, Villigen PSI, Switzerland
// 2015 - 2016, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin, Germany
// 2017 - 2020, Christof Metzger-Kraus
// All rights reserved
//
// This file is part of OPAL.
......@@ -220,8 +222,6 @@ private:
std::set<ParticleMatterInteractionHandler*> activeParticleMatterInteractionHandlers_m;
bool particleMatterStatus_m;
unsigned long totalParticlesInSimulation_m;
/********************** END VARIABLES ***********************************/
void kickParticles(const BorisPusher &pusher);
......
......@@ -288,6 +288,8 @@ Option::Option(const std::string &name, Option *parent):
Attributes::setReal(itsAttr[MTSSUBSTEPS], mtsSubsteps);
Attributes::setReal(itsAttr[REMOTEPARTDEL], remotePartDel);
Attributes::setReal(itsAttr[REPARTFREQ], repartFreq);
Attributes::setReal(itsAttr[MINBINEMITTED], minBinEmitted);
Attributes::setReal(itsAttr[MINSTEPFORREBIN], minStepForRebin);
Attributes::setReal(itsAttr[REBINFREQ], rebinFreq);
Attributes::setBool(itsAttr[RHODUMP], rhoDump);
Attributes::setBool(itsAttr[EBDUMP], ebDump);
......@@ -428,6 +430,14 @@ void Option::execute() {
repartFreq = int(Attributes::getReal(itsAttr[REPARTFREQ]));
}
if (itsAttr[MINBINEMITTED]) {
minBinEmitted = int(Attributes::getReal(itsAttr[MINBINEMITTED]));
}
if (itsAttr[MINSTEPFORREBIN]) {
minStepForRebin = int(Attributes::getReal(itsAttr[MINSTEPFORREBIN]));
}
if(itsAttr[REBINFREQ]) {
rebinFreq = int(Attributes::getReal(itsAttr[REBINFREQ]));
}
......
......@@ -692,31 +692,22 @@ void Bend2D::findBendStrength() {
if (error < tolerance)
return;
double fieldStep = estimateFieldAdjustmentStep(actualBendAngle);
double fieldStep = std::copysign(1.0, fieldAmplitude_m) * std::abs(estimateFieldAdjustmentStep(actualBendAngle));
double amplitude1 = fieldAmplitude_m;
double amplitude2 = amplitude1;
double bendAngle1 = actualBendAngle;
double amplitude2 = fieldAmplitude_m + fieldStep;
fieldAmplitude_m = amplitude2;
double bendAngle2 = calculateBendAngle();
if (std::abs(bendAngle1) > std::abs(angle_m)) {
while (std::abs(bendAngle2) > std::abs(angle_m)) {
amplitude1 = amplitude2;
bendAngle1 = bendAngle2;
amplitude2 += fieldStep;
fieldAmplitude_m = amplitude2;
bendAngle2 = calculateBendAngle();
}
} else {
while (std::abs(bendAngle2) < std::abs(angle_m)) {
amplitude1 = amplitude2;
bendAngle1 = bendAngle2;
amplitude2 += fieldStep;
fieldAmplitude_m = amplitude2;
bendAngle2 = calculateBendAngle();
double bendAngle2 = bendAngle1;
int stepSign = std::abs(bendAngle1) > std::abs(angle_m) ? -1: 1;
while (true) {
amplitude1 = amplitude2;
bendAngle1 = bendAngle2;
amplitude2 += stepSign * fieldStep;
fieldAmplitude_m = amplitude2;
bendAngle2 = calculateBendAngle();
if (stepSign * (std::abs(bendAngle2) - std::abs(angle_m)) > 0) {
break;
}
}
......@@ -1689,4 +1680,32 @@ std::array<double,2> Bend2D::getExitFringeFieldLength() const {
extFFL[1] = ( exitParameter2_m-exitParameter1_m ); //before edge
return extFFL;
}
ElementBase::BoundingBox Bend2D::getBoundingBoxInLabCoords() const {
CoordinateSystemTrafo toBegin = getEdgeToBegin() * csTrafoGlobal2Local_m;
CoordinateSystemTrafo toEnd = getEdgeToEnd() * csTrafoGlobal2Local_m;
std::vector<Vector_t> outline = getOutline();
BoundingBox bb;
bb.lowerLeftCorner = std::numeric_limits<double>::max();
bb.upperRightCorner = std::numeric_limits<double>::lowest();
Vector_t dY(0, 0.5 * getFullGap(), 0);
for (int i : {-1, 1}) {
for (const Vector_t & vec: outline) {
Vector_t vecInLabCoords = csTrafoGlobal2Local_m.transformFrom(vec + i * dY);
for (unsigned int d = 0; d < 3; ++ d) {
if (vecInLabCoords(d) < bb.lowerLeftCorner(d)) {
bb.lowerLeftCorner(d) = vecInLabCoords(d);
}