Skip to content

Commit

Permalink
Simple Ring modification to solve energy loss computation bug (#710)
Browse files Browse the repository at this point in the history
* Remove U0 and damping elements from SimpleQuantDiffPass

* Add simple radiation pass method

* Modify SimpleQuantDiff Element and add SimpleRadiation element

* Fix tests

* Renamed SimpleRadPass. Modified error in lattice_object passmethod tuples, changed radiation to longdict in element, modified list of disable_6d passmethods in get_energy_loss and modified fast ring

* Update len of simple_ring in test_physics which is now different

* Add damping matrix

* Add radiative element

* Add dispersion and dispersion '

* Fix normalisation order which was problem with dispersion

* remove alpha and beta, change damp_mat to keep only diagonals, add build attributes, add exponent instead of linear

* Fix tests

* remove alpha and beta all over, but keeping dispersion

* bug in mexall

* Add dispersion to linear matrix

* Add optional for disp

---------

Co-authored-by: Lee Carver <[email protected]>
Co-authored-by: Lee Carver <[email protected]>
  • Loading branch information
3 people committed Dec 4, 2023
1 parent 4b4c363 commit d95aa5c
Show file tree
Hide file tree
Showing 7 changed files with 251 additions and 66 deletions.
53 changes: 6 additions & 47 deletions atintegrators/SimpleQuantDiffPass.c
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,11 @@ struct elem
double betay;
double sigma_xp;
double sigma_yp;
double U0;
double EnergyLossFactor;
};

void SimpleQuantDiffPass(double *r_in,
double sigma_xp, double sigma_yp, double espread,
double taux, double tauy, double tauz, double EnergyLossFactor,
double taux, double tauy, double tauz,
pcg32_random_t* rng, int num_particles)

{
Expand All @@ -36,28 +34,15 @@ void SimpleQuantDiffPass(double *r_in,
}

if(!atIsNaN(r6[0])) {
if(taux!=0.0) {
r6[1] -= 2*r6[1]/taux;
}
if(sigma_xp!=0.0) {
r6[1] += 2*sigma_xp*sqrt(1/taux)*randnorm[0];
}
if(sigma_yp!=0.0) {
r6[3] += 2*sigma_yp*sqrt(1/tauy)*randnorm[1];
}
if(tauy!=0.0) {
r6[3] -= 2*r6[3]/tauy;
}
if(espread!=0.0) {
r6[4] += 2*espread*sqrt(1/tauz)*randnorm[2];
}
if(tauz!=0.0) {
r6[4] -= 2*r6[4]/tauz;
}

if(EnergyLossFactor>0.0) {
r6[4] -= EnergyLossFactor;
}
}
}
}
Expand All @@ -68,7 +53,7 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
{
/* if (ElemData) {*/
if (!Elem) {
double emitx, emity, espread, taux, tauy, tauz, betax, betay, U0, EnergyLossFactor;
double emitx, emity, espread, taux, tauy, tauz, betax, betay;
emitx=atGetDouble(ElemData,"emitx"); check_error();
emity=atGetDouble(ElemData,"emity"); check_error();
espread=atGetDouble(ElemData,"espread"); check_error();
Expand All @@ -77,7 +62,6 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
tauz=atGetDouble(ElemData,"tauz"); check_error();
betax=atGetDouble(ElemData,"betax"); check_error();
betay=atGetDouble(ElemData,"betay"); check_error();
U0=atGetDouble(ElemData,"U0"); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->emitx=emitx;
Expand All @@ -90,19 +74,8 @@ ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
Elem->betay=betay;
Elem->sigma_xp=sqrt(emitx/betax);
Elem->sigma_yp=sqrt(emity/betay);
Elem->U0=U0;
Elem->EnergyLossFactor=U0/Param->energy;
}
SimpleQuantDiffPass(r_in, Elem->sigma_xp, Elem->sigma_yp, Elem->espread, Elem->taux, Elem->tauy, Elem->tauz, Elem->EnergyLossFactor, Param->thread_rng, num_particles);
/* }
else {
atFree(Elem->T1);
atFree(Elem->T2);
atFree(Elem->R1);
atFree(Elem->R2);
atFree(Elem->EApertures);
atFree(Elem->RApertures);
}*/
SimpleQuantDiffPass(r_in, Elem->sigma_xp, Elem->sigma_yp, Elem->espread, Elem->taux, Elem->tauy, Elem->tauz, Param->thread_rng, num_particles);
return Elem;
}

