Commit 71fe312b authored by frey_m's avatar frey_m
Browse files

SigmaGenerator: no template class, split into *.h and *.cpp file

parent f0154eb5
set (_SRCS
Distribution.cpp
LaserProfile.cpp
SigmaGenerator.cpp
)
include_directories (
......
......@@ -40,6 +40,7 @@
#include "Utilities/Options.h"
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
#include "Physics/Physics.h"
#include "AbstractObjects/OpalData.h"
......
......@@ -8,6 +8,7 @@
#include "Distribution/Distribution.h"
#include "Distribution/SigmaGenerator.h"
#include "Distribution/ClosedOrbitFinder.h"
#include "AbsBeamline/SpecificElementVisitor.h"
#include <cmath>
......@@ -212,16 +213,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
......@@ -1278,18 +1269,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]),
......@@ -1375,13 +1366,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) {
......
......@@ -279,7 +279,6 @@ private:
void eraseTOrZDist();
void eraseBGzDist();
// void printSigma(SigmaGenerator<double,unsigned int>::matrix_type& M, Inform& out);
void addDistributions();
void applyEmissionModel(double lowEnergyLimit, double &px, double &py, double &pz, std::vector<double> &additionalRNs);
void applyEmissModelAstra(double &px, double &py, double &pz, std::vector<double> &additionalRNs);
......
//
// Class: SigmaGenerator.h
// The SigmaGenerator class uses the class <b>ClosedOrbitFinder</b> to get the parameters(inverse bending radius, path length
// field index and tunes) to initialize the sigma matrix.
// The main function of this class is <b>match(double, unsigned int)</b>, where it iteratively tries to find a matched
// distribution for given
// emittances, energy and current. The computation stops when the L2-norm is smaller than a user-defined tolerance. \n
// In default mode it prints all space charge maps, cyclotron maps and second moment matrices. The orbit properties, i.e.
// tunes, average radius, orbit radius, inverse bending radius, path length, field index and frequency error, are printed
// as well.
//
// Copyright (c) 2014, 2018, Matthias Frey, Cristopher Cortes, ETH Zürich
// All rights reserved
//
// Implemented as part of the Semester thesis by Matthias Frey
// "Matched Distributions in Cyclotrons"
//
// Some adaptations done as part of the Bachelor thesis by Cristopher Cortes
// "Limitations of a linear transfer map method for finding matched distributions in high intensity cyclotrons"
//
// 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 "SigmaGenerator.h"
#include "AbstractObjects/OpalData.h"
#include "Utilities/Options.h"
#include "Utilities/Options.h"
#include "Utilities/OpalException.h"
#include "Utilities/Util.h"
#include "matrix_vector_operation.h"
#include "ClosedOrbitFinder.h"
#include "MapGenerator.h"
#include <cmath>
#include <fstream>
#include <functional>
#include <iomanip>
#include <iterator>
#include <limits>
#include <list>
#include <numeric>
#include <sstream>
#include <boost/numeric/odeint/stepper/runge_kutta4.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
extern Inform *gmsg;
SigmaGenerator::SigmaGenerator(double I,
double ex,
double ey,
double ez,
double E,
double m,
const Cyclotron* cycl,
unsigned int N,
unsigned int Nsectors,
unsigned int truncOrder,
bool write)
: I_m(I)
, wo_m(cycl->getRfFrequ(0)*1E6/cycl->getCyclHarm()*2.0*Physics::pi)
, E_m(E)
, gamma_m(E/m+1.0)
, gamma2_m(gamma_m*gamma_m)
, nh_m(cycl->getCyclHarm())
, beta_m(std::sqrt(1.0-1.0/gamma2_m))
, m_m(m)
, niterations_m(0)
, converged_m(false)
, Emin_m(cycl->getFMLowE())
, Emax_m(cycl->getFMHighE())
, N_m(N)
, nSectors_m(Nsectors)
, nStepsPerSector_m(N/cycl->getSymmetry())
, nSteps_m(N)
, error_m(std::numeric_limits<double>::max())
, truncOrder_m(truncOrder)
, write_m(write)
, sigmas_m(nStepsPerSector_m)
, rinit_m(0.0)
, prinit_m(0.0)
{
// set emittances (initialization like that due to old compiler version)
// [ex] = [ey] = [ez] = pi*mm*mrad --> [emittance] = mm mrad
emittance_m[0] = ex * Physics::pi;
emittance_m[1] = ey * Physics::pi;
emittance_m[2] = ez * Physics::pi;
// minimum beta*gamma
double minGamma = Emin_m / m_m + 1.0;
double bgam = std::sqrt(minGamma * minGamma - 1.0);
// normalized emittance (--> multiply with beta*gamma)
// [emittance] = mm mrad
emittance_m[0] *= bgam;
emittance_m[1] *= bgam;
emittance_m[2] *= bgam;
// Define the Hamiltonian
Series::setGlobalTruncOrder(truncOrder_m);
// infinitesimal elements
x_m = Series::makeVariable(0);
px_m = Series::makeVariable(1);
y_m = Series::makeVariable(2);
py_m = Series::makeVariable(3);
z_m = Series::makeVariable(4);
delta_m = Series::makeVariable(5);
H_m = [&](double h, double kx, double ky) {
return 0.5*px_m*px_m + 0.5*kx*x_m*x_m - h*x_m*delta_m +
0.5*py_m*py_m + 0.5*ky*y_m*y_m +
0.5*delta_m*delta_m/gamma2_m;
};
Hsc_m = [&](double sigx, double sigy, double sigz) {
// convert m from MeV/c^2 to eV*s^{2}/m^{2}
double m = m_m * 1.0e6 / (Physics::c * Physics::c);
// formula (57)
double lam = Physics::two_pi*Physics::c / (wo_m * nh_m); // wavelength, [lam] = m
double K3 = 3.0 * /* physics::q0 */ 1.0 * I_m * lam / (20.0 * std::sqrt(5.0) * Physics::pi * Physics::epsilon_0 * m *
Physics::c * Physics::c * Physics::c * beta_m * beta_m * gamma_m * gamma2_m); // [K3] = m
double milli = 1.0e-3;
// formula (30), (31)
// [sigma(0,0)] = mm^{2} --> [sx] = [sy] = [sz] = mm
// multiply with 0.001 to get meter --> [sx] = [sy] = [sz] = m
double sx = std::sqrt(std::abs(sigx)) * milli;
double sy = std::sqrt(std::abs(sigy)) * milli;
double sz = std::sqrt(std::abs(sigz)) * milli;
double tmp = sx * sy; // [tmp] = m^{2}
double f = std::sqrt(tmp) / (3.0 * gamma_m * sz); // [f] = 1
double kxy = K3 * std::abs(1.0 - f) / ((sx + sy) * sz); // [kxy] = 1/m
double Kx = kxy / sx;
double Ky = kxy / sy;
double Kz = K3 * f / (tmp * sz);
return -0.5 * Kx * x_m * x_m
-0.5 * Ky * y_m * y_m
-0.5 * Kz * z_m * z_m * gamma2_m;
};
}
bool SigmaGenerator::match(double accuracy,
unsigned int maxit,
unsigned int maxitOrbit,
Cyclotron* cycl,
double denergy,
double rguess,
bool full)
{
/* compute the equilibrium orbit for energy E_
* and get the the following properties:
* - inverse bending radius h
* - step sizes of path ds
* - tune nuz
*/
try {
if ( !full )
nSteps_m = nStepsPerSector_m;
// object for space charge map and cyclotron map
MapGenerator<double,
unsigned int,
Series,
Map,
Hamiltonian,
SpaceCharge> mapgen(nSteps_m);
// compute cyclotron map and space charge map for each angle and store them into a vector
std::vector<matrix_type> Mcycs(nSteps_m), Mscs(nSteps_m);
container_type h(nSteps_m), r(nSteps_m), ds(nSteps_m), fidx(nSteps_m);
ClosedOrbitFinder<double, unsigned int,
boost::numeric::odeint::runge_kutta4<container_type> > cof(m_m, N_m, cycl, false, nSectors_m);
if ( !cof.findOrbit(accuracy, maxitOrbit, E_m, denergy, rguess) ) {
throw OpalException("SigmaGenerator::match()",
"Closed orbit finder didn't converge.");
}
cof.computeOrbitProperties(E_m);
// properties of one turn
std::pair<double,double> tunes = cof.getTunes();
double ravg = cof.getAverageRadius(); // average radius
double angle = cycl->getPHIinit();
container_type h_turn = cof.getInverseBendingRadius(angle);
container_type r_turn = cof.getOrbit(angle);
container_type ds_turn = cof.getPathLength(angle);
container_type fidx_turn = cof.getFieldIndex(angle);
container_type peo = cof.getMomentum(angle);
// write properties to file
if (write_m)
writeOrbitOutput_m(tunes, ravg, cof.getFrequencyError(),
r_turn, peo, h_turn, fidx_turn, ds_turn);
// write to terminal
*gmsg << "* ----------------------------" << endl
<< "* Closed orbit info:" << endl
<< "*" << endl
<< "* average radius: " << ravg << " [m]" << endl
<< "* initial radius: " << r_turn[0] << " [m]" << endl
<< "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
<< "* frequency error: " << cof.getFrequencyError()*100 <<" [ % ] "<< endl
<< "* horizontal tune: " << tunes.first << endl
<< "* vertical tune: " << tunes.second << endl
<< "* ----------------------------" << endl << endl;
// copy properties
std::copy_n(r_turn.begin(), nSteps_m, r.begin());
std::copy_n(h_turn.begin(), nSteps_m, h.begin());
std::copy_n(fidx_turn.begin(), nSteps_m, fidx.begin());
std::copy_n(ds_turn.begin(), nSteps_m, ds.begin());
rinit_m = r[0];
prinit_m = peo[0];
std::string fpath = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"maps"
});
if (!boost::filesystem::exists(fpath)) {
boost::filesystem::create_directory(fpath);
}
// initialize sigma matrices (for each angle one) (first guess)
initialize(tunes.second,ravg);
// for writing
std::ofstream writeMturn, writeMcyc, writeMsc;
if (write_m) {
std::string energy = float2string(E_m);
std::string fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"maps",
"OneTurnMapsForEnergy" + energy + "MeV.dat"
});
writeMturn.open(fname, std::ios::app);
fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"maps",
"SpaceChargeMapPerAngleForEnergy" + energy + "MeV_iteration_0.dat"
});
writeMsc.open(fname, std::ios::app);
fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"maps",
"CyclotronMapPerAngleForEnergy" + energy + "MeV.dat"
});
writeMcyc.open(fname, std::ios::app);
}
// calculate only for a single sector (a nSector_-th) of the whole cyclotron
for (unsigned int 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);
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);
writeMatrix(writeMcyc, Mcycs[i]);
writeMatrix(writeMsc, Mscs[i]);
}
if (write_m)
writeMsc.close();
// one turn matrix
mapgen.combine(Mscs,Mcycs);
matrix_type Mturn = mapgen.getMap();
writeMatrix(writeMturn, Mturn);
// (inverse) transformation matrix
sparse_matrix_type R(6, 6), invR(6, 6);
// new initial sigma matrix
matrix_type newSigma(6,6);
// for exiting loop
bool stop = false;
double weight = 0.5;
while (error_m > accuracy && !stop) {
// decouple transfer matrix and compute (inverse) tranformation matrix
decouple(Mturn,R,invR);
// construct new initial sigma-matrix
newSigma = updateInitialSigma(Mturn, R, invR);
// compute new sigma matrices for all angles (except for initial sigma)
updateSigma(Mscs,Mcycs);
// compute error
error_m = L2ErrorNorm(sigmas_m[0],newSigma);
// write new initial sigma-matrix into vector
sigmas_m[0] = weight*newSigma + (1.0-weight)*sigmas_m[0];
// compute new space charge maps
for (unsigned int i = 0; i < nSteps_m; ++i) {
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);
if (write_m) {
std::string energy = float2string(E_m);
std::string fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"maps",
"SpaceChargeMapPerAngleForEnergy" + energy + "MeV_iteration_"
+ std::to_string(niterations_m + 1)
+ ".dat"
});
writeMsc.open(fname, std::ios::out);
}
writeMatrix(writeMsc, Mscs[i]);
if (write_m) {
writeMsc.close();
}
}
// construct new one turn transfer matrix M
mapgen.combine(Mscs,Mcycs);
Mturn = mapgen.getMap();
writeMatrix(writeMturn, Mturn);
// check if number of iterations has maxit exceeded.
stop = (niterations_m++ > maxit);
}
// store converged sigma-matrix
sigma_m.resize(6,6,false);
sigma_m.swap(newSigma);
// returns if the sigma matrix has converged
converged_m = error_m < accuracy;
// Close files
if (write_m) {
writeMturn.close();
writeMcyc.close();
}
} catch(const std::exception& e) {
std::cerr << e.what() << std::endl;
}
if ( converged_m && write_m ) {
// write tunes
std::string fname = Util::combineFilePath({
OpalData::getInstance()->getAuxiliaryOutputDirectory(),
"MatchedDistributions.dat"
});
std::ofstream writeSigmaMatched(fname, std::ios::app);
std::array<double,3> emit = this->getEmittances();
writeSigmaMatched << "* Converged (Ex, Ey, Ez) = (" << emit[0]
<< ", " << emit[1] << ", " << emit[2]
<< ") pi mm mrad for E= " << E_m << " (MeV)"
<< std::endl << "* Sigma-Matrix " << std::endl;
for(unsigned int i = 0; i < sigma_m.size1(); ++ i) {
writeSigmaMatched << std::setprecision(4) << std::setw(11)
<< sigma_m(i,0);
for(unsigned int j = 1; j < sigma_m.size2(); ++ j) {
writeSigmaMatched << " & " << std::setprecision(4)
<< std::setw(11) << sigma_m(i,j);
}
writeSigmaMatched << " \\\\" << std::endl;
}
writeSigmaMatched << std::endl;
writeSigmaMatched.close();
}
return converged_m;
}
void SigmaGenerator::eigsolve_m(const matrix_type& Mturn,
sparse_matrix_type& R)
{
typedef gsl_matrix* gsl_matrix_t;
typedef gsl_vector_complex* gsl_cvector_t;
typedef gsl_matrix_complex* gsl_cmatrix_t;
typedef gsl_eigen_nonsymmv_workspace* gsl_wspace_t;
typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
gsl_cvector_t evals = gsl_vector_complex_alloc(6);
gsl_cmatrix_t evecs = gsl_matrix_complex_alloc(6, 6);
gsl_wspace_t wspace = gsl_eigen_nonsymmv_alloc(6);
gsl_matrix_t M = gsl_matrix_alloc(6, 6);
// go to GSL
for (unsigned int i = 0; i < 6; ++i){
for (unsigned int j = 0; j < 6; ++j) {
gsl_matrix_set(M, i, j, Mturn(i,j));
}
}
/*int status = */gsl_eigen_nonsymmv(M, evals, evecs, wspace);
// if ( !status )
// throw OpalException("SigmaGenerator::eigsolve_m()",
// "Couldn't perform eigendecomposition!");
/*status = *///gsl_eigen_nonsymmv_sort(evals, evecs, GSL_EIGEN_SORT_ABS_ASC);
// if ( !status )
// throw OpalException("SigmaGenerator::eigsolve_m()",
// "Couldn't sort eigenvalues and eigenvectors!");
// go to UBLAS
for( unsigned int i = 0; i < 6; i++){
gsl_vector_complex_view evec_i = gsl_matrix_complex_column(evecs, i);
for(unsigned int j = 0;j < 6; j++){
gsl_complex zgsl = gsl_vector_complex_get(&evec_i.vector, j);
complex_t z(GSL_REAL(zgsl), GSL_IMAG(zgsl));
R(i,j) = z;
}
}
// Sorting the Eigenvectors
// This is an arbitrary threshold that has worked for me. (We should fix this)
double threshold = 10e-12;
bool isZdirection = false;
std::vector<complex_vector_type> zVectors{};
std::vector<complex_vector_type> xyVectors{};
for(unsigned int i = 0; i < 6; i++){
complex_t z = R(i,0);
if(std::abs(z) < threshold) z = 0.;
if(z == 0.) isZdirection = true;
complex_vector_type v(6);
if(isZdirection){
for(unsigned int j = 0;j < 6; j++){
complex_t z = R(i,j);
v(j) = z;
}
zVectors.push_back(v);
}
else{
for(unsigned int j = 0; j < 6; j++){
complex_t z = R(i,j);
v(j) = z;
}
xyVectors.push_back(v);
}
isZdirection = false;
}
//if z-direction not found, then the system does not behave as expected
if(zVectors.size() != 2)
throw OpalException("SigmaGenerator::eigsolve_m()",
"Couldn't find the vertical Eigenvectors.");