Improved performance of PropagationBP

[reichherzerp]

Description This pull request improves the behavior of the PropagationBP module in CRPropa for fixed step sizes, resulting in a speed improvement of up to a factor of 3. In addition, I have implemented a test to verify the correct behavior of my changes in the module. I also manually checked the output files and compared them to the current master.

Performance comparison: simulation time with grid-less turbulence: master: 75.05±0.11 s this PR: 25.43±0.02 s

simulation time with grid turbulence: master: 82.04±0.23 s this PR: 57.16±0.22 s

Code to reproduce:

import crpropa as crp

sim = crp.ModuleList()

# magnetic turbulence
turbulence_method = 'PW' # or 'grid'
b_field = crp.MagneticFieldList()
turbulence_spectrum = crp.SimpleTurbulenceSpectrum(1e-6*crp.gauss, crp.pc, 100*crp.pc)
if turbulence_method == 'grid':
    max_trajectory = 2e8*crp.pc
    grid_properties = crp.GridProperties(crp.Vector3d(0.,0.,0.), 256, crp.pc/2)
    turbulence = crp.SimpleGridTurbulence(turbulence_spectrum, grid_properties, 1)
else:  
    max_trajectory = 2e6*crp.pc
    turbulence = crp.PlaneWaveTurbulence(turbulence_spectrum, 10**3, 1)
b_field.addField(turbulence)

# propagation
prop_bp = crp.PropagationBP(b_field, crp.pc)
sim.add(prop_bp)
sim.add(crp.MaximumTrajectoryLength(max_trajectory))

# output
output = crp.TextOutput('output_'+turbulence_method+'.txt', crp.Output.Trajectory3D)

# observer
obs = crp.Observer()
n_obs = 3
dist = max_trajectory/n_obs
time_observer = crp.ObserverTimeEvolution(0, dist, n_obs)
obs.add(time_observer)
obs.setDeactivateOnDetection(False)
obs.onDetection(output)
sim.add(obs)

# source
source = crp.Source()
source.add(crp.SourcePosition(crp.Vector3d(0)))
source.add(crp.SourceParticleType(crp.nucleusId(1, 1)))
source.add(crp.SourceEnergy(1e16*crp.eV))
source.add(crp.SourceDirection(crp.Vector3d(1,0,0)))

# run simulation
sim.setShowProgress(True)
sim.run(source, 1, True)

[lukasmerten]

So basically with this PR simply the whole error calculation part is skipped as it would not be used anyway, right?

[reichherzerp]

Yes exactly, currently, we compute the propagation step (h) and the two small helper propagation steps (h/2) to get the error estimate that we do not need when using fixed steps. By removing the two helper steps, we save 2/3 of these time-consuming computations.

[lukasmerten]

I can roughly confirm the speed improvement. Using the above script with the SIMD optimization for the plane wave approach ("Fast waves"), decreases the performance a bit, but it still is a significant improvement.

So if you fix the typo and comment on my other question, I can approve this PR.

Nice performance improvement for "no cost" by the way!

[lukasmerten]

Looks good to me now.