Commit 5b3a63e3 authored by cortes_c's avatar cortes_c
Browse files

Fixed a bug in Sigma generator. In class Distribution, multiple simulations...

Fixed a bug in Sigma generator. In class Distribution, multiple simulations are allowed without throwing an mpi - abort, when Option cloTuneOnly is active. ClosedOrbitFinder shouldn't have significant changes.
parent f1aa535b
......@@ -511,7 +511,7 @@ bool ClosedOrbitFinder<Value_type, Size_type, Stepper>::findOrbit(value_type acc
dydt[0] = y[0] * y[1] * invptheta;
// Gordon, formula (5b)
dydt[1] = ptheta - y[0] * bint;
// Gordon, formulas (9a) and (9b)
// Gordon, formulas (9a) and (9b)
for (size_type i = 2; i < 5; i += 2) {
dydt[i] = (y[1] * y[i] + y[0] * p2 * y[i+1] * invptheta * invptheta) * invptheta;
dydt[i+1] = - y[1] * y[i+1] * invptheta - (bint + y[0] * brint) * y[i];
......
......@@ -1366,8 +1366,8 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
*/
if (Options::cloTuneOnly)
throw OpalException("Do only CLO and tune calculation","... ");
//if (Options::cloTuneOnly)
// throw OpalException("Do only CLO and tune calculation","... ");
auto sigma = siggen->getSigma();
......
......@@ -471,6 +471,7 @@ template<typename Value_type, typename Size_type>
* - tune nuz
*/
// object for space cyclotron map
maxit = 1000;
try{
MapGenerator<value_type,size_type,Series,Map,Hamiltonian,SpaceCharge> mapgen(nStepsPerSector_m);
......@@ -531,25 +532,27 @@ template<typename Value_type, typename Size_type>
return false;
}
if(Options::cloTuneOnly)
throw OpalException("Do only CLO and tune calculation","... ");
// computes the cyclotron maps for each angle and stores them into a vector
std::vector<matrix_type> Mcycs(nStepsPerSector_m);
// calculate only for a single sector (a nSector_-th) of the whole cyclotron
for (size_type i = 0; i < nStepsPerSector_m; ++i)
Mcycs[i] = mapgen.generateMap(H_m(h[i],h[i]*h[i]+fidx[i],-fidx[i]),ds[i],truncOrder_m);
if(Options::cloTuneOnly){
//throw OpalException("Do only CLO and tune calculation","... ");
*gmsg << "Doing just close orbit and tune calculation only" << endl;
}
else{
// computes the cyclotron maps for each angle and stores them into a vector
std::vector<matrix_type> Mcycs(nStepsPerSector_m);
// calculate only for a single sector (a nSector_-th) of the whole cyclotron
for (size_type i = 0; i < nStepsPerSector_m; ++i)
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
// Set to search for the permutation which converged to a matched distribution upon tracking
// Note that 4 permutations converge to the same matched distribution
for(size_type permutation = 0; permutation < 1; permutation++){
findMatchDistribution(accuracy,h,fidx,ds,tunes,ravg,maxit,Mcycs,harmonic,2);
// search for a match distribution for each permutation
// Set to search for the permutation which converged to a matched distribution upon tracking
// Note that 4 permutations converge to the same matched distribution
for(size_type permutation = 0; permutation < 1; permutation++){
findMatchDistribution(accuracy,h,fidx,ds,tunes,ravg,maxit,Mcycs,harmonic,2);
}
std::cout << "OPAL> * Number of matched distributions: " <<matchedSigmas_m.size() << std::endl;
std::cout << "OPAL> * Last matched distribution information: "<<std::endl;
}
std::cout << "OPAL> * Number of matched distributions: " <<matchedSigmas_m.size() << std::endl;
std::cout << "OPAL> * Last matched distribution information: "<<std::endl;
}catch(const std::exception& e) {
std::cerr << e.what() << std::endl;
}
......@@ -557,7 +560,7 @@ template<typename Value_type, typename Size_type>
if ( 0 < matchedSigmas_m.size()) return (0 < matchedSigmas_m.size());
//Returns an identity matrix, if no matched distribution was found
else{
*gmsg << "No matched distribution found, returning unity matrix instead.";
*gmsg << "No matched distribution found, returning unity matrix instead or CLOTUNEONLY option"<<endl;
matrix_type One(6,6);
for(size_type i = 0; i < 6; i++){
for(size_type j = 0; j < 6; j++) One(i,j) = 0;
......@@ -789,6 +792,8 @@ template<typename Value_type, typename Size_type>
// for exiting the loop if maxit is reached
bool stop = false;
// checks if the eigenvector matrix E and E^-1 could be built
bool isBuild = false;
// initialize sigma matrices (for each angle one) (first guess)
initialize(tunes.second,ravg);
......@@ -812,10 +817,11 @@ template<typename Value_type, typename Size_type>
* and updateInitialSigma was changed according to Eq. 18 Sec. 2.4
* Frey Semester thesis
*/
// Returns the Eigenvalues of the Twiss-Matrix,
// Sets the Eigenvalues of the Twiss-Matrix,
// orders them and computes R and invR
setEigenVectors(Mturn,R,invR,permutation);
// Also checks if the Eigenvector matrix E is invertible
isBuild = setEigenVectors(Mturn,R,invR,permutation);
// computes the new sigma matrix
newSigma = updateInitialSigma(R,invR);
......@@ -840,9 +846,9 @@ template<typename Value_type, typename Size_type>
// construct new one turn transfer matrix M
mapgen.combine(Mscs,Mcycs);
Mturn = mapgen.getMap();
// check if number of iterations has maxit exceeded.
stop = (niterations_m++ > maxit);
// check if number of iterations has maxit exceeded or if the eigenvector matrix could be built
stop = (niterations_m++ > maxit) || !isBuild;
// check for convergence
converged_m = error_m < accuracy;
......@@ -915,20 +921,17 @@ template<typename Value_type, typename Size_type>
gsl_matrix_complex_free(eigenVectors);
gsl_matrix_free(MTurngsl);
// Sorts the Eigenvectors such that
// Re(eig_x) < Re(eig_y) < Re(eig_z)
// Sorts the Eigenvectors
std::vector<size_type> permutation = permutations_m[perm];
value_type threshold = 10e-12;
std::vector<complex_vector_type> zVectors{};
std::vector<complex_vector_type> xyVectors{};
bool isZdirection = false;
std::vector<complex_vector_type> zVectors;
std::vector<complex_vector_type> xyVectors;
for(size_type i = 0; i < 6; i++){
complex_number_type z = R(i,0);
if(std::abs(z) < 1e-13) z = 0.;
if(std::abs(z) < threshold) z = 0.;
if(z == 0.) isZdirection = true;
complex_vector_type v(6);
if(isZdirection){
for(size_type j = 0;j < matrixDimension; j++){
......@@ -944,7 +947,17 @@ template<typename Value_type, typename Size_type>
}
xyVectors.push_back(v);
}
isZdirection = false;
isZdirection = false;
}
//z-direction not found, then the system does not have the expected form
if(zVectors.size() != 2){
*gmsg << "Could not found the vertical eigenvectors of the system"<<endl;
complex_matrix_type Copy = R;
for(size_type i = 0; i < 6; i++){
for(size_type j = 0; j < 6; j++) R(i,j) = Copy(j,i);
}
InvertMatrix(R,invR);
return false;
}
// Norms the Eigenvectors
for(size_type i = 0; i < 4; i++){
......
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