Murua Steppers (and Progress in General)

[jsuagee]

Small Report of Current Work:

I have written 3 working Nystrom steppers, and 2 promising Nystrom steppers that are still in development, and now I would like to be able to evaluate their performance in a full stack environment with geometry. Specifically, I would like to be able to measure the global error vs. the number of function calls to the Mag Field for each of these steppers.

Plan:

• Stepper1 will be the Nystrom stepper that I want to evaluate.

• Stepper2 will be ClassicalRK4 (or some other reliable stepper).

What I am going to do is generate an almost exact solution trajectory with stepper2 by using a really small step size.

Then I will generate a trajectory for stepper2 using a starting moderate step size (which is an independent variable).

Then I will compare the points on trajectory2 with the points on trajectory1 using an interpolant of either stepper1 or stepper2. What I have been focusing on now is using an interpolant of stepper1 (the Nystrom stepper I am testing). But I could do it the other way around and use an interpolant for stepper2 instead [Somnath: Do you have yet a reliable RK stepper with interpolation working?].

(I have to use interpolants because the arclength values of the trajectory points of the two steppers will most likely not coincide when using the full stack.)

There is one main challenge to this plan: Many steps that the stepper performs will be rejected higher up in the stack. (Much of these rejected steps come from calculating the intersection points of the trajectory of the stepper with volumes of the geometry.) In order to record only steps that were not rejected I have to modify several of the higher up classes in the stack to alert the stepper when the steps are rejected.

I created a separate class to handle the recording of steps, StepTracker. I have only so far been testing it with the program testPropagateMagField. I have had to modify that program slightly in order to create a StepTracker object and pass a pointer to it to the other classes in the stack (every class besides G4PropagatorInField, which is outside of the scope of the project). Every class object - G4ChordFinder, G4MagIntegratorDriver, and the stepper - have references to the StepTracker object. The StepTracker object has several member functions which the other stack objects can call to alert it when steps are accepted/rejected.

(It turns out that this whole scheme is possible without modifying G4PropagatorInField because steps are only accepted within a call to G4ChordFinder::AdvanceChordLimited(), which is called by G4PropagatorInField: ComputeStep() and no where else. So I can alert the StepTracker when we enter and when we leave the scope of AdvanceChordLimited(), and to only record steps in that scope. I just finally figured this out yesterday after looking at a lot of the code to handle calculating the intersection points of the trajectory with Volumes of the geometry.)

So I am still testing this idea to see if it will work. When I tried it on the testPropagateMagField program there were no calls to the stepper outside of G4ChordFinder::AdvanceChordLimited(), which means that the trajectory did not encounter any possible intersection points. (Is this possible John?) I plan to test the StepTracker on the NTST program because I am sure that there are possible intersection points that occur there.

Previous Work:

I have created a template wrapper classe for the Nystrom Steppers, MagIntegratorStepper_byTime, whose only real job is to convert back and forth between momentum and velocity coordinates and scale the step length back and forth betwen arclength and time units. It functions as an intermediary between the rest of the stack, which was written with momentum coordinates and arc length in mind and the Nystrom stepper, which can only work when using velocity coordinates and when integrating with respect to time. So currently there is some extra overhead (roughly 6-12 extra multiplications per step call). Perhaps this can be avoided by tweaking the stack. Or maybe by a careful inspection of the code in the rest of the stack, this extra overhead can be avoided with by just changing the input to the top level function G4PropagatorInField::ComputeStep().

I had to write a separate template wrapper class for the Magnetic field. This also acts to scale each Right hand side evaluation. This template class overrides G4MagneticField::EvaluateRhsGivenB(). So this also introduces extra overhead of 1 multiplication per ComputeRightHandSide() call from the stepper (we have to perform a division by the relativistic mass of the particle).

Also, I have as yet been assuming constant relativistic mass (equivalently zero Electric Field).

I have also written several python routines for the ipython notebook, which can be used to produce plots (both 2d and 3d).

The two Nystrom methods still under development are almost done except for their interpolation methods, and their DistChord() methods (which use the interpolants).

[jsuagee]

I believe I got the above plan to work. I pushed the code, but I still have to clean it up. I only pushed it now because I was afraid to loose it in case my computer was damaged.

I would have had this much faster I think, but I made a typo in a flag: "#define Tracking" instead of "#define TRACKING". This simple mistake cost me several hours of frustration.

I will clean it up tomorrow (and document), and then try to move towards getting it to work with the NTST program.

[jsuagee]

