The OPAL Discussion Forum

["Finn O'Shea" <finn.oshea AT nusano.com>]

Hi Jochem,

Thanks for the workarounds!  I am afraid I have found something else that seems like a bug.  

After a bend, the reference particle does not appear to encounter subsequent bends.  Here is what happens when my beamline has two bends (MDL_S.01 and TLA_B.01) with HAPERT=0.2:

The beam doesn't see the "middle" (body?) of the second bend magnet.  When I increase the HAPERT=2.2 the beam sees more of the downstream edge field:

But it still doesn't seem to see the middle of the bend and so it isn't bent though nearly a large enough angle.  The 2.2 meter half aperture is more than large enough to accommodate the bent trajectory of the beam in the second magnet, but the reference particle never sees it.

I'm defining the TLA_B.01 magnet using ELEMEDGE = and the _3D.opal file returns the following for that magnet (trimmed for clarity):

TLA_B.01 : SBEND, ANGLE = 25.0 / 180 * PI,              FMAPFN = "1DPROFILE1-DEFAULT",          X = -0.0909291338, Z = 3.1350046446, THETA = 345.0 * PI / 180,         L = 2 * bend_radius * sin(bend_angle / 2.0),  HAPERT = 2.2;

When I remove the second magnet and perform the simulation again (this time plotting ref_px), I find that the reference particle is not bent through the correct angle:

ref_px is almost exactly 1/2 of the value that it should be (the bend is a little too strong).  When I increase the charge from 2 to 3, the value is almost exactly 1/3 what it should be (again, a little too strong).  When I set the charge to 1, it is exactly what it should be.  It is pretty clear from this that the EM force on the particles due to the magnetic field in the SBEND is not being multiplied by the charge of the particle.  If I set the DESIGNENERGY for the SBEND to be the "correct" value and not the correction that you have given me, I get the correct value for ref_px when Q = 2*e.

If I define the magnets without your adjusted formulas for the design energy, I get what I expect for ref_px:

I guess what is happening is that the particle pusher was designed for singly charged particles and the work around to make it work for different charge values was to scale all the fields with the charge.  That would explain the magnetic field formula, the bending magnet strength issue and (I think) the bunch charge growth.    The last one is tricky for me to understand, because it isn't clear to me if the bunch charge should grow with Q or Q^2 for the space charge forces because in this case the charge is both the source of the fields and then there is the charge of the pushed particle to account for.  The manual says that the space charge pusher is a different algorithm from the external field pusher, so maybe one of the Qs is taken up in there somehow?

Finn

On Fri, May 22, 2020 at 6:18 AM Snuverink Jochem (PSI) <jochem.snuverink AT psi.ch> wrote:

Hi Finn,

Many thanks for the detailed email.

Unless I misunderstood something, I think you discovered two bugs in Opal-T.

1) qi_m, the charge per macroparticle is defined in TrackRun.cpp (last line):

    // Return charge per macroparticle.

    return beam->getCharge() * beam->getCurrent() / (beam->getFrequency()*1.0e6) / numberOfParticles;

The factor beam->getCharge() shouldn't be there.

2) As you discovered, the charge (only its sign) is not taken into account when calculating the bending field (in the 1DPROFILE1-DEFAULT case). This is the case both for RBend and SBend.

I'll make two issues in the issue tracker.

As a workaround for each case:

1) You can halve the current, since the current only used to calculate the bunch charge.

2) You can try to reduce the design energy of the bend such to halve beta*gamma:

new_design_energy = m0 / 2 * ( sqrt ( ( design_energy / m0 + 1 ) ^ 2 + 3 ) - 1 )

Best regards,

Jochem

On 22 May 2020, at 02:20, finn.oshea AT nusano.com wrote:

I've been trying to get a simulation going with alpha particles and the
various printouts suggest that charge isn't being handled consistently in
Opal-t.

I had a long version of this message with more output from Opal-t, but sympa
would not let me send it.

Starting with a FODO example that Andreas sent me, I can change the charge of
the bunch by changing the charge of the particles, in this case simply
changing from protons to carbon atoms increases the bunch charge by a factor
of 12, even though the current has not been changed.

