PartBunch::get_bounds can produce NaNs

Can you please check if OPAL-1.6 shows the same behaviour?

And yes I agree, OPAL should catch that nicer if it is an input file (user) error

The problem here isn't get_bounds but that one or several particles have coordinates that are NaN

The simulation was working on the 1.6 branch before I upgraded to the 2.0 syntax, but I'll check on Monday on 1.6 to be sure.

So, yes this is working on 1.6 both single- and multi-particle.

For 2.0 the only relevant change (for single particle) I can see is this line from commit 7f73ff9

- DistSP: DISTRIBUTION, DISTRIBUTION = fromfile, FNAME = "_spdistribution_",
+ DistSP: DISTRIBUTION, TYPE = fromfile, FNAME = "_spdistribution_",

with spdist.opal the same being:

1
0.0         0.0           0.0           0.0          0.0           0.00

The std output is the same until the bunch information, for 1.6:

OPAL> * Phase space dump -1 (global frame) at integration step 0
OPAL> * T = 0 ns, Live Particles: 1
OPAL> * E = 71.6 MeV, beta * gamma = 0.39805
OPAL> * Bunch position: R =  2037 mm, Theta = 110 Deg, Z = 0 mm
OPAL> * Local Azimuth = 112.36 Deg, Local Elevation = 0 Deg
OPAL> 
OPAL> * *********************** Bunch information in global frame: ***********************
OPAL> * ************** B U N C H ********************************************************* 
OPAL> * NP              =   1
OPAL> * Qtot            =    1.60218e-10 [nC]       Qi    =  1.60218e-10 [nC]
OPAL> * Ekin            =    7.16000e+01 [MeV]      dEkin =  0.00000e+00 [MeV]
OPAL> * rmax            = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [m]
OPAL> * rmin            = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [m]
OPAL> * rms beam size   = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [m]
OPAL> * rms momenta     = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [beta gamma]
OPAL> * mean position   = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [m]
OPAL> * mean momenta    = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [beta gamma]
OPAL> * rms emittance   = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) (not normalized)
OPAL> * rms correlation = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 )
OPAL> * hr              = ( -1.00000e+00 , -1.00000e+00 , -1.00000e+00 ) [m]
OPAL> * dh              =    1.00000e-10 [%]
OPAL> * t               =    0.00000e+00 [s]        dT    =  5.48426e-11 [s]
OPAL> * spos            =    0.00000e+00 [m]
OPAL> * **********************************************************************************

For master (Qtot, rmin, rmax different):

OPAL> * Integration step 0 (no phase space dump for <= 2 particles)
OPAL> * T = 0 ns, Live Particles: 1
OPAL> * E = 71.6 MeV, beta * gamma = 0.39805
OPAL> * Bunch position: R =  2037 mm, Theta = 110 Deg, Z = 0 mm
OPAL> * Local Azimuth = 112.36 Deg, Local Elevation = 0 Deg
OPAL> 
OPAL> * *********************** Bunch information in global frame: ***********************
OPAL> * ************** B U N C H ********************************************************* 
OPAL> * NP              = 1
OPAL> * Qtot            =        0.000 [fC]         Qi    =        0.000 [fC]
OPAL> * Ekin            =      71.600 [MeV]         dEkin =        0.000 [eV]
OPAL> * rmax            = (     -0.69670 ,      1.91415 ,      0.00000 ) [m]
OPAL> * rmin            = (     -0.69670 ,      1.91415 ,      0.00000 ) [m]
OPAL> * rms beam size   = (      0.00000 ,      0.00000 ,      0.00000 ) [um]
OPAL> * rms momenta     = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [beta gamma]
OPAL> * mean position   = (      0.00000 ,      0.00000 ,      0.00000 ) [um]
OPAL> * mean momenta    = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) [beta gamma]
OPAL> * rms emittance   = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 ) (not normalized)
OPAL> * rms correlation = (  0.00000e+00 ,  0.00000e+00 ,  0.00000e+00 )
OPAL> * hr              = (      0.00001 ,      0.00001 ,      0.00001 ) [um]
OPAL> * dh              =   1.00000e-10 [%]
OPAL> * t               =        0.000 [fs]         dT    =       54.843 [ps]
OPAL> * spos            =        0.000 [um]
OPAL> * **********************************************************************************

