Skip to content

Forbid negative lattice or element energies#1066

Merged
swhite2401 merged 3 commits intomasterfrom
forbid_negative_energy
Mar 28, 2026
Merged

Forbid negative lattice or element energies#1066
swhite2401 merged 3 commits intomasterfrom
forbid_negative_energy

Conversation

@swhite2401
Copy link
Copy Markdown
Contributor

This PR solves issue #1039, atelem.c is reformatted to forbid negative energies, all passmethods requiring energy field are modified to correctly report errors.

Comment thread atintegrators/atelem.c Outdated
atError("Energy must be positive. Check lattice or passmethod parameters.");
return 0.0; /* Never reached but makes the compiler happy */
}
}
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In python, the energy is taken only from the lattice, never from the element (when reading .mat files, the energy field is removed from all elements). So you need different versions for Matlab and python. The python version should be:

double atEnergy(double ringenergy, double elemenergy)
{
    if (ringenergy > 0.0)
        return ringenergy;
    else {
        atError("Energy must be positive. Check lattice or passmethod parameters.");
        return 0.0;
    }
}

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This makes element_pass crash for passmethods requiring energy, any idea?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Starting from your branch, I set this in the python section of atelem.c:

double atEnergy(double ringenergy, double elemenergy)
{
    if (ringenergy > 0.0)
        return ringenergy;
    else {
        atError("Energy not defined.");
        return 0.0;
    }
}

double atGamma(double ringenergy, double elemenergy, double rest_energy)
{
    if (ringenergy > 0.0) {
        if (rest_energy == 0.0)
            return 1.0E-9 * ringenergy / __E0;
        else
            return ringenergy / rest_energy;
    }
    else {
        atError("Energy not defined.");
        return 0.0;
    }
}

And I get the expected result (no crash):

(pyat) laurent@garou at % python
Python 3.10.20 | packaged by conda-forge | (main, Mar  5 2026, 17:06:34) [Clang 19.1.7 ] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import at
>>> import numpy as np
>>> dip = at.Dipole('D', 2.0, 1.0, k=0.0, PassMethod='BndMPoleSymplectic4RadPass')
>>> pout = dip.track(np.zeros(6))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/laurent/dev/libraries/at/pyat/at/tracking/track.py", line 419, in element_track
    rout = _element_pass(element, r_in, **kwargs)
  File "/Users/laurent/dev/libraries/at/pyat/at/tracking/utils.py", line 109, in wrapper
    return func(lattice, r_in, *args, **kwargs)
  File "/Users/laurent/dev/libraries/at/pyat/at/tracking/track.py", line 93, in _element_pass
    return _elempass(element, r_in, **kwargs)
ValueError: Energy not defined.
>>> 

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Strange... with the RFCavity it did.... let me try again

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I confirm that for me it crashes for dipole as well..., ringenergy is not defined in an element so in fact I don't know how if can work for you....

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It does because you give the kw argument energy to the track() function. Shouldn't we allow the user to alternatively provide the Energy element attribute as in the example I provided above?

I don't think this costs anything in terms of efficiency since this second check is ran only when the ringenergy is not defined which basically never happens.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It' a question of logic:

Defining the energy in element attributes introduces the risk of having different energies in different elements, that's why when starting pyat we removed it systematically. Of course this does not prevent from explicitly adding such an attribute, like any other custom attribute.

But when tracking comes in, I still think that the energy is related with the tracked particle rather than with the element itself. I very much prefer giving it as a keyword in Lattice.track or Element.track rather than as an element attribute… Similarly I would not set the particle as a element attribute, while it is accepted as a keyword in Lattice.track and Element.track. Even if it is still not used…

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok understood as long as this is clearly defined this is fine for me, is there a way to set a meaningful error for element_pass, for instance if the ringenergy is really undefined (NaN)?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The error message could be something like:
"Energy needs to be defined and positive. Check lattice parameters and tracking options."

Strictly speaking, the energy should be greater than the particle rest energy, but we may forget that…

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I have changed it to something similar

Comment thread atintegrators/atelem.c
Comment on lines +354 to +362
double atGamma(double ringenergy, double elemenergy, double rest_energy)
{
double energy = atEnergy(ringenergy, elemenergy);
if (rest_energy == 0.0)
return 1.0E-9 * energy / __E0;
else
return energy / rest_energy;
}

Copy link
Copy Markdown
Contributor

@lfarv lfarv Mar 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For python, the macro as it is is enough (and more efficient).

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not so sure about this, how about a simulation of energy ramp where the energy is set to negative during the ramp? I think we still need the energy check in atGamma()

atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
check_error();
}
Energy = atEnergy(Param->energy, Energy); check_error();
Copy link
Copy Markdown
Contributor

@lfarv lfarv Mar 26, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The test as it is was implemented by @oscarxblanco in #1016 to solve the same issue as you do.
What you write here is exactly what I suggested in #1016 (comment) but was not implemented. So it's ok for me

@swhite2401
Copy link
Copy Markdown
Contributor Author

I have implemented the requested change

Copy link
Copy Markdown
Contributor

@lfarv lfarv left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That looks good !

@swhite2401 swhite2401 merged commit 82eb1a1 into master Mar 28, 2026
22 checks passed
@swhite2401 swhite2401 deleted the forbid_negative_energy branch March 30, 2026 08:05
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants