Commit 64d9a708 authored by adelmann's avatar adelmann 🎗

Merge branch 'OPAL-maps' into 'master'

Opal maps

See merge request !37
parents d5af9f6a 611ab378
......@@ -51,7 +51,7 @@ PROJECT_BRIEF = "OPAL"
# and the maximum width should not exceed 200 pixels. Doxygen will copy the logo
# to the output directory.
PROJECT_LOGO =
PROJECT_LOGO = ./amas.png
# The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) path
# into which the generated documentation will be written. If a relative path is
......@@ -1572,7 +1572,7 @@ EXTRA_SEARCH_MAPPINGS =
# If the GENERATE_LATEX tag is set to YES doxygen will generate LaTeX output.
# The default value is: YES.
GENERATE_LATEX = NO
GENERATE_LATEX = YES
# The LATEX_OUTPUT tag is used to specify where the LaTeX docs will be put. If a
# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
......@@ -1675,7 +1675,7 @@ PDF_HYPERLINKS = YES
# The default value is: YES.
# This tag requires that the tag GENERATE_LATEX is set to YES.
USE_PDFLATEX = NO
USE_PDFLATEX = YES
# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \batchmode
# command to the generated LaTeX files. This will instruct LaTeX to keep running
......@@ -2328,4 +2328,4 @@ GENERATE_LEGEND = YES
# The default value is: YES.
# This tag requires that the tag HAVE_DOT is set to YES.
DOT_CLEANUP = YES
\ No newline at end of file
DOT_CLEANUP = YES
A note on documentation etc. for collegues starting work on OPAL:
1) The full documentation gor OPAL 2.x.x is on the Wiki: https://gitlab.psi.ch/OPAL/src/wikis/Manual/home
2) Also the C++ source code must be documented please read: http://en.wikipedia.org/wiki/Doxygen and document
your source code accordingly.
3) We have to do some artwork for our paper(s) and documentation. Please use only TikZ and see the power of
this package (and many examples to start) at http://www.texample.net/tikz/examples/ .
%% objective - bibliography file for femaxx reference card
%% author - benedikt oswald
%% modified - 2008 jun 11, creation, benedikt oswald
@Article{arbenzetal2001,
author = {Arbenz, P. and Geus, R. and Adam, S.},
title = {Solving Maxwell eigenvalue problems for accelerating cavities},
journal = {Physical Review Special Topics - Accelerators and Beams},
year = {2001},
volume = {4},
pages = {022001-1:022001-10},
annote = {Abstract: we investigate algorithms for computing steady state electromagnetic waves
in cavities. The Maxwell equations for the strength of the electric field are solved by
a mixed method with quadratic finite edge (N{\'e}d{\'e}lec) elements for the field
values and corresponding node-based finite elements for the Lagrange multiplier, This
approach avoids so-called spurious modes which are introduced if the divergence-free
condition for the electric field is not treated properly. To compute a few of the smalles
positive eigenvalues and corresponding eigenmodes of the resulting large sparse-matrix
eigenvalue problems, two algorithms have been used: the implicitly restarted Lanczos
algorithm and the Jacobi-Davidson algorithm, both with shit-and-invert spectral transformation.
Two-level hierarchical basis preconditioners have been employed for the iterative solution
of the resulting systems of equations.}
}
@PhdThesis{geus2002,
author = {Geus, R.},
title = {The Jacobi-Davidson Algorithm for solving large sparse symmetric eigenvalue problems with application to the design of accelerator cavities},
school = {Swiss Federal Institute of Technology Zurich},
year = {2002},
OPTkey = {Diss. ETH NO. 14734},
type = {},
address = {},
month = {},
note = {},
annote = {}
}
This diff is collapsed.
This diff is collapsed.
......@@ -22,6 +22,10 @@
#include "Steppers/BorisPusher.h"
#include "Structure/DataSink.h"
#include "Classic/FixedAlgebra/FTps.h"
#include "Classic/FixedAlgebra/FVps.h"
#include "Algorithms/OrbitThreader.h"
#include "Algorithms/IndexMap.h"
#include "AbsBeamline/AlignWrapper.h"
......@@ -52,6 +56,12 @@
#include "Elements/OpalBeamline.h"
#include <cmath>
#define PSdim 6
#define PHIL_WRITE 1
class BMultipoleField;
template <class T, unsigned Dim>
......@@ -124,8 +134,9 @@ public:
const std::vector<unsigned long long> &maxSTEPS,
double zstart,
const std::vector<double> &zstop,
const std::vector<double> &dt);
const std::vector<double> &dt,
const int& truncOrder);
virtual ~ThickTracker();
......@@ -216,12 +227,56 @@ public:
void autophaseCavities(const BorisPusher &pusher);
void findStartPosition(const BorisPusher &pusher);
typedef FTps<double, PSdim> series_t;
typedef FVps<double, PSdim> map_t;
typedef FMatrix<double, PSdim, PSdim> fMatrix_t;
typedef FMatrix<std::complex<double>, PSdim, PSdim> cfMatrix_t;
series_t x, px, y, py, delta; //z
struct structMapTracking {
std::string elementName;
map_t elementMap;
std::size_t nSlices;
double elementPos;
double stepSize;
};
series_t createHamiltonian(std::shared_ptr<Component> element, double& stepSize, std::size_t& nSlices);
void defMapTrackingElement(std::shared_ptr<Component> element, structMapTracking& elSrct, std::list<structMapTracking>& mBL);
void fillDrift(std::list<structMapTracking>& mapBeamLine,double& elementPos, double& undefSpace);
void setHamiltonianDrift(series_t& H, double& beta0, double& gamma0);
void setHamiltonianSBend(series_t& H, double& beta0, double& gamma0, double& q, double& h, double& K0 );
void setHamiltonianRBend(series_t& H, double& beta0, double& gamma0, double& q, double& h, double& K0 );
void setHamiltonianQuadrupole(series_t& H, double& beta0, double& gamma0, double& q, double& K1 );
void eigenDecomp(fMatrix_t& M, cfMatrix_t& eigenVal, cfMatrix_t& eigenVec, cfMatrix_t& invEigenVec);
void linTAnalyze(fMatrix_t& tMatrix);
void linSigAnalyze();
cfMatrix_t getBlockDiagonal(fMatrix_t& M, cfMatrix_t& eigenVecM, cfMatrix_t& invEigenVecM);
void dumpStats(long long step, bool psDump, bool statDump);
void writePhaseSpace(const long long step, bool psDump, bool statDump);
void setTime();
private:
// Not implemented.
ThickTracker();
ThickTracker(const ThickTracker &);
void operator=(const ThickTracker &);
void trackParticles_m(
#ifdef PHIL_WRITE
std::ofstream& outfile,
#endif
const std::list<structMapTracking>& mapBeamLine);
// Fringe fields for entrance and exit of magnetic elements.
void applyEntranceFringe(double edge, double curve,
......@@ -254,7 +309,9 @@ private:
CoordinateSystemTrafo referenceToLabCSTrafo_m;
bool globalEOL_m;
int truncOrder_m;
};
inline void ThickTracker::visitAlignWrapper(const AlignWrapper &wrap) {
......@@ -366,4 +423,4 @@ inline void ThickTracker::visitCyclotronValley(const CyclotronValley &cv) {
itsOpalBeamline_m.visit(cv, *this, itsBunch_m);
}
#endif // OPAL_ThickTracker_HH
\ No newline at end of file
#endif // OPAL_ThickTracker_HH
......@@ -71,7 +71,8 @@ Bend::Bend():
cosEntranceAngle_m(1.0),
sinEntranceAngle_m(0.0),
tanEntranceAngle_m(0.0),
tanExitAngle_m(0.0) {
tanExitAngle_m(0.0),
nSlices_m(1){
setElType(isDipole);
......@@ -111,7 +112,8 @@ Bend::Bend(const Bend &right):
cosEntranceAngle_m(right.cosEntranceAngle_m),
sinEntranceAngle_m(right.sinEntranceAngle_m),
tanEntranceAngle_m(right.tanEntranceAngle_m),
tanExitAngle_m(right.tanExitAngle_m) {
tanExitAngle_m(right.tanExitAngle_m),
nSlices_m(right.nSlices_m){
setElType(isDipole);
......@@ -151,7 +153,8 @@ Bend::Bend(const std::string &name):
cosEntranceAngle_m(1.0),
sinEntranceAngle_m(0.0),
tanEntranceAngle_m(0.0),
tanExitAngle_m(0.0) {
tanExitAngle_m(0.0),
nSlices_m(1){
setElType(isDipole);
......@@ -1851,4 +1854,14 @@ void Bend::setupFringeWidths()
{
widthEntranceFringe_m = 2 * std::min(entranceParameter3_m - entranceParameter1_m, aperture_m.second[0]) + aperture_m.second[0];
widthExitFringe_m = 2 * std::min(exitParameter3_m - exitParameter1_m, aperture_m.second[0]) + aperture_m.second[0];
}
\ No newline at end of file
}
//set the number of slices for map tracking
void Bend::setNSlices(const std::size_t& nSlices) {
nSlices_m = nSlices;
}
//get the number of slices for map tracking
std::size_t Bend::getNSlices() const {
return nSlices_m;
}
......@@ -131,6 +131,14 @@ public:
virtual CoordinateSystemTrafo getBeginToEnd() const;
virtual bool isInside(const Vector_t &r) const;
//set number of slices for map tracking
void setNSlices(const std::size_t& nSlices);
//set number of slices for map tracking
std::size_t getNSlices() const;
protected:
void setMessageHeader(const std::string & header);
double getStartField() const;
......@@ -294,6 +302,8 @@ private:
CoordinateSystemTrafo computeAngleTrafo_m;
double maxAngle_m;
std::size_t nSlices_m;
};
......@@ -415,4 +425,4 @@ Vector_t Bend::transformToExitRegion(const Vector_t &R) const {
return toExitRegion_m.transformTo(R);
}
#endif // CLASSIC_BEND_H
\ No newline at end of file
#endif // CLASSIC_BEND_H
......@@ -28,19 +28,22 @@ extern Inform *gmsg;
// ------------------------------------------------------------------------
Drift::Drift():
Component()
Component(),
nSlices_m(1)
{ }
Drift::Drift(const Drift &right):
Component(right)
Component(right),
nSlices_m(right.nSlices_m)
{ }
Drift::Drift(const std::string &name):
Component(name) {
Component(name),
nSlices_m(1)
{ }
}
Drift::~Drift()
{ }
......@@ -56,6 +59,17 @@ void Drift::initialise(PartBunchBase<double, 3> *bunch, double &startField, doub
startField_m = startField;
}
//set the number of slices for map tracking
void Drift::setNSlices(const std::size_t& nSlices) {
nSlices_m = nSlices;
}
//get the number of slices for map tracking
std::size_t Drift::getNSlices() const {
return nSlices_m;
}
void Drift::finalise() {
}
......@@ -70,4 +84,4 @@ void Drift::getDimensions(double &zBegin, double &zEnd) const {
ElementBase::ElementType Drift::getType() const {
return DRIFT;
}
\ No newline at end of file
}
......@@ -54,12 +54,19 @@ public:
virtual void getDimensions(double &zBegin, double &zEnd) const;
//set number of slices for map tracking
void setNSlices(const std::size_t& nSlices); // Philippe was here
//set number of slices for map tracking
std::size_t getNSlices() const; // Philippe was here
private:
double startField_m;
std::size_t nSlices_m;
// Not implemented.
void operator=(const Drift &);
};
#endif // CLASSIC_Drift_HH
\ No newline at end of file
#endif // CLASSIC_Drift_HH
......@@ -46,7 +46,8 @@ Multipole::Multipole():
SkewComponents(1, 0.0),
SkewComponentErrors(1, 0.0),
max_SkewComponent_m(1),
max_NormalComponent_m(1) {
max_NormalComponent_m(1),
nSlices_m(1) {
setElType(isMultipole);
}
......@@ -58,7 +59,8 @@ Multipole::Multipole(const Multipole &right):
SkewComponents(right.SkewComponents),
SkewComponentErrors(right.SkewComponentErrors),
max_SkewComponent_m(right.max_SkewComponent_m),
max_NormalComponent_m(right.max_NormalComponent_m) {
max_NormalComponent_m(right.max_NormalComponent_m),
nSlices_m(right.nSlices_m) {
setElType(isMultipole);
}
......@@ -70,7 +72,8 @@ Multipole::Multipole(const std::string &name):
SkewComponents(1, 0.0),
SkewComponentErrors(1, 0.0),
max_SkewComponent_m(1),
max_NormalComponent_m(1) {
max_NormalComponent_m(1),
nSlices_m(1) {
setElType(isMultipole);
}
......@@ -161,6 +164,16 @@ void Multipole::setSkewComponent(int n, double v, double vError) {
}
}
//set the number of slices for map tracking
void Multipole::setNSlices(const std::size_t& nSlices) {
nSlices_m = nSlices;
}
//get the number of slices for map tracking
std::size_t Multipole::getNSlices() const {
return nSlices_m;
}
//ff
// radial focussing term
void Multipole::addKR(int i, double t, Vector_t &K) {
......
......@@ -96,6 +96,12 @@ public:
size_t getMaxNormalComponentIndex() const;
size_t getMaxSkewComponentIndex() const;
//set number of slices for map tracking
void setNSlices(const std::size_t& nSlices);
//set number of slices for map tracking
std::size_t getNSlices() const;
bool isFocusing(unsigned int component) const;
......@@ -137,6 +143,7 @@ private:
std::vector<double> SkewComponentErrors;
int max_SkewComponent_m;
int max_NormalComponent_m;
std::size_t nSlices_m;
};
inline
......@@ -159,4 +166,4 @@ size_t Multipole::getMaxSkewComponentIndex() const {
return SkewComponents.size();
}
#endif // CLASSIC_Multipole_HH
\ No newline at end of file
#endif // CLASSIC_Multipole_HH
......@@ -525,6 +525,9 @@ public:
/// timer for selfField calculation
IpplTimings::TimerRef selfFieldTimer_m;
// get 2nd order momentum matrix
FMatrix<double, 2 * Dim, 2 * Dim> getSigmaMatrix();
protected:
IpplTimings::TimerRef boundpTimer_m;
IpplTimings::TimerRef boundpBoundsTimer_m;
......@@ -676,4 +679,4 @@ protected:
#include "PartBunchBase.hpp"
#endif
\ No newline at end of file
#endif
......@@ -2155,7 +2155,7 @@ Inform &PartBunchBase<T, Dim>::print(Inform &os) {
os << "* ************** B U N C H ********************************************************* \n";
os << "* NP = " << getTotalNum() << "\n";
os << "* Qtot = " << std::setw(17) << Util::getChargeString(std::abs(sum(Q))) << " "
<< "Qi = " << std::setw(17) << Util::getChargeString(std::abs(qi_m)) << "\n";
<< "Qi = " << std::setw(17) << Util::getChargeString(std::abs(qi_m)) << "\n";
os << "* Ekin = " << std::setw(17) << Util::getEnergyString(eKin_m) << " "
<< "dEkin = " << std::setw(17) << Util::getEnergyString(dE_m) << "\n";
os << "* rmax = " << Util::getLengthString(rmax_m, 5) << "\n";
......@@ -2560,4 +2560,24 @@ void PartBunchBase<T, Dim>::setBeamFrequency(double f) {
periodLength_m = Physics::c / f;
}
#endif
\ No newline at end of file
template <class T, unsigned Dim>
FMatrix<double, 2 * Dim, 2 * Dim> PartBunchBase<T, Dim>::getSigmaMatrix() {
const double N = static_cast<double>(this->getTotalNum());
Vektor<double, 2*Dim> rpmean;
for (unsigned int i = 0; i < Dim; i++) {
rpmean(2*i)= rmean_m(i);
rpmean((2*i)+1)= pmean_m(i);
}
FMatrix<double, 2 * Dim, 2 * Dim> sigmaMatrix = moments_m / N;
for (unsigned int i = 0; i < 2 * Dim; i++) {
for (unsigned int j = 0; j <= i; j++) {
sigmaMatrix[i][j] = moments_m(i, j) - rpmean(i) * rpmean(j);
sigmaMatrix[j][i] = sigmaMatrix[i][j];
}
}
return sigmaMatrix;
}
#endif
......@@ -97,7 +97,7 @@ protected:
double charge; // Particle charge.
double mass; // Particle mass.
double beta; // particle velocity divided by c.
double gamma; // particle energy divided by particle mas
double gamma; // particle energy divided by particle mass
};
......
......@@ -303,4 +303,4 @@ inline
bool OpalBeamline::containsSource() {
return containsSource_m;
}
#endif // OPAL_BEAMLINE_H
\ No newline at end of file
#endif // OPAL_BEAMLINE_H
......@@ -76,6 +76,9 @@ OpalBend::OpalBend(const char *name, const char *help):
itsAttr[GREATERTHANPI] = Attributes::makeBool
("GREATERTHANPI",
"-- not supported any more --");
itsAttr[NSLICES] = Attributes::makeReal
("NSLICES",
"The number of slices/ steps for this element in Map Tracking", 1);
registerRealAttribute("ANGLE");
registerRealAttribute("K0L");
......@@ -99,6 +102,7 @@ OpalBend::OpalBend(const char *name, const char *help):
registerRealAttribute("HAPERT");
registerRealAttribute("ROTATION");
registerRealAttribute("DESIGNENERGY");
registerRealAttribute("NSLICES");
}
......@@ -115,4 +119,4 @@ void OpalBend::print(std::ostream &os) const {
OpalElement::print(os);
}
\ No newline at end of file
}
......@@ -48,6 +48,7 @@ public:
ROTATION, // Magnet rotation about z axis.
DESIGNENERGY, // the design energy of the particles
GREATERTHANPI, // Boolean flag set to true if bend angle is greater
NSLICES, // The number of slices / steps per element for map tracking
// than 180 degrees.
SIZE // Total number of attributes.
};
......
......@@ -41,8 +41,14 @@ OpalDrift::OpalDrift():
// registerRealAttribute("LENGTH");
itsAttr[GEOMETRY] = Attributes::makeString
("GEOMETRY", "BoundaryGeometry for Drifts");
registerStringAttribute("GEOMETRY");
itsAttr[NSLICES] = Attributes::makeReal
("NSLICES",
"The number of slices/ steps for this element in Map Tracking", 1);
registerStringAttribute("GEOMETRY");
registerRealAttribute("NSLICES");
registerOwnership();
setElement(new DriftRep("DRIFT"));
......@@ -83,6 +89,7 @@ void OpalDrift::update() {
DriftRep *drf = static_cast<DriftRep *>(getElement());
drf->setElementLength(Attributes::getReal(itsAttr[LENGTH]));
drf->setNSlices(Attributes::getReal(itsAttr[NSLICES]));
if(itsAttr[WAKEF] && owk_m == NULL) {
owk_m = (OpalWake::find(Attributes::getString(itsAttr[WAKEF])))->clone(getOpalName() + std::string("_wake"));
owk_m->initWakefunction(*drf);
......@@ -103,4 +110,4 @@ void OpalDrift::update() {
// Transmit "unknown" attributes.
OpalElement::updateUnknown(drf);
}
\ No newline at end of file
}
......@@ -35,6 +35,7 @@ public:
enum {
GEOMETRY = COMMON, // geometry of boundary, one more enum member besides the common ones in OpalElement.
NSLICES, // The number of slices / steps per element for map tracking
SIZE
};
......
......@@ -51,11 +51,16 @@ OpalQuadrupole::OpalQuadrupole():
itsAttr[DK1S] = Attributes::makeReal
("DK1S", "Normalised skew quadrupole coefficient error in m^(-2)");
itsAttr[NSLICES] = Attributes::makeReal
("NSLICES",
"The number of slices/ steps for this element in Map Tracking", 1);
registerRealAttribute("K1");
registerRealAttribute("DK1");
registerRealAttribute("K1S");
registerRealAttribute("DK1S");
registerRealAttribute("NSLICES");
registerOwnership();
......@@ -147,6 +152,7 @@ void OpalQuadrupole::update() {
quad->setField(field);
quad->setNormalComponent(2, Attributes::getReal(itsAttr[K1]), Attributes::getReal(itsAttr[DK1]));
quad->setSkewComponent(2, Attributes::getReal(itsAttr[K1S]), Attributes::getReal(itsAttr[DK1S]));
quad->setNSlices(Attributes::getReal(itsAttr[NSLICES]));
if(itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
parmatint_m = (ParticleMatterInteraction::find(Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION])))->clone(getOpalName() + std::string("_parmatint"));
......@@ -156,4 +162,4 @@ void OpalQuadrupole::update() {
// Transmit "unknown" attributes.
OpalElement::updateUnknown(quad);
}
\ No newline at end of file
}
......@@ -36,6 +36,7 @@ public:
DK1, // The normal quadupole coefficient error.
K1S, // The skew quadrupole coefficient.
DK1S, // The skew quadrupole coefficient error.
NSLICES, // The number of slices / steps per element for map tracking
SIZE
};
......
......@@ -125,6 +125,9 @@ void OpalRBend::update() {
}
geometry.setBendAngle(angle);
// Define number of slices for map tracking
bend->setNSlices(Attributes::getReal(itsAttr[NSLICES]));
// Define pole face angles.
bend->setEntryFaceRotation(Attributes::getReal(itsAttr[E1]));
bend->setExitFaceRotation(Attributes::getReal(itsAttr[E2]));
......@@ -229,4 +232,4 @@ void OpalRBend::update() {
// Transmit "unknown" attributes.
OpalElement::updateUnknown(bend);
}
\ No newline at end of file
}
......@@ -126,6 +126,8 @@ void OpalSBend::update() {
} else {
geometry = PlanarArcGeometry(angle);
}
// Define number of slices for map tracking
bend->setNSlices(Attributes::getReal(itsAttr[NSLICES]));
// Define pole face angles.
bend->setEntryFaceRotation(Attributes::getReal(itsAttr[E1]));
......@@ -236,4 +238,4 @@ void OpalSBend::update() {
// Transmit "unknown" attributes.
OpalElement::updateUnknown(bend);
}
\ No newline at end of file
}
......@@ -57,7 +57,9 @@ Track::Track(BeamSequence *u, const PartData &ref, const std::vector<double> & d
stepsPerTurn(stepsperturn),
zstart(zStart),
zstop(zStop),
timeIntegrator(timeintegrator) {
timeIntegrator(timeintegrator),
truncOrder(1)
{
if(nslices > 0) {
if(!OpalData::getInstance()->hasSLBunchAllocated())
OpalData::getInstance()->setSLPartBunch(new EnvelopeBunch(&ref));
......
......@@ -94,6 +94,9 @@ public:
// 3 --- AMTS
int timeIntegrator;
/// Trunction order for map tracking
int truncOrder;
private:
// Not implemented.
......
......@@ -48,6 +48,7 @@ namespace {
ZSTOP, // Defines a z-location [m], after which the simulation stops when the last particles passes
STEPSPERTURN, // Return the timsteps per revolution period. ONLY available for OPAL-cycl.
TIMEINTEGRATOR, // the name of time integrator
MAP_ORDER, // Truncation order of maps for ThickTracker (default: 1 (linear))
SIZE
};
}
......@@ -77,6 +78,9 @@ TrackCmd::TrackCmd():
("ZSTOP", "Defines a z-location [m], after which the simulation stops when the last particles passes");
itsAttr[TIMEINTEGRATOR] = Attributes::makeString
("TIMEINTEGRATOR", "Name of time integrator to be used", "RK-4");
itsAttr[MAP_ORDER] = Attributes::makeReal
("MAP_ORDER", "Truncation order of maps for ThickTracker (default: 1, i.e. linear)", 1);
registerOwnership(AttributeHandler::COMMAND);
AttributeHandler::addAttributeOwner("TRACK", AttributeHandler::COMMAND, "RUN");
......@@ -199,7 +203,12 @@ void TrackCmd::execute() {
}
// Execute track block.
Track::block = new Track(use, beam->getReference(), dt, maxsteps, stepsperturn, zstart, zstop, timeintegrator, nslices, t0, getDTSCINIT(), getDTAU());
Track::block = new Track(use, beam->getReference(), dt, maxsteps,
stepsperturn, zstart, zstop,
timeintegrator, nslices, t0, getDTSCINIT(), getDTAU());
Track::block->truncOrder = (int)Attributes::getReal(itsAttr[MAP_ORDER]);
Track::block->parser.run();
// Clean up.
......
......@@ -470,9 +470,10 @@ void TrackRun::setupThickTracker()
<< "statistics dump frequency " << Options::statDumpFreq << " w.r.t. the time step." << endl;