Commit a41bd0b0 authored by kraus's avatar kraus
Browse files

Merge remote-tracking branch 'origin/master' into...

Merge remote-tracking branch 'origin/master' into 509-orbitthreader-computemaximalimplicitdrift-not-yielding-correct-results
parents 3c7ea342 05397703
......@@ -71,7 +71,7 @@ double CavityAutophaser::getPhaseAtMaxEnergy(const Vector_t &R,
<< std::left << std::setw(68) << std::setfill('*') << ss.str()
<< std::setfill(' ') << endl);
if (!apVeto) {
double initialEnergy = Util::getEnergy(P, itsReference_m.getM()) * 1e-6;
double initialEnergy = Util::getKineticEnergy(P, itsReference_m.getM()) * 1e-6;
double AstraPhase = 0.0;
double designEnergy = element->getDesignEnergy();
......@@ -190,7 +190,7 @@ double CavityAutophaser::guessCavityPhase(double t) {
return orig_phi;
}
Phimax = element->getAutoPhaseEstimate(getEnergyMeV(refP),
Phimax = element->getAutoPhaseEstimate(Util::getKineticEnergy(refP, itsReference_m.getM()) * 1e-6,
t,
itsReference_m.getQ(),
itsReference_m.getM() * 1e-6);
......@@ -287,7 +287,7 @@ double CavityAutophaser::track(Vector_t /*R*/,
out);
rfc->setPhasem(initialPhase);
double finalKineticEnergy = Util::getEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
double finalKineticEnergy = Util::getKineticEnergy(Vector_t(0.0, 0.0, pe.first), itsReference_m.getM() * 1e-6);
return finalKineticEnergy;
}
\ No newline at end of file
......@@ -28,8 +28,6 @@ private:
const double dt,
const double phase,
std::ofstream *out = NULL) const;
double getEnergyMeV(const Vector_t &P);
double getMomentum(double kineticEnergyMeV);
const PartData &itsReference_m;
std::shared_ptr<Component> itsCavity_m;
......@@ -39,14 +37,4 @@ private:
};
inline
double CavityAutophaser::getEnergyMeV(const Vector_t &P) {
return itsReference_m.getM() * 1e-6 * (sqrt(dot(P,P) + 1) - 1);
}
inline
double CavityAutophaser::getMomentum(double kineticEnergyMeV) {
return sqrt(std::pow(kineticEnergyMeV / (itsReference_m.getM() * 1e-6) + 1, 2) - 1);
}
#endif
\ No newline at end of file
......@@ -1160,7 +1160,7 @@ void ParallelTTracker::findStartPosition(const BorisPusher &pusher) {
selectDT();
if ((back_track && itsOpalBeamline_m.containsSource()) ||
Util::getEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
Util::getKineticEnergy(itsBunch_m->RefPartP_m, itsBunch_m->getM()) < 1e-3) {
double gamma = 0.1 / itsBunch_m->getM() + 1.0;
double beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
itsBunch_m->RefPartP_m = itsBunch_m->toLabTrafo_m.rotateTo(beta * gamma * Vector_t(0, 0, 1));
......
......@@ -494,22 +494,22 @@ ElementBase::ElementType RFCavity::getType() const {
double RFCavity::getAutoPhaseEstimateFallback(double E0, double t0, double q, double mass) {
const double dt = 1e-13;
const double p0 = Util::getP(E0, mass);
const double p0 = Util::getBetaGamma(E0, mass);
const double origPhase =getPhasem();
double dphi = Physics::pi / 18;
double phi = 0.0;
setPhasem(phi);
std::pair<double, double> ret = trackOnAxisParticle(E0 / mass, t0, dt, q, mass);
std::pair<double, double> ret = trackOnAxisParticle(p0, t0, dt, q, mass);
double phimax = 0.0;
double Emax = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
double Emax = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
phi += dphi;
for (unsigned int j = 0; j < 2; ++ j) {
for (unsigned int i = 0; i < 36; ++ i, phi += dphi) {
setPhasem(phi);
ret = trackOnAxisParticle(p0, t0, dt, q, mass);
double Ekin = Util::getEnergy(Vector_t(0.0, 0.0, ret.first), mass);
double Ekin = Util::getKineticEnergy(Vector_t(0.0, 0.0, ret.first), mass);
if (Ekin > Emax) {
Emax = Ekin;
phimax = phi;
......@@ -635,7 +635,7 @@ double RFCavity::getAutoPhaseEstimate(const double &E0, const double &t0, const
}
double cosine_part = 0.0, sine_part = 0.0;
double p0 = std::sqrt((E0 / mass + 1) * (E0 / mass + 1) - 1);
double p0 = Util::getBetaGamma(E0, mass);
cosine_part += scale_m * std::cos(frequency_m * t0) * F[0];
sine_part += scale_m * std::sin(frequency_m * t0) * F[0];
......@@ -675,11 +675,11 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
const double zend = length_m + startField_m;
Vector_t z(0.0, 0.0, zbegin);
double dz = 0.5 * p(2) / std::sqrt(1.0 + dot(p, p)) * cdt;
double dz = 0.5 * p(2) / Util::getGamma(p) * cdt;
Vector_t Ef(0.0), Bf(0.0);
if (out) *out << std::setw(18) << z[2]
<< std::setw(18) << Util::getEnergy(p, mass)
<< std::setw(18) << Util::getKineticEnergy(p, mass)
<< std::endl;
while (z(2) + dz < zend && z(2) + dz > zbegin) {
z /= cdt;
......@@ -700,7 +700,7 @@ std::pair<double, double> RFCavity::trackOnAxisParticle(const double &p0,
t += dt;
if (out) *out << std::setw(18) << z[2]
<< std::setw(18) << Util::getEnergy(p, mass)
<< std::setw(18) << Util::getKineticEnergy(p, mass)
<< std::endl;
}
......
......@@ -123,8 +123,7 @@ void PartBunch::computeSelfFields(int binNumber) {
if(fs_m->hasValidSolver()) {
/// Mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// Scatter charge onto space charge grid.
this->Q *= this->dt;
......@@ -328,12 +327,14 @@ void PartBunch::computeSelfFields(int binNumber) {
}
void PartBunch::resizeMesh() {
if (fs_m->getFieldSolverType() != "SAAMG") {
return;
}
double xmin = fs_m->solver_m->getXRangeMin();
double xmax = fs_m->solver_m->getXRangeMax();
double ymin = fs_m->solver_m->getYRangeMin();
double ymax = fs_m->solver_m->getYRangeMax();
double zmin = fs_m->solver_m->getZRangeMin();
double zmax = fs_m->solver_m->getZRangeMax();
if(xmin > rmin_m[0] || xmax < rmax_m[0] ||
ymin > rmin_m[1] || ymax < rmax_m[1]) {
......@@ -354,14 +355,14 @@ void PartBunch::resizeMesh() {
boundp();
get_bounds(rmin_m, rmax_m);
}
Vector_t mymin = Vector_t(xmin, ymin , zmin);
Vector_t mymax = Vector_t(xmax, ymax , zmax);
for (int i = 0; i < 3; i++)
hr_m[i] = (mymax[i] - mymin[i])/nr_m[i];
Vector_t origin = Vector_t(0.0, 0.0, 0.0);
// update the mesh origin and mesh spacing hr_m
fs_m->solver_m->resizeMesh(origin, hr_m, rmin_m, rmax_m, dh_m);
getMesh().set_meshSpacing(&(hr_m[0]));
getMesh().set_origin(mymin);
getMesh().set_origin(origin);
rho_m.initialize(getMesh(),
getFieldLayout(),
......@@ -384,8 +385,7 @@ void PartBunch::computeSelfFields() {
if(fs_m->hasValidSolver()) {
//mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
//scatter charges onto grid
this->Q *= this->dt;
......@@ -510,8 +510,7 @@ void PartBunch::computeSelfFields_cycl(double gamma) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......@@ -639,8 +638,7 @@ void PartBunch::computeSelfFields_cycl(int bin) {
if(fs_m->hasValidSolver()) {
/// mesh the whole domain
if(fs_m->getFieldSolverType() == "SAAMG")
resizeMesh();
resizeMesh();
/// scatter particles charge onto grid.
this->Q.scatter(this->rho_m, this->R, IntrplCIC_t());
......
......@@ -411,6 +411,8 @@ public:
// virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
virtual FieldLayout_t &getFieldLayout() = 0;
virtual void resizeMesh() { };
/*
* Wrapped member functions of IpplParticleBase
*/
......
......@@ -6,6 +6,7 @@
#include <string>
#include <cstring>
#include <limits>
#include <sstream>
#include <type_traits>
#include <functional>
......@@ -28,16 +29,23 @@ namespace Util {
}
inline
double getEnergy(Vector_t p, double mass) {
double getKineticEnergy(Vector_t p, double mass) {
return (getGamma(p) - 1.0) * mass;
}
inline
double getP(double E, double mass) {
double gamma = E / mass + 1;
return std::sqrt(std::pow(gamma, 2.0) - 1.0);
double getBetaGamma(double Ekin, double mass) {
double value = std::sqrt(std::pow(Ekin / mass + 1.0, 2) - 1.0);
if (value < std::numeric_limits<double>::epsilon())
value = std::sqrt(2 * Ekin / mass);
return value;
}
inline
double convertMomentumeVToBetaGamma(double p, double mass) {
return p / mass;
}
inline
std::string getTimeString(double time, unsigned int precision = 3) {
std::string timeUnit(" [ps]");
......@@ -162,7 +170,7 @@ namespace Util {
}
Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName = "");
std::string toUpper(const std::string &str);
std::string combineFilePath(std::initializer_list<std::string>);
......
set (_SRCS
Distribution.cpp
LaserProfile.cpp
RealDiracMatrix.cpp
SigmaGenerator.cpp
)
include_directories (
......@@ -12,10 +14,10 @@ add_opal_sources(${_SRCS})
set (HDRS
ClosedOrbitFinder.h
Distribution.h
Harmonics.h
LaserProfile.h
MapGenerator.h
matrix_vector_operation.h
RealDiracMatrix.h
SigmaGenerator.h
)
......
......@@ -40,6 +40,7 @@
#include "Utilities/Options.h"
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
#include "Physics/Physics.h"
#include "AbstractObjects/OpalData.h"
......
......@@ -7,7 +7,7 @@
// OPAL is licensed under GNU GPL version 3.
#include "Distribution/Distribution.h"
#include "Distribution/SigmaGenerator.h"
#include "Distribution/ClosedOrbitFinder.h"
#include "AbsBeamline/SpecificElementVisitor.h"
#include <cmath>
......@@ -17,7 +17,6 @@
#include <string>
#include <vector>
#include <numeric>
#include <limits>
// IPPL
#include "DataSource/DataConnect.h"
......@@ -212,16 +211,6 @@ Distribution::~Distribution() {
delete laserProfile_m;
}
/*
void Distribution::printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out) {
for (int i=0; i<M.size1(); ++i) {
for (int j=0; j<M.size2(); ++j) {
*gmsg << M(i,j) << " ";
}
*gmsg << endl;
}
}
*/
/**
* Calculate the local number of particles evenly and adjust node 0
......@@ -885,14 +874,6 @@ void Distribution::chooseInputMomentumUnits(InputMomentumUnitsT::InputMomentumUn
}
double Distribution::converteVToBetaGamma(double valueIneV, double massIneV) {
double value = std::copysign(std::sqrt(std::pow(std::abs(valueIneV) / massIneV + 1.0, 2) - 1.0), valueIneV);
if (std::abs(value) < std::numeric_limits<double>::epsilon())
value = std::copysign(std::sqrt(2 * std::abs(valueIneV) / massIneV), valueIneV);
return value;
}
void Distribution::createDistributionBinomial(size_t numberOfParticles, double massIneV) {
setDistParametersBinomial(massIneV);
......@@ -1079,9 +1060,9 @@ void Distribution::createDistributionFromFile(size_t /*numberOfParticles*/, doub
saveProcessor = 0;
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
P(0) = converteVToBetaGamma(P(0), massIneV);
P(1) = converteVToBetaGamma(P(1), massIneV);
P(2) = converteVToBetaGamma(P(2), massIneV);
P(0) = Util::convertMomentumeVToBetaGamma(P(0), massIneV);
P(1) = Util::convertMomentumeVToBetaGamma(P(1), massIneV);
P(2) = Util::convertMomentumeVToBetaGamma(P(2), massIneV);
}
pmean_m += P;
......@@ -1278,18 +1259,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
bool writeMap = true;
typedef SigmaGenerator<double, unsigned int> sGenerator_t;
sGenerator_t *siggen = new sGenerator_t(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
CyclotronElement,
Nint,
Nsectors,
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap);
std::unique_ptr<SigmaGenerator> siggen = std::unique_ptr<SigmaGenerator>(
new SigmaGenerator(I_m,
Attributes::getReal(itsAttr[Attrib::Distribution::EX])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::EY])*1E6,
Attributes::getReal(itsAttr[Attrib::Distribution::ET])*1E6,
E_m*1E-6,
massIneV*1E-6,
CyclotronElement,
Nint,
Nsectors,
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
writeMap));
if (siggen->match(accuracy,
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
......@@ -1297,7 +1278,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
CyclotronElement,
denergy,
rguess,
false, full)) {
full)) {
std::array<double,3> Emit = siggen->getEmittances();
......@@ -1311,62 +1292,18 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
*gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
if (siggen->getSigma()(i,j) < 10e-12){
if (std::abs(siggen->getSigma()(i,j)) < 1.0e-15) {
*gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
}
else{
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
}
*gmsg << " \\\\" << endl;
}
/*
Now setup the distribution generator
Units of the Sigma Matrix:
spatial: mm
moment: rad
*/
double gamma = E_m / massIneV + 1.0;
double beta = std::sqrt(1.0 - 1.0 / (gamma * gamma));
auto sigma = siggen->getSigma();
// change units from mm to m
for (unsigned int i = 0; i < 6; ++ i)
for (unsigned int j = 0; j < 6; ++ j) sigma(i, j) *= 1e-6;
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::CreateMatchedGaussDistribution()",
"Negative value on the diagonal of the sigma matrix.");
}
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1))*beta*gamma;
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3))*beta*gamma;
sigmaR_m[1] = std::sqrt(sigma(4, 4));
//p_l^2 = [(delta+1)*beta*gamma]^2 - px^2 - pz^2
double pl2 = (std::sqrt(sigma(5,5)) + 1)*(std::sqrt(sigma(5,5)) + 1)*beta*gamma*beta*gamma -
sigmaP_m[0]*sigmaP_m[0] - sigmaP_m[2]*sigmaP_m[2];
double pl = std::sqrt(pl2);
sigmaP_m[1] = gamma*(pl - beta*gamma);
correlationMatrix_m(1, 0) = sigma(0, 1) / (std::sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (std::sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (std::sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (std::sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (std::sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (std::sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (std::sqrt(sigma(1, 1) * sigma(5, 5)));
createDistributionGauss(numberOfParticles, massIneV);
generateMatchedGauss(siggen->getSigma(), numberOfParticles, massIneV);
// update injection radius and radial momentum
CyclotronElement->setRinit(siggen->getInjectionRadius() * 1.0e3);
......@@ -1375,13 +1312,9 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
else {
*gmsg << "* Not converged for " << E_m*1E-6 << " MeV" << endl;
delete siggen;
throw OpalException("Distribution::CreateMatchedGaussDistribution",
"didn't find any matched distribution.");
}
delete siggen;
}
void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV) {
......@@ -1480,7 +1413,7 @@ void Distribution::createOpalT(PartBunchBase<double, 3> *beam,
// This is PC from BEAM
double deltaP = Attributes::getReal(itsAttr[Attrib::Distribution::OFFSETP]);
if (inputMoUnits_m == InputMomentumUnitsT::EV) {
deltaP = converteVToBetaGamma(deltaP, beam->getM());
deltaP = Util::convertMomentumeVToBetaGamma(deltaP, beam->getM());
}
avrgpz_m = beam->getP()/beam->getM() + deltaP;
......@@ -2421,6 +2354,136 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
gsl_matrix_free(corMat);
}
void Distribution::generateMatchedGauss(const SigmaGenerator::matrix_t& sigma,
size_t numberOfParticles, double massIneV)
{
/* This particle distribution generation is based on a
* symplectic method described in
* https://arxiv.org/abs/1205.3601
*/
/* Units of the Sigma Matrix:
* spatial: m
* moment: rad
*
* Attention: The vertical and longitudinal direction must be interchanged!
*/
for (unsigned int i = 0; i < 3; ++ i) {
if ( sigma(2 * i, 2 * i) < 0 || sigma(2 * i + 1, 2 * i + 1) < 0 )
throw OpalException("Distribution::generateMatchedGauss()",
"Negative value on the diagonal of the sigma matrix.");
}
double bgam = Util::getBetaGamma(E_m, massIneV);
/*
* only used for printing
*/
// horizontal
sigmaR_m[0] = std::sqrt(sigma(0, 0));
sigmaP_m[0] = std::sqrt(sigma(1, 1)) * bgam;
// longitudinal
sigmaR_m[1] = std::sqrt(sigma(4, 4));
sigmaP_m[1] = std::sqrt(sigma(5, 5)) * bgam;
// vertical
sigmaR_m[2] = std::sqrt(sigma(2, 2));
sigmaP_m[2] = std::sqrt(sigma(3, 3)) * bgam;
correlationMatrix_m(1, 0) = sigma(0, 1) / (std::sqrt(sigma(0, 0) * sigma(1, 1)));
correlationMatrix_m(3, 2) = sigma(2, 3) / (std::sqrt(sigma(2, 2) * sigma(3, 3)));
correlationMatrix_m(5, 4) = sigma(4, 5) / (std::sqrt(sigma(4, 4) * sigma(5, 5)));
correlationMatrix_m(4, 0) = sigma(0, 4) / (std::sqrt(sigma(0, 0) * sigma(4, 4)));
correlationMatrix_m(4, 1) = sigma(1, 4) / (std::sqrt(sigma(1, 1) * sigma(4, 4)));
correlationMatrix_m(5, 0) = sigma(0, 5) / (std::sqrt(sigma(0, 0) * sigma(5, 5)));
correlationMatrix_m(5, 1) = sigma(1, 5) / (std::sqrt(sigma(1, 1) * sigma(5, 5)));
inputMoUnits_m = InputMomentumUnitsT::NONE;
/*
* decouple horizontal and longitudinal direction
*/
// extract horizontal and longitudinal directions
RealDiracMatrix::matrix_t A(4, 4);
A(0, 0) = sigma(0, 0);
A(1, 1) = sigma(1, 1);
A(2, 2) = sigma(4, 4);
A(3, 3) = sigma(5, 5);
A(0, 1) = sigma(0, 1);
A(0, 2) = sigma(0, 4);
A(0, 3) = sigma(0, 5);
A(1, 0) = sigma(1, 0);
A(2, 0) = sigma(4, 0);
A(3, 0) = sigma(5, 0);
A(1, 2) = sigma(1, 4);
A(2, 1) = sigma(4, 1);
A(1, 3) = sigma(1, 5);
A(3, 1) = sigma(5, 1);
A(2, 3) = sigma(4, 5);
A(3, 2) = sigma(5, 4);
RealDiracMatrix rdm;
RealDiracMatrix::sparse_matrix_t R1 = rdm.diagonalize(A);
RealDiracMatrix::vector_t variances(8);
for (int i = 0; i < 4; ++i) {
variances(i) = std::sqrt(A(i, i));
}
/*
* decouple vertical direction
*/
A *= 0.0;
A(0, 0) = 1;
A(1, 1) = 1;
A(2, 2) = sigma(2, 2);
A(3, 3) = sigma(3, 3);
A(2, 3) = sigma(2, 3);
A(3, 2) = sigma(3, 2);
RealDiracMatrix::sparse_matrix_t R2 = rdm.diagonalize(A);
for (int i = 0; i < 4; ++i) {
variances(4 + i) = std::sqrt(A(i, i));
}
int saveProcessor = -1;
const int myNode = Ippl::myNode();
const int numNodes = Ippl::getNodes();
const bool scalable = Attributes::getBool(itsAttr[Attrib::Distribution::SCALABLE]);
RealDiracMatrix::vector_t p1(4), p2(4);
for (size_t i = 0; i < numberOfParticles; i++) {
for (int j = 0; j < 4; j++) {
p1(j) = gsl_ran_gaussian(randGen_m, 1.0) * variances(j);
p2(j) = gsl_ran_gaussian(randGen_m, 1.0) * variances(4 + j);
}
p1 = boost::numeric::ublas::prod(R1, p1);
p2 = boost::numeric::ublas::prod(R2, p2);
// Save to each processor in turn.