Commit 7359476d authored by frey_m's avatar frey_m
Browse files

Matched-Gauss: Fixed overflow bug

parent 17c678c5
...@@ -1306,7 +1306,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1306,7 +1306,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
"'CYCLOTRON' definition."); "'CYCLOTRON' definition.");
} }
int Nint = 2440; int Nint = 1440;
double scaleFactor = 1.0; double scaleFactor = 1.0;
bool writeMap = true; bool writeMap = true;
...@@ -1337,7 +1337,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub ...@@ -1337,7 +1337,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
CyclotronElement->getPHIinit(), CyclotronElement->getPHIinit(),
Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]), Attributes::getReal(itsAttr[Attrib::Distribution::RGUESS]),
Attributes::getString(itsAttr[Attrib::Distribution::FMTYPE]), Attributes::getString(itsAttr[Attrib::Distribution::FMTYPE]),
false)) { false, false)) {
std::array<double,3> Emit = siggen->getEmittances(); std::array<double,3> Emit = siggen->getEmittances();
......
...@@ -155,9 +155,10 @@ class SigmaGenerator ...@@ -155,9 +155,10 @@ class SigmaGenerator
* @param guess value of radius for closed orbit finder * @param guess value of radius for closed orbit finder
* @param type specifies the magnetic field format (e.g. CARBONCYCL) * @param type specifies the magnetic field format (e.g. CARBONCYCL)
* @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder. * @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.
* @param full match over full turn not just single sector
*/ */
bool match(value_type accuracy, size_type maxit, size_type maxitOrbit, value_type angle, bool match(value_type accuracy, size_type maxit, size_type maxitOrbit, value_type angle,
value_type guess, const std::string& type, bool harmonic); value_type guess, const std::string& type, bool harmonic, bool full);
/// Returns True if an Inverse Matrix was found /// Returns True if an Inverse Matrix was found
/* /*
...@@ -221,6 +222,10 @@ class SigmaGenerator ...@@ -221,6 +222,10 @@ class SigmaGenerator
size_type N_m; size_type N_m;
/// Number of integration steps per sector (--> also: number of maps per sector) /// Number of integration steps per sector (--> also: number of maps per sector)
size_type nStepsPerSector_m; size_type nStepsPerSector_m;
/// Number of integration steps in total
size_type nSteps_m;
/// Error of computation /// Error of computation
value_type error_m; value_type error_m;
...@@ -278,6 +283,7 @@ class SigmaGenerator ...@@ -278,6 +283,7 @@ class SigmaGenerator
* @param maxit is the maximum number of iterations (the program stops if within this value no stationary * @param maxit is the maximum number of iterations (the program stops if within this value no stationary
* distribution was found) * distribution was found)
* @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.(not implemented yet) * @param harmonic is a boolean. If "true" the harmonics are used instead of the closed orbit finder.(not implemented yet)
* @param full match over full turn
* @param permutation is the permutation of the order of the Eigenvalues (default: 0) * @param permutation is the permutation of the order of the Eigenvalues (default: 0)
* *
*/ */
...@@ -388,12 +394,14 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type e ...@@ -388,12 +394,14 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type e
Emin_m(Emin), Emax_m(Emax), nSector_m(nSector), N_m(N), nStepsPerSector_m(N/nSector), Emin_m(Emin), Emax_m(Emax), nSector_m(nSector), N_m(N), nStepsPerSector_m(N/nSector),
error_m(std::numeric_limits<value_type>::max()), error_m(std::numeric_limits<value_type>::max()),
fieldmap_m(fieldmap), truncOrder_m(truncOrder), write_m(write), fieldmap_m(fieldmap), truncOrder_m(truncOrder), write_m(write),
scaleFactor_m(scaleFactor), sigmas_m(nStepsPerSector_m),permutations_m{{0,1,2,3},{1,2,3,0},{0,1,3,2},{1,0,2,3},{0,2,3,1},{0,3,2,1}, scaleFactor_m(scaleFactor), sigmas_m(0),permutations_m{{0,1,2,3},{1,2,3,0},{0,1,3,2},{1,0,2,3},{0,2,3,1},{0,3,2,1},
{1,0,3,2},{1,3,2,0},{2,1,0,3},{2,0,1,3},{2,3,0,1},{2,3,1,0},{3,1,0,2},{3,2,0,1}, {3,2,1,0},{3,0,1,2}} {1,0,3,2},{1,3,2,0},{2,1,0,3},{2,0,1,3},{2,3,0,1},{2,3,1,0},{3,1,0,2},{3,2,0,1}, {3,2,1,0},{3,0,1,2}}
// The permutations up converge to a matched distribution, but the first four give the 4 unique solutions // The permutations up converge to a matched distribution, but the first four give the 4 unique solutions
// Following permutations threw an OPAL error, whose reason I still don't know. // Following permutations threw an OPAL error, whose reason I still don't know.
//{0,2,1,3},{0,3,1,2},{1,2,0,3},{1,3,0,2},{2,1,3,0},{2,0,3,1},{3,1,2,0},{3,0,2,1} //{0,2,1,3},{0,3,1,2},{1,2,0,3},{1,3,0,2},{2,1,3,0},{2,0,3,1},{3,1,2,0},{3,0,2,1}
{ {
nSteps_m = N_m;
// set emittances (initialization like that due to old compiler version) // set emittances (initialization like that due to old compiler version)
// [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad // [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad
emittance_m[0] = ex * Physics::pi; emittance_m[0] = ex * Physics::pi;
...@@ -453,7 +461,7 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type e ...@@ -453,7 +461,7 @@ SigmaGenerator<Value_type, Size_type>::SigmaGenerator(value_type I, value_type e
value_type Kx = kxy / sx; value_type Kx = kxy / sx;
value_type Ky = kxy / sy; value_type Ky = kxy / sy;
value_type Kz = K3 * f / (tmp * sz); value_type Kz = K3 * f / (tmp * sz);
return -0.5 * Kx * x_m * x_m return -0.5 * Kx * x_m * x_m
-0.5 * Ky * y_m * y_m -0.5 * Ky * y_m * y_m
-0.5 * Kz * z_m * z_m * gamma2_m; -0.5 * Kz * z_m * z_m * gamma2_m;
...@@ -467,7 +475,8 @@ template<typename Value_type, typename Size_type> ...@@ -467,7 +475,8 @@ template<typename Value_type, typename Size_type>
value_type angle, value_type angle,
value_type rguess, value_type rguess,
const std::string& type, const std::string& type,
bool harmonic) bool harmonic,
bool full)
{ {
/* compute the equilibrium orbit for energy E /* compute the equilibrium orbit for energy E
* and get the the following properties: * and get the the following properties:
...@@ -477,10 +486,13 @@ template<typename Value_type, typename Size_type> ...@@ -477,10 +486,13 @@ template<typename Value_type, typename Size_type>
*/ */
try{ try{
if ( !full )
nSteps_m = nStepsPerSector_m;
// object for space cyclotron map // object for space cyclotron map
MapGenerator<value_type,size_type,Series,Map,Hamiltonian,SpaceCharge> mapgen(nStepsPerSector_m); MapGenerator<value_type,size_type,Series,Map,Hamiltonian,SpaceCharge> mapgen(nSteps_m);
container_type h(nStepsPerSector_m), r(nStepsPerSector_m), ds(nStepsPerSector_m), fidx(nStepsPerSector_m); container_type h(nSteps_m), r(nSteps_m), ds(nSteps_m), fidx(nSteps_m);
value_type ravg = 0.0; value_type ravg = 0.0;
std::pair<value_type,value_type> tunes; std::pair<value_type,value_type> tunes;
...@@ -521,13 +533,18 @@ template<typename Value_type, typename Size_type> ...@@ -521,13 +533,18 @@ template<typename Value_type, typename Size_type>
value_type deg_step = N_m / 360.0; value_type deg_step = N_m / 360.0;
// compute starting point of computation // compute starting point of computation
size_type start = deg_step * angle; size_type start = deg_step * angle;
// copy properties of the length of one sector (--> nStepsPerSector_m) // copy properties of the length of one sector (--> nSteps_m)
size_type step = 0;
std::copy_n(r_turn.begin()+start,nStepsPerSector_m, r.begin()); for (size_type i = 0; i < nSteps_m; ++i) {
std::copy_n(h_turn.begin()+start,nStepsPerSector_m, h.begin());
std::copy_n(fidx_turn.begin()+start,nStepsPerSector_m, fidx.begin()); step = (start + i) % r_turn.size();
std::copy_n(ds_turn.begin()+start,nStepsPerSector_m, ds.begin());
r[i] = r_turn[step];
h[i] = h_turn[step];
fidx[i] = fidx_turn[step];
ds[i] = ds_turn[step];
}
if(write_m){ if(write_m){
writeOrbitOutput_m(tunes,ravg, cof.getFrequencyError(),r_turn,peo,h_turn,fidx_turn,ds_turn); writeOrbitOutput_m(tunes,ravg, cof.getFrequencyError(),r_turn,peo,h_turn,fidx_turn,ds_turn);
...@@ -544,12 +561,12 @@ template<typename Value_type, typename Size_type> ...@@ -544,12 +561,12 @@ template<typename Value_type, typename Size_type>
} }
else{ else{
// computes the cyclotron maps for each angle and stores them into a vector // computes the cyclotron maps for each angle and stores them into a vector
std::vector<matrix_type> Mcycs(nStepsPerSector_m); std::vector<matrix_type> Mcycs(nSteps_m);
// calculate only for a single sector (a nSector_-th) of the whole cyclotron // calculate only for a single sector (a nSector_-th) of the whole cyclotron
for (size_type i = 0; i < nStepsPerSector_m; ++i) for (size_type i = 0; i < nSteps_m; ++i)
Mcycs[i] = mapgen.generateMap(H_m(h[i],h[i]*h[i]+fidx[i],-fidx[i]),ds[i],truncOrder_m); Mcycs[i] = mapgen.generateMap(H_m(h[i],h[i]*h[i]+fidx[i],-fidx[i]),ds[i],truncOrder_m);
// search for a match distribution for each permutation // search for a match distribution for each permutation
// Set to search for the permutation which converged to a matched distribution upon tracking // Set to search for the permutation which converged to a matched distribution upon tracking
// Note that 4 permutations converge to the same matched distribution // Note that 4 permutations converge to the same matched distribution
...@@ -757,7 +774,8 @@ void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz, value_typ ...@@ -757,7 +774,8 @@ void SigmaGenerator<Value_type, Size_type>::initialize(value_type nuz, value_typ
sigma(5,5) = invAB * (ex / (B * Omega) + ez / (A * omega)) * mega; sigma(5,5) = invAB * (ex / (B * Omega) + ez / (A * omega)) * mega;
// fill in initial guess of the sigma matrix (for each angle the same guess) // fill in initial guess of the sigma matrix (for each angle the same guess)
sigmas_m.resize(nStepsPerSector_m); sigmas_m.resize(nSteps_m);
for (typename std::vector<matrix_type>::iterator it = sigmas_m.begin(); it != sigmas_m.end(); ++it) { for (typename std::vector<matrix_type>::iterator it = sigmas_m.begin(); it != sigmas_m.end(); ++it) {
*it = sigma; *it = sigma;
} }
...@@ -782,12 +800,12 @@ template<typename Value_type, typename Size_type> ...@@ -782,12 +800,12 @@ template<typename Value_type, typename Size_type>
bool harmonic, bool harmonic,
size_type permutation){ size_type permutation){
// object for space charge map // object for space charge map
MapGenerator<value_type,size_type,Series,Map,Hamiltonian,SpaceCharge> mapgen(nStepsPerSector_m); MapGenerator<value_type,size_type,Series,Map,Hamiltonian,SpaceCharge> mapgen(nSteps_m);
value_type const_ds = 0.0; value_type const_ds = 0.0;
try { try {
// stores computed space charge map for each angle // stores computed space charge map for each angle
std::vector<matrix_type> Mscs(nStepsPerSector_m); std::vector<matrix_type> Mscs(nSteps_m);
// (inverse) transformation matrix // (inverse) transformation matrix
complex_matrix_type R(6,6), invR(6,6); complex_matrix_type R(6,6), invR(6,6);
...@@ -805,11 +823,12 @@ template<typename Value_type, typename Size_type> ...@@ -805,11 +823,12 @@ template<typename Value_type, typename Size_type>
initialize(tunes.second,ravg); initialize(tunes.second,ravg);
// calculate only for a single sector (a nSector_-th) of the whole cyclotron // calculate only for a single sector (a nSector_-th) of the whole cyclotron
for (size_type i = 0; i < nStepsPerSector_m; ++i) { for (size_type i = 0; i < nSteps_m; ++i) {
if (!harmonic) { if (!harmonic) {
Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),ds[i],truncOrder_m); Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),ds[i],truncOrder_m);
} else { } else {
Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),const_ds,truncOrder_m); Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),const_ds,truncOrder_m);
} }
} }
...@@ -838,12 +857,12 @@ template<typename Value_type, typename Size_type> ...@@ -838,12 +857,12 @@ template<typename Value_type, typename Size_type>
// computes error // computes error
error_m = L2ErrorNorm(sigmas_m[0],newSigma); error_m = L2ErrorNorm(sigmas_m[0],newSigma);
// write new initial sigma-matrix into vector // write new initial sigma-matrix into vector
sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0]; sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];
// compute new space charge maps // compute new space charge maps
for (size_type i = 0; i < nStepsPerSector_m; ++i) { for (size_type i = 0; i < nSteps_m; ++i) {
if (!harmonic) { if (!harmonic) {
Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),ds[i],truncOrder_m); Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),sigmas_m[i](2,2),sigmas_m[i](4,4)),ds[i],truncOrder_m);
} else { } else {
...@@ -1054,7 +1073,7 @@ void SigmaGenerator<Value_type, Size_type>::updateSigma(const std::vector<matrix ...@@ -1054,7 +1073,7 @@ void SigmaGenerator<Value_type, Size_type>::updateSigma(const std::vector<matrix
} }
// initial sigma is already computed // initial sigma is already computed
for (size_type i = 1; i < nStepsPerSector_m; ++i) { for (size_type i = 1; i < nSteps_m; ++i) {
// transfer matrix for one angle // transfer matrix for one angle
M = boost::numeric::ublas::prod(Mscs[i - 1],Mcycs[i - 1]); M = boost::numeric::ublas::prod(Mscs[i - 1],Mcycs[i - 1]);
// transfer the matrix sigma // transfer the matrix sigma
...@@ -1139,7 +1158,7 @@ template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type ...@@ -1139,7 +1158,7 @@ template<typename Value_type, typename Size_type> void SigmaGenerator<Value_type
writeSigma.open("data/SigmaPerAngleForEnergy"+p+"_"+energy+"MeV.dat",std::ios::app); writeSigma.open("data/SigmaPerAngleForEnergy"+p+"_"+energy+"MeV.dat",std::ios::app);
matrix_type M = boost::numeric::ublas::matrix<value_type>(6,6); matrix_type M = boost::numeric::ublas::matrix<value_type>(6,6);
// initial sigma is already computed // initial sigma is already computed
for (size_type i = 1; i < nStepsPerSector_m; ++i) { for (size_type i = 1; i < nSteps_m; ++i) {
// transfer matrix for one angle // transfer matrix for one angle
M = boost::numeric::ublas::prod(Mscs[i - 1],Mcycs[i - 1]); M = boost::numeric::ublas::prod(Mscs[i - 1],Mcycs[i - 1]);
// transfer the matrix sigma // transfer the matrix sigma
......
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