- Defined in File PlaneWaveTurbulence.h
PlaneWaveTurbulence: public TurbulentField¶
Interpolation-free turbulent magnetic field based on the GJ99 and TD13 papers.
This class provides a turbulent magnetic field that is generated as described by (Giacalone and Jokipii, 1999) and (Tautz and Dosch, 2013). Instead of using an inverse Fourier transform to generate the field on a grid which then needs to be interpolated to obtain in-between values this method only generates the wave modes making up the turbulent magnetic field ahead of time. At run time, when the field’s value at a particular position is required, these plane waves are then evaluated analytically at that position. This guarantees that the resulting field is completely free of divergence, reproduces the mean field strength accurately, and does not suffer from other interpolation-induced problems. The disadvantage is that the number of wave modes is drastically smaller when compared with initTurbulence, which might have physical ramifications on the particles propagating through the field. Furthermore, depending on the exact setting, the implementation is somewhat slower (see Schlegel et al. for details).
Using the SIMD optimization
In order to mitigate some of the performance impact that is inherent in this method of field generation, an optimized version is provided. According to our tests (see the paper above), this version runs 20-30x faster than the baseline implementation and matches the speed of trilinear interpolation on a grid at a bit less than 100 wavemodes. To do this, it uses special CPU instructions called AVX, which are unfortunately not supported by every CPU. On Linux, you can check whether your CPU supports AVX by using the
lscpucommand, and searching the flags section for the string “avx”. Alternatively, if you are building on the same CPU that CRPropa will run on, you can have the compiler detect this for you automatically (see below).
An additional speedup of about 1.33 can be achieved if the CPU supports the FMA extension in addition to AVX. Again, you can either check for this manually, or have the compiler figure it out for you.
Note** that the optimized and non-optimized implementations to not return the exact same results. In fact, since the effective wave numbers used by the optimized implementation are very slightly different from those used by the non-optimized version (a difference smaller than the precision of a double, but nevertheless relevant at some point), the wavemodes go out of phase for large distances from the origin, and the fields are no longer comparable at all.
If you are building on the same machine that the code will run on:
- In cmake: enable the FAST_WAVES flag.
- Also in cmake: set SIMD_EXTENSIONS to “native”. The compiler will automatically detect support for your CPU and run the build with the appropriate settings.
- Generate files and exit cmake, then build.
- If your CPU does not support the necessary extensions, the build will fail with an error telling you so. In this case, you won’t be able to use the optimization; go back into cmake, disable FAST_WAVES, and build again.
- If the build runs through without errors, the code is built with the optimization.
If you are building on a different machine from the one where the code will run:
- Figure out which extensions your target machine supports, using
lscpu | grep avxfor AVX, and
lscpu | grep fmafor FMA. Run these commands on your target machine, not on your build machine. If your target machine does not run Linux, you may have to use different commands.
- If the CPU does not support AVX, you can stop here, and build normally. The CPU does not support the necessary extensions, and will not run code that is compiled using them.
- If the CPU does support AVX, continue. In cmake, enable the FAST_WAVES flag.
- Also in cmake, set SIMD_EXTENSIONS to the appropriate value. If your target CPU supports only AVX, but not FMA, set it to “avx”. If it supports both, set it to “avx+fma”.
- Generate and exit cmake, then build.
- If you made a mistake and used the wrong flags, you should see an illegal instruction error when trying to import CRPropa, or (at the latest) when trying to call
PlaneWaveTurbulence(const TurbulenceSpectrum &spectrum, int Nm = 64, int seed = 0)¶
Create a new instance of PlaneWaveTurbulence with the specified parameters. This generates all of the wavemodes according to the given parameters.
Nm: number of wavemodes that will be used when computing the field. A higher value will give a more accurate representation of the turbulence, but increase the runtime for getField.
seed: can be used to seed the random number generator used to generate the field. This works just like in initTurbulence: a seed of 0 will lead to a randomly initialized RNG.
getField(const Vector3d &pos) const
Evaluates the field at the given position.
Theoretical runtime is O(Nm), where Nm is the number of wavemodes.
const TurbulenceSpectrum &