The tracking then is for 1.6:

OPAL> * --------------------------------- Start tracking ------------------------------------ *
OPAL> * Cavity MAIN2 Phase= -5.2776 [deg] transit time factor=  0.96939 dE= 0.45362 [MeV] E_kin= 72.054 [MeV] Time dep freq = 1
OPAL> * Cavity MAIN3 Phase= -5.2096 [deg] transit time factor=  0.96956 dE= 0.46962 [MeV] E_kin= 72.523 [MeV] Time dep freq = 1
OPAL>

and for master:

OPAL> * --------------------------------- Start tracking ------------------------------------ *
OPAL > RING: particle 0 out of the field map boundary!
OPAL > RING: Coords: ( -nan , -nan , -nan )
OPAL> * SPT: The particle was lost at step 0!
Error> 
Error> *** User error detected by function "ParallelCyclotronTracker"
Error> *** in line 311 of file "Ring.in" at end of statement:
Error>      RUN,METHOD="CYCLOTRON-T",BEAM=BEAM1,FIELDSOLVER=FS1,DISTRIBUTION=DISTSP;
Error>     The particle is out of the region of interest.
Error>

Let me know if there is anything else I can check.

added ~48 label

The difference in Qtot, rmin, rmax different are most likely only in the notation (scientific vs. fixed). Do you have the input file on merlin somewhere?

added Bug label

okay. by the way I noticed that hr is also different by a factor 10^12 (if the units are correct).

I'll check if this still exists and put the input file on Merlin.

Okay seems still to be there. As mentioned in the original report, this can be reproduced with the master of https://gitlab.psi.ch/AMAS-BDModels/PSI-Ring (careful 5GB(!)), load setupEnv.sh and runOpal.py --nobatch.

I have also put a version on Merlin at /gpfs/home/snuverink_j/PSI-Ring-master/PSI-Ring

@adelmann assumes a problem in the input file, so I attach it here (for single particle): Ring.in

Compared to 1.6 the main change that I see is the distribution command which changed its syntax:

DistSP: DISTRIBUTION, TYPE = fromfile, FNAME = "spdist.opal",
	INPUTMOUNITS = EV;

I have reduced the input file to the essentials: Ring.in

I have also checked the regression test RingCyclotron (which is by itself working for me) and changed the distribution from gauss to binomial (like in the example) and I get the same error.

The diff wrt to the test:

< Dist1:DISTRIBUTION, TYPE=binomial,
---
> Dist1:DISTRIBUTION, TYPE=gauss,

Also I modified it to run single particle and it gives the same error.

The single particle distribution: spdist.opal

The diff wrt to the test:

< DistSP: DISTRIBUTION, TYPE = fromfile, FNAME = "spdist.opal",
< 	INPUTMOUNITS = EV;
< 
< Fs1:FIELDSOLVER, FSTYPE=NONE, MX=32, MY=32, MT=16,
---
> Fs1:FIELDSOLVER, FSTYPE=FFT, MX=32, MY=32, MT=16,
83c80
< beam1: BEAM, PARTICLE=PROTON, pc=P0, SPACECHARGE=true, NPART=1, BCURRENT=1.0E-3, CHARGE=1.0, BFREQ= frequency;
---
> beam1: BEAM, PARTICLE=PROTON, pc=P0, SPACECHARGE=true, NPART=1e5, BCURRENT=1.0E-3, CHARGE=1.0, BFREQ= frequency;
88,89c85
< 
< RUN, method = "CYCLOTRON-T", beam=beam1, fieldsolver=Fs1, distribution=DistSP;
---
>  run, method = "CYCLOTRON-T", beam=beam1, fieldsolver=Fs1, distribution=Dist1;

Therefore, I think we can conclude it is not an input file problem.

Even more can we conclude it is the Binomial distribution?

It also happens for the single particle case.

I think we should wait fixing this, until the integrators are fixed (@frey_m)

I finished for the moment #124 (closed). You can proceed ...

Just to say that this is still there with the current master (db76ed47). The way I would probably debug this is to find the commit that caused this behaviour, which I might do at some point. Open to suggestions.

The problem comes most probable from the Trimcoils when in single particle tracking mode.

@snuverink_j can you confirm?

The regression test RingCyclotron has no trim coils and shows the same issue in single particle mode: RingCyclotron.in. Would that contradict it?

Ok but the RingCyclotron regressiontest works with the master, however with works at the moment means the test is running

Yes, but it doesn't when one changes the distribution from gauss to either single particle or binomial (#110 (comment 1916)).

Ok so we have a problem in single particle mode and when using the binomial distribution.

@snuverink_j I ran the RingCyclotron test with binomial distribution and in single particle mode (only serial supported) and both cases worked fine. Could you check?

I checked with current master version and for me it works fine now as well. However, the original still gives problems. I'll try to break it down to see where the difference is. As far as I can see now it is in the difference in the cyclotron definition.

When I just add trim coils to the RingCyclotron test, I get it again, now also with gaussian beam.

The output changed a bit and it might happen in a different place now, but I still get the NaNs. I noticed that when I switch to the LF-2 integrator it seems all fine.

Okay after some searching I think I found the culprit in Cyclotron::applyTrimCoil():

Printout of some variables:

OPAL> r 2.13 z -8.4701e-22
OPAL> slope 1
OPAL> tcr1 2131 tcr2 2132.1
OPAL> magnitude 0
OPAL> part1 1.0325e-64 part2 inf
OPAL> part3 -2.1397e-65 part4 -nan
OPAL> dr -nan
OPAL> br -nan bz -14.716

It seems r is in meters, while tcr1 is in mm. This change should be made consistently for all trim coils.

In addition of the unit fix, it would perhaps be good to protect the method against infinites.

p.s. I think another issue was fixed in the mean time with the Distribution since we had some NaNs without trim coils too.

Yes I worked on the distribution particularly on the binomial distribution.

I've cleaned the code a little bit, added a few checks and convert the units from mm to m in the code. For the user doesn't change anything except that I removed the attributes TCR1V, TCR2V, MBTCV and SLPTCV and changed the type of TCR1, TCR2, MBTC and SLPTC to real arrays.

With these changes the adjusted version of RingCyclotron runs smoothly.

@adelmann in the code there was (and still is) the statement

if (r < tcr2) {
//do the computation
}

I've changed it to (it's not actually written like this any more but the meaning is still the same)

if (r < tcr2 && r > tcr1) {
//do the computation
}

Could you look at this line and decide whether this makes sense?

Perfect, makes very much sense to me.

Sorry everyone, the m/mm problem is likely a remnant from when I changed the CyclotronTracker to use m internally during the Erice meeting.

@kraus: Many thanks for the fixes! I confirm this solves the issue. I had to add a missing include in order to compile (e56cc0f7).

closed

Hello, I am going to ask a couple questions. First 3 trim cells are a bit strange in that coils 1, 2, and 3 are nested, in that 1 sits inside of 2 and 3, does this impact the new conditions for the trim coils. Second if you change the code to what was listed above, are the condition for the trim coils still in place for if the beam is outside the trim coil? As that has a critical impact on getting the orbits right as the magnetic field reverses.

Sorry I don't understand you question. Where do you see 3 trim coils in the code? The old condition said that we don't compute the field if the particle is outside the outer limit. The one I added means that we also don't compute the field if it is inside the inner limit. Since I'm not a cyclotron expert I don't know whether this second condition makes sense...?