Since the total bunch charge is defined through the relationship BCURRENT = Q
* BFREQ and neither BCURRENT nor BFREQ has changed, Q shouldn't either, but it
does.

This problem also seems to have crept in to the computation of the field
strength of an SBEND in a way that seems doubly incorrect to me.

Here is the input to Opal-t for a particle of charge +2 (I've deleted some
lines that do not matter, this is not the same beam as shown above):
OPAL>
OPAL> value: {EDES,GAMMA,BETA,P0,BRHO}
={0.05064,1.01359,0.163181,0.616503,1.02822}
OPAL>

  54 REAL switch_angle = 15.0 * Pi / 180.0; // 15-degree switch angle
  55 REAL switch_radius = 1; // 1 meter bending radius
  56
  57 MDL_S.01 : SBEND, ANGLE = switch_angle,
  58           FMAPFN = "1DPROFILE1-DEFAULT",
  59           ELEMEDGE = 2.664,
  60           DESIGNENERGY = design_energy,
  61           L = 2 * switch_radius * sin(switch_angle / 2.0),
  62           E1 = 0, E2 = 0, // make a "true" sector magnet

P0 and BRHO give the particle momentum and rigidity.

However, when I look at the output of the SBEND magnet (info level 2) I see
that the strength of the magnet is too large by a factor of 2:
SBend [2]> Reference Trajectory Properties
SBend [2]>
======================================================================
SBend [2]>
SBend [2]> Bend angle magnitude:    0.261799 rad (15 degrees)
SBend [2]> Entrance edge angle:     0 rad (0 degrees)
SBend [2]> Exit edge angle:         0 rad (0 degrees)
SBend [2]> Bend design radius:      1 m
SBend [2]> Bend design energy:      5.0655e+07 eV
SBend [2]>
SBend [2]> Bend Field and Rotation Properties
SBend [2]>
======================================================================
SBend [2]>
SBend [2]> Field amplitude:         2.05674 T
SBend [2]> Field index:  0
SBend [2]> Rotation about z axis:   0 rad (0 degrees)
SBend [2]>
SBend [2]> Reference Trajectory Properties Through Bend Magnet with Fringe
Fields
SBend [2]>
======================================================================
SBend [2]>
SBend [2]> Reference particle is bent: 0.2618 rad (15 degrees) in x plane
SBend [2]> Reference particle is bent: 0 rad (0 degrees) in y plane

Where I've deleted some lines that I don't think matter.  Changing the charge
of the particle does not change the output.  The field strength is still
2.05674 T.  When I look in to the code that produces this text (Bend2D.cpp) it
appears to get that value from the following function (line 1070 of the file):

void Bend2D::setBendStrength() {

   // Estimate bend field magnitude.
   double mass = RefPartBunch_m->getM();
   double gamma = designEnergy_m / mass + 1.0;
   double betaGamma = std::sqrt(std::pow(gamma, 2.0) - 1.0);
   double charge = RefPartBunch_m->getQ();

   fieldAmplitude_m = ((charge / std::abs(charge)) * betaGamma * mass /
                       (Physics::c * designRadius_m));


   // Search for angle if initial guess is not good enough.
   findBendStrength(mass, gamma, betaGamma, charge);
}

This function is interesting for two reasons:
(1) The formula being used is q*P0 / (rho*c) where it should be P0/(q*rho*c)
(2) For some reason the charge is divided by its absolute value, perhaps the
author just wanted the sign of the charge
This formula will always return the wrong value for particles of charge
greater than 1 (in an absolute sense).

Nevertheless, I tracked down the getQ() function to the file PartBunchBase.hpp
(which also happens to be the source of the BUNCH output producing the
erroneous bunch charge value), so perhaps that is the root of the problem.
But it does appear that the incorrect charge output is not just a superficial
print bug, but rather the value for the bunch charge really is incorrect in
the variable qi_m.  I can't figure out where qi_m comes from, I've never used
C++.

Thanks for any assistance,
Finn

--

Senior Scientist

Nusano, Inc

28575 Livingston Ave

Valencia, CA 91355

Office: 1 (424) 293-3174