Estimating error was complicated. Here are 2 plots of some error estimates. Position error is self explanatory. The X-axes is how far the particle has traveled and the Y-axes is the magnitude of the error in arc length units (mm). I am comparing them all against a run of the Fine45 stepper with step size really small (0.01 mm) so that it approximates really well the true trajectory (I will call this run of Fine45 the "Base Trajectory"). I chose this because I could get the interpolant to work with the Base Trajectory. (If you don't use interpolation, and instead pick the closest point in Base Trajectory, you will get a much less informative plot since the position error is on the scale of 10^-12 mm but the step size in the Base Trajectory is of the order 0.01 mm.)

I generated the Base Trajectory by using the compareByTime program, which is a modified testH program (only the stepper, not the full stack). That way I could have active control of the step size and keep it down to 0.01 mm. Of course all the starting parameters were made the same.

Also, I turned the geometry off. The only volume in tPMF is the World Box. If you turn on the geometry you get spikes in the error plots when you cross into a new volume. The spikes are large, but this is only because the Base Trajectory was calculated without the full stack and it doesn't deviate from its course when it passes into a new volume.

The steppers are undifferentiated in the momentum error plot. I don't know exactly why this happens (maybe it is an error on my part, or maybe there is some other reason). The plot is not a straight line though, there is a little bit of wiggle. I also computed everything in terms of velocity and then scaled by the mass back to momentum coordinates.

I chose a quadropole mag field, but I will construct the same plots for the uniform mag field.

The last thing is the plot of the required function evaluations. This doesn't make much sense to include without the full geometry in place, but I also did the test with the full geometry in tPMF, and the plot looked pretty much similar (although there were only a few intersections because the geometry of tPMF is not that complicated). Also, the function evaluations recorded for the Murua5459 stepper assume that you are using interpolation inside of the DistChord() method. I didn't have time to program this yet, so I just kept track of how many function calls were coming from DistChord() and I disregarded them in the sum. This is just to see what's possible in the optimal sense.

In order to get Murua5459::DistChord() working with interpolation, I will have to solve a small linear system to find the right coefficients, which I can do but it will take a day probably (given my track record).

Overall the two new steppers, Murua and Fine45, don't do that bad in this simple test. I wish I would have been able to get this working 2 or 3 weeks ago but there were all these little hard to see bugs in the design which took that long to get out.

[jsuagee]

There is a new tag "alpha-03". That tag is a mess, but I wanted to push it just in case of WWIII and my computer exploded. I'm working on cleaning it up.

[jonapost]

Thanks for the progress report and the plots.
A few comments:

• It would be great to see how curved the path is. This can be seen by plotting the dot product of the initial momentum with the momentum along the curve. (You can subtract 1 from this in case it is a very small value difference.)
• I do not understand the second plot at all - the deviation of the momentum. Could you please recheck it ? It seems that either the baseline stepper is slowly going wrong (at a larger rate than the deviations of the rest) or something else is going wrong.\
• Could you please re-do plot #3 (number of calls vs length) ? Unless you give a large length for the integration, its result is just the factor of different number of stages! All steppers likely succeed for all steps!

[jsuagee]

[jonapost]

Also plot #3 is more interesting if you vary the accuracy required - and use that accuracy as the 'x' axis. ( Once the length of the arc is more interesting - say deviations of at least 1.5 * pi.

[jsuagee]

The first 3 plots are for a Uniform field (I also included a 3dplot of the trajectory). Compute Step Size is 2000 mm. This is with all the geometry in tPMF, so intersections are checked for etc... Error would be lower in general without corrections for intersections with volumes.

Whatever was causing the dubious momentum error plots is not occurring anymore:

Trajectory (from tPMF). As you can see it makes at least one full rotation. The Radius of curvature is 970.4907... mm:

These similar plots but for the Quadropole Mag field:

The number of function calls used by the Murua stepper shoots up after around 4000 mm, which is when the trajectory starts to double back on itself. I haven't had a chance to think about why that happens?

[jonapost]

Thanks for the new plots. Using a more challenging initial distance makes it more interesting to my eyes at least.

The plots for the uniform magnetic field look consistent - the errors are very good for the Murua stepper, and generally grow with some noise for other steppers.

A few questions:

• It is not immediately clear what the step sizes were for the different steppers. Would it be possible to make a note of them - or even try to display it in a new version of the diagrams?
• how are the errors estimated for the points between the endpoints of the integration steps?

[jonapost]

The errors look very large in the quadrupole field.

Even more important:

• What was the error tolerance used for these runs (and for the uniform field)?

[jsuagee]

Well the errors are estimated by Fine's 4th order interpolation scheme (you construct an interpolated point with the same time/arclength value as the point you are trying to estimate the error for). StepTracker will have stored the arclength and time values of each point that we want to estimate, so this is possible.

Yes it is possible to plot the different step size's chosen by the full stack, for each of the different steppers.

You know, I don't know what the error tolerance was for these runs. It was whatever was specified in the version of tPMF that I was using (I didn't change those values).

[jsuagee]

The step size chosen for the baseline solution is currently 0.1 mm (it was 0.01 mm before). The stepper chosen for the base line solution was Fine's RKN45 method. This was chosen only because I had a ready interpolation scheme working for it.