Expand All @@ -117,7 +90,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
double *r_in;
const mxArray *ElemData = prhs[0];
int num_particles = mxGetN(prhs[1]);
double emitx, emity, espread, taux, tauy, tauz, betax, betay, sigma_xp, sigma_yp, U0, EnergyLossFactor;
double emitx, emity, espread, taux, tauy, tauz, betax, betay, sigma_xp, sigma_yp;

emitx=atGetDouble(ElemData,"emitx"); check_error();
emity=atGetDouble(ElemData,"emity"); check_error();
Expand All @@ -127,19 +100,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
tauz=atGetDouble(ElemData,"tauz"); check_error();
betax=atGetDouble(ElemData,"betax"); check_error();
betay=atGetDouble(ElemData,"betay"); check_error();
U0=atGetDouble(ElemData,"U0"); check_error();
sigma_xp=sqrt(emitx/betax);
sigma_yp=sqrt(emity/betay);
EnergyLossFactor=U0/6e9;
if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix");
/* ALLOCATE memory for the output array of the same size as the input */
plhs[0] = mxDuplicateArray(prhs[1]);
r_in = mxGetDoubles(plhs[0]);
SimpleQuantDiffPass(r_in, sigma_xp, sigma_yp, espread, taux, tauy, tauz, EnergyLossFactor, &pcg32_global, num_particles);
SimpleQuantDiffPass(r_in, sigma_xp, sigma_yp, espread, taux, tauy, tauz, &pcg32_global, num_particles);
}
else if (nrhs == 0) {
/* list of required fields */
plhs[0] = mxCreateCellMatrix(9,1);
plhs[0] = mxCreateCellMatrix(8,1);
mxSetCell(plhs[0],0,mxCreateString("emitx"));
mxSetCell(plhs[0],1,mxCreateString("emity"));
mxSetCell(plhs[0],2,mxCreateString("espread"));
Expand All @@ -148,19 +119,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
mxSetCell(plhs[0],5,mxCreateString("tauz"));
mxSetCell(plhs[0],6,mxCreateString("betax"));
mxSetCell(plhs[0],7,mxCreateString("betay"));
mxSetCell(plhs[0],8,mxCreateString("U0"));


/*
if (nlhs>1) {
plhs[1] = mxCreateCellMatrix(6,1);
mxSetCell(plhs[1],0,mxCreateString("T1"));
mxSetCell(plhs[1],1,mxCreateString("T2"));
mxSetCell(plhs[1],2,mxCreateString("R1"));
mxSetCell(plhs[1],3,mxCreateString("R2"));
mxSetCell(plhs[1],4,mxCreateString("RApertures"));
mxSetCell(plhs[1],5,mxCreateString("EApertures"));
}*/
}
else {
mexErrMsgIdAndTxt("AT:WrongArg","Needs 0 or 2 arguments");
Expand Down
153 changes: 153 additions & 0 deletions atintegrators/SimpleRadiationPass.c
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@

#include "atelem.c"
#include "atrandom.c"

struct elem
{
double *damp_mat_diag;
double dispx;
double dispxp;
double dispy;
double dispyp;
double U0;
double EnergyLossFactor;
};

void SimpleRadiationPass(double *r_in,
double *damp_mat_diag, double dispx, double dispxp,
double dispy, double dispyp, double EnergyLossFactor, int num_particles)

{
double *r6;
int c, i;
double x,xp,y,yp,z,dpp;

#pragma omp parallel for if (num_particles > OMP_PARTICLE_THRESHOLD*10) default(shared) shared(r_in,num_particles) private(c,r6)
for (c = 0; c<num_particles; c++) { /*Loop over particles */
r6 = r_in+c*6;

if(!atIsNaN(r6[0])) {

dpp = r6[4];
x = r6[0] - dispx*dpp;
xp = r6[1] - dispxp*dpp;

y = r6[2] - dispy*dpp;
yp = r6[3] - dispyp*dpp;

z = r6[5];

if(damp_mat_diag[0]!=1.0) {
x *= damp_mat_diag[0];
}

if(damp_mat_diag[1]!=1.0) {
xp *= damp_mat_diag[1];
}

if(damp_mat_diag[2]!=1.0) {
y *= damp_mat_diag[2];
}

if(damp_mat_diag[3]!=1.0) {
yp *= damp_mat_diag[3];
}

if(damp_mat_diag[4]!=1.0) {
dpp *= damp_mat_diag[4];
}

if(damp_mat_diag[5]!=1.0) {
z *= damp_mat_diag[5];
}

r6[4] = dpp;
r6[5] = z;
r6[0] = x + dispx*dpp;
r6[1] = xp + dispxp*dpp;
r6[2] = y + dispy*dpp;
r6[3] = yp + dispyp*dpp;

r6[4] -= EnergyLossFactor;
}
}
}

#if defined(MATLAB_MEX_FILE) || defined(PYAT)
ExportMode struct elem *trackFunction(const atElem *ElemData,struct elem *Elem,
double *r_in, int num_particles, struct parameters *Param)
{
/* if (ElemData) {*/
if (!Elem) {
double U0, EnergyLossFactor;
double dispx, dispxp, dispy, dispyp;
double *damp_mat_diag;

damp_mat_diag=atGetDoubleArray(ElemData,"damp_mat_diag"); check_error();
U0=atGetDouble(ElemData,"U0"); check_error();
dispx=atGetOptionalDouble(ElemData,"dispx",0.0); check_error();
dispxp=atGetOptionalDouble(ElemData,"dispxp",0.0); check_error();
dispy=atGetOptionalDouble(ElemData,"dispy",0.0); check_error();
dispyp=atGetOptionalDouble(ElemData,"dispyp",0.0); check_error();

Elem = (struct elem*)atMalloc(sizeof(struct elem));
Elem->U0=U0;
Elem->dispx=dispx;
Elem->dispxp=dispxp;
Elem->dispy=dispy;
Elem->dispyp=dispyp;
Elem->EnergyLossFactor=U0/Param->energy;
Elem->damp_mat_diag=damp_mat_diag;
}
SimpleRadiationPass(r_in, Elem->damp_mat_diag, Elem->dispx, Elem->dispxp, Elem->dispy, Elem->dispyp, Elem->EnergyLossFactor, num_particles);
return Elem;
}

MODULE_DEF(SimpleRadiationPass) /* Dummy module initialisation */

#endif /*defined(MATLAB_MEX_FILE) || defined(PYAT)*/

#if defined(MATLAB_MEX_FILE)
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
if (nrhs == 2) {
double *r_in;
const mxArray *ElemData = prhs[0];
int num_particles = mxGetN(prhs[1]);
double U0, EnergyLossFactor;
double dispx, dispxp, dispy, dispyp;
double *damp_mat_diag;

damp_mat_diag=atGetDoubleArray(ElemData,"damp_mat_diag"); check_error();
dispx=atGetOptionalDouble(ElemData,"dispx",0.0); check_error();
dispy=atGetOptionalDouble(ElemData,"dispy",0.0); check_error();
dispxp=atGetOptionalDouble(ElemData,"dispxp",0.0); check_error();
dispyp=atGetOptionalDouble(ElemData,"dispyp",0.0); check_error();
U0=atGetDouble(ElemData,"U0"); check_error();
EnergyLossFactor=U0/6e9;
if (mxGetM(prhs[1]) != 6) mexErrMsgIdAndTxt("AT:WrongArg","Second argument must be a 6 x N matrix");
/* ALLOCATE memory for the output array of the same size as the input */
plhs[0] = mxDuplicateArray(prhs[1]);
r_in = mxGetDoubles(plhs[0]);
SimpleRadiationPass(r_in, damp_mat_diag, dispx, dispxp, dispy, dispyp, EnergyLossFactor, num_particles);
}
else if (nrhs == 0) {
/* list of required fields */
plhs[0] = mxCreateCellMatrix(2,1);
mxSetCell(plhs[0],0,mxCreateString("damp_mat_diag"));
mxSetCell(plhs[0],1,mxCreateString("U0"));

if (nlhs>1) {
/* list of optional fields */
plhs[1] = mxCreateCellMatrix(4,1); /* No optional fields */
mxSetCell(plhs[1],0,mxCreateString("dispx"));
mxSetCell(plhs[1],1,mxCreateString("dispxp"));
mxSetCell(plhs[1],2,mxCreateString("dispy"));
mxSetCell(plhs[1],3,mxCreateString("dispyp"));
}
}
else {
mexErrMsgIdAndTxt("AT:WrongArg","Needs 0 or 2 arguments");
}
}
#endif /*defined(MATLAB_MEX_FILE)*/
59 changes: 55 additions & 4 deletions pyat/at/lattice/elements.py
Original file line number Diff line number Diff line change
Expand Up @@ -1102,7 +1102,10 @@ def __init__(self, family_name: str, m66=None, **kwargs):
class SimpleQuantDiff(_DictLongtMotion, Element):
"""
Linear tracking element for a simplified quantum diffusion,
radiation damping and energy loss
radiation damping and energy loss.
Note: The damping times are needed to compute the correct
kick for the emittance. Radiation damping is NOT applied.
"""
_BUILD_ATTRIBUTES = Element._BUILD_ATTRIBUTES
default_pass = {False: 'IdentityPass', True: 'SimpleQuantDiffPass'}
Expand All @@ -1111,7 +1114,7 @@ def __init__(self, family_name: str, betax: float = 1.0,
betay: float = 1.0, emitx: float = 0.0,
emity: float = 0.0, espread: float = 0.0,
taux: float = 0.0, tauy: float = 0.0,
tauz: float = 0.0, U0: float = 0.0,
tauz: float = 0.0,
**kwargs):
"""
Args:
Expand All @@ -1126,7 +1129,6 @@ def __init__(self, family_name: str, betax: float = 1.0,
taux: Horizontal damping time [turns]
tauy: Vertical damping time [turns]
tauz: Longitudinal damping time [turns]
U0: Energy Loss [eV]
Default PassMethod: ``SimpleQuantDiffPass``
"""
Expand Down Expand Up @@ -1156,12 +1158,61 @@ def __init__(self, family_name: str, betax: float = 1.0,
if espread > 0.0:
assert tauz > 0.0, 'if espread is given, tauz must be non zero'

self.U0 = U0
self.betax = betax
self.betay = betay
super(SimpleQuantDiff, self).__init__(family_name, **kwargs)


class SimpleRadiation(_DictLongtMotion, Radiative, Element):
"""Simple radiation damping and energy loss"""
_BUILD_ATTRIBUTES = Element._BUILD_ATTRIBUTES + ['taux', 'tauy', 'tauz']
default_pass = {False: 'IdentityPass', True: 'SimpleRadiationPass'}

def __init__(self, family_name: str,
taux: float = 0.0, tauy: float = 0.0,
tauz: float = 0.0, U0: float = 0.0,
**kwargs):
"""
Args:
family_name: Name of the element
Optional Args:
taux: Horizontal damping time [turns]
tauy: Vertical damping time [turns]
tauz: Longitudinal damping time [turns]
U0: Energy loss per turn [eV]
Default PassMethod: ``SimpleRadiationPass``
"""
kwargs.setdefault('PassMethod', self.default_pass[True])

assert taux >= 0.0, 'taux must be greater than or equal to 0'
if taux == 0.0:
dampx = 1
else:
dampx = numpy.exp(-1/taux)

assert tauy >= 0.0, 'tauy must be greater than or equal to 0'
if tauy == 0.0:
dampy = 1
else:
dampy = numpy.exp(-1/tauy)

assert tauz >= 0.0, 'tauz must be greater than or equal to 0'
if tauz == 0.0:
dampz = 1
else:
dampz = numpy.exp(-1/tauz)

self.U0 = U0

self.damp_mat_diag = numpy.array([dampx, dampx,
dampy, dampy,
dampz, dampz])

super(SimpleRadiation, self).__init__(family_name, **kwargs)


class Corrector(LongElement):
"""Corrector element"""
_BUILD_ATTRIBUTES = LongElement._BUILD_ATTRIBUTES + ['KickAngle']
Expand Down
6 changes: 4 additions & 2 deletions pyat/at/lattice/lattice_object.py
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,8 @@
('collective_pass', elt.Collective, 'auto'),
('diffusion_pass', elt.QuantumDiffusion, 'auto'),
('energyloss_pass', elt.EnergyLoss, 'auto'),
('simplequantdiff_pass', elt.SimpleQuantDiff, 'auto')
('simplequantdiff_pass', elt.SimpleQuantDiff, 'auto'),
('simpleradiation_pass', elt.SimpleRadiation, 'auto')
),
True: (
('cavity_pass', elt.RFCavity, 'auto'),
Expand All @@ -63,7 +64,8 @@
('collective_pass', elt.Collective, 'auto'),
('diffusion_pass', elt.QuantumDiffusion, 'auto'),
('energyloss_pass', elt.EnergyLoss, 'auto'),
('simplequantdiff_pass', elt.SimpleQuantDiff, 'auto')
('simplequantdiff_pass', elt.SimpleQuantDiff, 'auto'),
('simpleradiation_pass', elt.SimpleRadiation, 'auto')
)
}

Expand Down
Loading

0 comments on commit d95aa5c

Please sign in to comment.