@kraus by the way the slopes are in kG/mm(?) and should be converted too I think.

Also some of the constants in Cyclotron::applyTrimCoil might have units. For example if slope is in kG/mm then x01, x02 should be in kG*mm^2. I think the original author (@adelmann?) should have a look at this.

So I guess what @pogue_n suggests is that we shouldn't cut when the particle is outside of the trim coil as there is still a remaining field that is relevant.

This is also reflected in (comments added by me):

if (2 * r < (tcr2 + tcr1)) { // inside the trim coil
...
} else { // outside the trim coil
...
}

I think the 3 trim coils @pogue_n is referring to are the first 3 coils in the PSI-Ring, and are an example that trim coils can be "nested". Therefore, the trim coil end values tcr2 are not necessarily descending (this is assumed in the current implementation).

ah sorry I just saw that the trim coils are sorted by their end value tcr2, so that should be okay then.

@snuverink_j the comments you added aren't quite right: (tcr1 +tcr2) / 2 is the center, thus if 2*r < tcr1 + tcr2 then the particle is inside the center.

reopened

@kraus yes you're right, sorry don't know what I was thinking there.

I don't think that the units are so wrong. From the formulas for P1 and P2 one can deduce that either

The first choice seems more probable but from a variable called slope I would expect that its unit is per length.

If the variables Amin, Amax1 and Amax2, x1, x2, h1, h2, justAnotherFudgeFactor are unitless, slope has unit length and magnitude has unit magnetic flux density then the unit of btr is magnetic flux density and the unit of dr is magnetic flux density per length. This would result in the desired units for br and bz.

@adelmann could you please check the variables, their units and their meaning and rename them such that someone else (such as me) understand them as well?

assigned to @pogue_n

@pogue_n you worked last on that. could you please do an sanity check and follow up what @kraus is suggesting?

Hi, I took another quick look. It appears to me that the slope is in units of dx^2/dB or how fast the slope is changing. Xo1 and Xo2 Look like they might be in B/dx as they essentially set the peak height and transition speed of the magnetic field as you are travelling through the virtual conductor with respect to the stated b field which is not ever reached unless Xo1 and x02 are adjusted properly. h looks to be the shape of the coil geometry, thus it shapes the profile of the field as it transition through the virtual conductor region on each side, such as the field on a dipole is a different profile inside and outside the coil - this effectively limits how it can tail off from the conductor. H would have to be dx/dB then.

assigned to @snuverink_j

@snuverink_j can you check that and update the code?

So if I understand it correctly slope = (d^2 B/dx^2 )^-1 ?

And Amin, Amax1 and Amax2 are unitless? What is there meaning?

So if the slope is (d^2 B/dx^2 )^-1 , then its unit mm^2 / kG ? .. Note that the manual actually says (chapter 8, ~p.96) : "The slope of the rising edge (kG/mm)"

@adelmann, @pogue_n : is there any paper that describes these formulas (where does the radial profile plot from #110 (comment 2213) come from)?

Right now the code gives for tcr1=3.29329 tcr2=3.49129 mb=0.014 slope=6 (trim coil 9 from the PSI-ring, @pogue_n: what values did you take for your plot?) basically a zero profile.

However, if I put tcr1/2 in mm I get:

This looks similar but I don't have the dip in the middle, is that expected?

Now if I put tcr1, tcr2 in meter (3.3 and 3.5), and I multiply slope by 0.001 (=0.006) and leaving everything else fixed, I get the same profile as before:

(in the plot the radius is multiplied by 1000 to make mm again)

That would suggest the slope is in mm/kG? Consequently h1, h2 have unit 1/kG, and x01, x02 unit kG, and Amax1, Amax2, Amin are unitless.

PartBunch::get_bounds can produce NaNs

Designs

Child items ...

Activity

PartBunch::get_bounds can produce NaNs

Relates to

Activity