Implementation of turbulent magnetic field from Tautz and Dosch 2013

[tfitoussi]

Dear all,

As far as I know CRPropa currently uses the algorithm of Giacalone & Jokipii 1999 to generate turbulent magnetic field. In 2012, R. C. Tautz pointed out that the implementation of a turbulent magnetic field with Giacalone & Jokipii 1999 is problematic.

"For the test case of isotropic turbulence, the three basic requirements for the turbulence generator are:

1. all three turbulent magnetic field components have the same magnitude on average;
2. the probability for the random orientation of the sample Fourier modes is isotropic;
3. the divergence of the turbulent magnetic field must be zero.

As shown recently, this was not possible using previously used approaches because at least one condition was always violated." (Tautz & Dosch 2013)

In Tautz & Dosch 2013 is presented a modified version of Giacalone & Jokipii 1999 that fulfills all three requirements. Moreover, it should start to have a turbulent behavior with only 8 Fourier mode and should not require more than 16 modes since results are almost indistinguishable with higher numbers of modes.

The low number of modes should allow a computation "on the fly" without using a grid. Moreover, this version does not require Fourier transformations.

I started a long time ago to develop the algorithm in this branch (TD13Field in magnetic field folder). I tested the implementation by reproducing the results of section II-C and in particular figure 2 from the article. Here are the pictures and script to produce it.

From this test, it seems that a correct turbulent magnetic field requires at least 64 modes not 8. Moreover, I do not find the same error on the total turbulent magnetic field strength. I guess there is a mistake somewhere in the code but I did not find it.

I have no more time to work on it now but if someone wants to take a look and contribute, he or she is welcome.

Cheers,

Thomas

[lukasmerten]

Hi Thomas

thanks for bringing up this topic. I know that @antoniusf and @LeanderSchlegel have looked into the turbulent quiet recently, maybe they can add to this discussion.

If we are planning to change the turbulence models we might think about combining this with @adundovi 's work on the bendover behavior of the turbulent cascade.

Best, Lukas

[antoniusf]

Hi, thanks for tagging me in!

Yes, we've also been working on a TD13 implementation for some time now (see here). I guess it's a bit unfortunate that we didn't talk earlier; on the other hand, we now have two implementations to compare! I'm actually rather grateful for this, since my grasp of the underlying physics is not the strongest. (I mostly just wrote the code, and made it fast.)

To give a bit of background, we implemented TD13 specifically to avoid grids and the associated problems with interpolation. I was aware that there were similarities between GJ99 and the code in initTurbulence, but I didn't know that they had the same flaws. I'm afraid that I'm out of my depth here.

We'd intended to PR our TD13 implementation once I'd gotten around to cleaning everything up a bit, but as you can tell, this has not happened yet. If you're comparing the implementations, I should also say that we swapped the “bent” spectrum out for a straight power law, to match the behavior of initTurbulence.

I'm not sure how to proceed from here, but I felt it'd be nice to jump in real quick and say what we've been doing. I'm also going to tag @eicheB @reichherzerp since they were part of the effort, and might be able to say more about the physics side of things. (If you wanted to, I know you're super busy!)

Anyway, thanks so much writing, and I hope we can make something nice out of this!

Cheers, Antonius

[adundovi]

@tfitoussi CRPropa doesn't use the same approach of generating turbulent fields as GJ99. In CRPropa, a so-called grid method is employed (turbulence is generated in k-space on a grid, and then inverted to the real space with FFT), while what GJ99, you and @antoniusf are mentioning we can call a plane-wave method. I quote GJ99 to define what the plane-wave method is:

We sum over a large number of plane waves with a spherically symmetric (and randomly distributed) propagation direction, and with random polarizations and phases. (...)

In basic practical considerations, the plane-wave method is limited by computation time, while the grid method by memory space, hence, they are complementary in that sense. However, because how random phases are constructed in the plane-wave method, there can be some theoretical issues. The grid method, as far as I know, doesn't have these issues. I still have to read Tautz's paper to see if there is something further than I know.

For the grid method, all three points are satisfied under the following conditions:

1. all three turbulent magnetic field components have the same magnitude on average;

If the L_max (the biggest scale of turbulence) is several times smaller than L_box

2. the probability for the random orientation of the sample Fourier modes is isotropic;

In the same case as above (or averaging everything over multiple realizations of the field if L_max ~ L_box).

3. the divergence of the turbulent magnetic field must be zero.

Yes, by construction in k-space ( div(B) = 0 => k . B(k) = 0 ).

I would like to see a plane-wave implementation in CRPropa because it would be easier to understand if some effects in results are due to numerics or not (if two complementary methods yield same results, the chance for numerical artefacts is naturally lower). Yet, in both cases, we would like to have as many modes as possible. More modes mean a wider inertial range of turbulence. If you study scatterings where a particle's gyroradius is larger than typical correlation scales, you don't profit from expanding the inertial range, and the mechanism which explains particle movement is well understood (diff. coeff. ~ r_g^2). On the other hand, the inertial range is where interesting phenomena can occur, especially when a regular component of the magnetic field is present... and we are limited both in memory, and computation time by roughly 3 orders of magnitude in the inertial range :-/

[tfitoussi]

Hi everyone,

@adundovi Thank you for the clarification !

@antoniusf Indeed, it would have been good to speak before, but we can compare our implementations. I took a quick look at yours but I had no time to compare precisely the differences. My implementation is clearly not optimized nowadays since I still have some issues apparently. Did you test the physics of your implementation ? At a certain point, I think we can either choose one or the other implementation, or merge both and do PR if everything look ok in term of physics.

As I said, I'm currently working on another topic but I can be involved time to time on it, if you want to work in that direction.

Cheers, Thomas

[antoniusf]

Hi @tfitoussi,

sorry for the late response, I've been away for an extended weekend.

I was mostly not involved with physics testing for the implementation; I do know that tests were made, but I can't tell you how they relate to the tests you did.

I appreciate the offer, but I sadly do not have that much time either! I should really be working on my master's thesis now. What I can do, in the next few days, is comparing the two implementations for mathematical equivalence. If they're fundamentally doing the same thing, merging them should be simple. However, when it comes to any questions regarding the underlying physics, I'm afraid I won't be of much help – insofar as the current issue requires dealing with those, I think that other people will have to answer them. I can still help with the software side of things, though.

Cheers, Antonius

[antoniusf]

Hey Thomas (and everyone else),

this took me some time to get around to, but I've finally managed to compare the two implementations. (For reference, here are the two files: yours, and ours.) If I've done my job correctly, the following should be an exhaustive list of differences on the algebraic level (meaning that I have ignored differences in the implementations as long as the results end up being computed with the same formula).

1. Most importantly, the two implementations use different amplitude distributions in k-space (the G_k). I haven't checked this, but I assume that @tfitoussi uses the distribution specified by the paper (line 49); in our version, we've switched that out for a straight power law (line 161), because we wanted to make it comparable to initTurbulence.
2. Our implementation uses ψ as the polarization vector, whereas @tfitoussi's uses ξ. The paper states that it does not matter which is chosen (and this is in fact easy to see), so this isn't really important.
3. Similarly, one implementation draws φ (the randomly generated azimuthal component) from [–π, π], while the other draws from [0, 2π]. This also does not matter.
4. Lastly, @tfitoussi's implementation uses Δk_i = k_i+1 – k_i, while ours uses Δk_i = k_i – k_i-1. As far as I can tell, this also does not matter since both of these are only a constant factor apart (recall that the k are log-spaced), and will thus get normalized out of the final A_k.

What this means, in conclusion, is that there are no relevant algebraical differences between the implementations, except for the power spectrum. Which of two (or both, or a completely customizable option) should be used, I cannot say. And as I said, there are also some differences in the specific ways that certain things are implemented (I like parts of @tfitoussi's logspace code better than mine, for example).

What happens next is mostly not up to me. I can't offer much help with the physics, and anything that relates to implementation details and code style would likely be discussed in an eventual PR anyways. I could continue to clean up the code as planned, but I'm not sure that this would be a useful thing to do right now.

– Antonius

[tfitoussi]

Dear @antoniusf,

Thank you for the great job of comparison. I had sadly no time to work on it on my side.

From what I can say, the computation of G_k as a powerlaw is a limit case of the general implementation in the article (i.e. k_n>>1 and q=0). For the rest I don't think indeed it changes a lot.

I don't think that cleaning the code is priority here and you have surely others on your side.

I will take look at the impact of these parameters in my one code and try to merge both if there is no major differences in the results.

Cheers, Thomas

[antoniusf]

Dear @tfitoussi,

that sounds great, thank you! Feel free to let me know if any questions come up with regards to my code. Also, I hadn't realized that the straight power-law was a limit of the general formulation, so that's good.

Thanks, Antonius

[tfitoussi]

Hi @antoniusf,

I made a fork of your branch here. I cleaned up a little the tree to merge the last update from the master branch of CRPropa and I made a few changes so that I can make comparison.

1. I replaced the definition with a coherence length by a minimal and maximal length like in initTurbulence
2. I added an option to choose between a pure Gk powerlaw and Gk as define in the paper (default)
3. I kept the SLEEF parts even if I don't use it and don't know it.

It can be use as follows:

Bfield = crpropa.TD13Field(Brms, Lmin, Lmax, s, q, Nmode, 0, True) 

True: set the powerlaw instead of the article definition of Gk.

I run the same tests I did with mine. The results are here with and without the powerlaw. On each graphic I put the result in the same condition with the initTurbulence grid.

From what I see, the results looks much better than with my code, especially with a low number of modes. There a turbulent behavior with only 16 modes and is equivalent to the version with the initTurbulence method.

If it is ok for everyone, we can make a pull request for the implementation.

Thomas

[antoniusf]

Hi @tfitoussi,

thank you so much for all your work. I feel a bit bad responding as slowly as I do, given how much you've done already.

I feel even worse that didn't tell you how terribly organized the repo is. Counterintuitively, the TD13 branch that you merged in does not actually contain our current implementation, because that sits in the field-perf branch (I was working on making it faster, hence the name). I'm really sorry for this, there is no way you could have anticipated this problem. I should have stated the branch name explicitly and not just linked to it.

However, here's the good news: There aren't actually that many differences! Firstly, the current version uses kmin and kmax instead of that weird correlation length thing, and you've already fixed that. Secondly, the old code used a weird normalization method (the current one is back to doing it properly). I think this was a test that we did at the time, but this probably shouldn't be in the final version (it wasn't in yours, either). So this should be removed, but it's only like two lines of code, and one big comment that can be deleted. And lastly, the current version has a very different optimized implementation (what corresponds to the SLEEF parts of the old version), but you haven't touched this part at all! So it should be simple to transplant in.

I don't want all of the work you've put into the merge to be wasted, and I think the simplest way to do this is to just fix the problems directly, no weird git magic. Maybe split up into two commits, one for fixing the normalization, and one for the new optimized code?

I'm offering to fix this, since I feel like it's kind of my fault. What do you think?

Sorry again, Antonius

[tfitoussi]

Hi @antoniusf,

Don't feel bad, I didn't answer for three weeks before '^^. It was clear since the beginning that none of us really had a lot of time to work on it, so we contribute time to time and it's ok for me. This week, I have more time since I'm waiting some other results to be computed, so I can answer quickly.

After I'm not a big guru with git neither. It's the first time I use it to collaborate with someone. The rest of the time, I'm the only developer and the only one to complain of my dirty stuff. So I would also prefer to avoid black magic git command except if it 100% transparent to me.

Then, I think it would be better to work on the same branch. Perhaps, you can clone mine and put the modifications you think necessary. When it's done, you tell me and I can pull it and so on. That way, I think we can avoid having 3 branches with 3 levels of development. When the code is clean enough then one or the other can ask for a pull request. After if you want to improve the code, you can keep working on it and ask a new PR.

Tell me if it's ok for you. Thomas

[antoniusf]

Hi @tfitoussi,

Thanks, that's nice to hear.

I think your approach sounds great! I'll pull down your branch and do the necessary modifications, then notify you when that's done.

– Antonius

[antoniusf]

Hi @tfitoussi!

Good news: Thanks to the CRPropa busy week, I finally got around to doing this. :) My branch is here, but I've also opened a PR against your TD13 branch, so you should be able to pull in my changes easily. I've checked the results of the merged version against my old code and found them to be identical. (Not exactly, since the small change in k generation has knock-on effects, but I think that's okay. When I change the k generation back, the results are bit-by-bit identical, which is a strong signal that nothing has gone wrong.)

I've also run clang-format with CRPropa's style file and cleaned up some old comments of mine.

Oh, @adundovi: It says here that you assigned this issue to yourself. What does this mean?

[tfitoussi]

Hi @antoniusf,

Great job ! I quickly checked the code and I still find the same results. So for me, it is OK for a PR on the master branch of CRPropa. Do you want to do it or I go ?

Cheers

[adundovi]

@antoniusf Nothing special, just that I'll try to take more responsibility to coordinate this issue as one of CRPropa developers.

First, I'm curious about the compatibility of the faster implementation with recent compilers?

[antoniusf]

@tfitoussi I'm not sure, what does @adundovi say? Also, if one of us opens a PR, can the other still add commits?

@adundovi Oh, okay, thank you so much! So, as for compiler compatibility, I think the only thing in the code that needs C++11 is the call to std::align. However, I've recently learned how to do this manually, so we should be able to remove this dependency and get the same result.¹ And of course, the compiler will have to provide AVX intrinsics, which are about a decade old by now. This should be independent of C++ version though, I think. Does this answer your question? I'm not sure if I understood it correctly.

¹ replacing std::align(32, 32, pointer, size) with (pointer + 32 - 1) & ~(32 - 1)

[antoniusf]

Also on the topic of AVX: I've also written a version that uses AVX's predecessor, SSE. While there are still some old processors in use that don't support AVX, SSE is old enough that it's probably supported by everything people use for simulations. However, that version is about 2x slower than the AVX version. (Though this is still a lot faster than non-optimized.) The code is very similar to the AVX version, but different enough that I don't see a way to merge the two. So including both would result in a fair amount of seemingly duplicate code, and I'm not sure this is worth it, given that AVX is pretty widespread, too.

[TobiasWinchen]

Please use C++11 align, as C++11 is now mandatory for CRPropa anyway

[antoniusf]

Okay, great!

[adundovi]

Yes, C++11 is now required in CRPropa. Thanks for explaining my compatibility question. I'll try to test this implementation myself, however, let see if we can generalize this to put it side by side to the grid turbulence functions since they provide the same logic. I'll think more about it, and return here with some better ideas.

[antoniusf]

Okay! I'm curious what you come up with, because right now I'm not quite sure I understand.

I think you could conceptualize grid computation as two steps: First, generation of the wavemode parameters, and second, evaluation of these wavemodes at certain points. In TD13, evaluation happens in getField for a single point. In initTurbulence, evaluation happens for an entire grid of points at once, in the call to executeInverseFFTInplace. The two functions do things that are very similar mathematically, but the implementations are optimized for very different sets of constraints. (Doing one thing in as little time as possible vs doing lots of things as efficiently as possible.)

This would leave generation as a place to generalize. I don't know much about the physical details (though TD13 seem to have put a fair bit of thought into those) but from a pure implementation perspective, there are differences here, too. initTurbulence works strictly on a grid, while TD13 distributes wavenumbers in logspace, which means that they have to apply different scaling factors to each wavemode to achieve the same effective densities. (In TD13, the density of discrete wavemodes decreases significantly for high wavenumbers, so the amplitudes are scaled up accordingly to compensate.) Also, initTurbulence currently uses complex numbers, while TD13 only uses reals (and adds on the phase separately). iirc, in TD13 this is important for performance, since evaluating a cosine is faster for real than for complex numbers.

At least, this is my view of the situation. Let me know what you think, or what your ideas were. I hope I wasn't too obscure in what I wrote.

[adundovi]

No, no. Of course, I understand this. :-)

I'm just saying that the grid version is placed outside of the magnetic field folder as in principle can be used elsewhere (not only for magnetic fields) - but at the moment I don't see if a vector field is of any use in the context of cosmic rays and CRPropa besides in the magnetic fields. Therefore, I would like to put those two methods side by side in the CRPropa code structure, either by putting the grid method in the magnetic field folder (preferably) or generalizing the wave method so that it is not directly seen as magnetic field and putting it outside. This is the first reason, more technical.

The second reason is that the spectrum of turbulence should be generalized in both methods systematically because the one used now in the current default grid method is questionable anyway as it doesn't respect the homogeneity requirement (see for example Batchelor's book on Homogeneous turbulence, chapter 3 - all theory we use in these synthetic models is based on those foundations) but introduce a sharp cut-off at L_max. I see that the T&D prescription respects it, so this is good. However, L_c in that case doesn't look correct as it depends also on q (the large scale index), not only s (the spectrum index). I'll give you the correct formula for L_c, although it can be easily calculated from the spectrum. The other related but minor reason is that we have a code duplication when calculating L_c. This should be shared between the two methods.

Let me study a bit your implementation for a day or two, and maybe I can comment something more and better than this.

[antoniusf]

Ooh, okay, thank you for explaining, that sounds like a really good idea! Now I understand what you meant by "generalize", and I agree that doing this for the spectrum (and L_c) sounds smart. I'm sorry for making assumptions, I should have waited for you to elaborate first.

Okay, I'm curious to hear what you come up with!

[adundovi]

Sorry for the delay. I have finally looked into the TD13 paper and your implementation of it in CRPropa, and I have several questions...

First and foremost, did anyone try to compare TD13 and the default one? From what I can tell, they cannot be easily compared because they use different spectra what causes different turbulent field. In specific, the comment says:

Usually, you'd want to use s=5/3 and q=0 here, which will be comparable to
     passing -11/3 to `initTurbulence`.

which will not lead to the same spectra since the default one is just a power law k^-s, while TD13 has (1+k)^-s. That makes a difference at small wavenumbers. To be even more specific:

Gk = pow(k, q) / pow(1 + k * k, (s + q) / 2.);

is not the same as:

b *= random.randNorm() * pow(k, alpha / 2);

even when you put q=0. Eh, now I see that there is the undocumented powerlaw parameter which uses the default definition. I stand corrected, however, this option should be dropped, and only the new spectra should be adopted.

Furthermore, since k^q / (1+k)^(s) is already there (which is just better than k^-s as I explained already in two comments earlier), it should be used with q=4 at least, and for 3D turbulence, this value should be the default; q=0 makes sense only for 1D turbulence, i.e., the slab model what was displayed in the TD13 paper as I can notice, but for the 3D case it is even more wrong than k^-s of the default model as it gives a crazy coherence length, which probably does not make sense.

As far as I can see, this implementation of TD13 misses an important parameter called the bend-over scale (labelled as l0 = 0.03 a. u. in the paper) which normalizes the wavenumber in the spectrum G(k). It is a number of dimension of length (L) and it determines the position of the energy scale. Implicitly, if I read correctly the code, it is set to 1 - which will not work if the other scales in the simulation are not comparable to it. It should be a user-definable parameter available in the constructor.

I'll look further into the implementation and the paper since they could be other details that I missed in the first read.

[adundovi]

I should have read this thread again from the beginning before posting my last comment. My first and second questions are already answered - I see that the spectrum issue is a by-product of merging two implementations, where the comment in the code left from only one of them. Moreover, I see that @tfitoussi performed a comparison test, yet, due to my last point (of the implicit bend-over scale l0=1), one doesn't see much of a difference between the two, which shouldn't be the case in a regular situation (when l0 ~ system size). Also, the current L_c works only for k^(-s), but not for k^q/(1+k)^-s...

The coherence length should be: L_c = \frac{4\pi}{s(s+2)} C(q, s) l_0 where C(q,s) = \Gamma((s+q)/2) / (2 \Gamma((s-1)/2) \Gamma((q+1)/2). Although, maybe ~4\pi can change depending on the taken spectrum normalization. It is obtained after integrating the correlation tensor R_ij=<B_i(x)B_j(x+r)> from 0 to \infty in r. However, when we have a finite spectrum, from k_min to k_max, and if l ~ 1, this formula will not work, and then the following crazy expression should be evaluated (with hyper-geometric functions :-/ ):

def corr_length_isotropic_finite(
    n: "spectral index", Lbo: "bendover scale", Lmax: "L max scale", Lmin: "L min scale"
) -> float:
    return (
        (2 * np.pi) ** n
        * (n - 1)
        / (2 * n * (n + 2))
        * (
            Lmax ** n
            / Lbo ** (n - 1)
            * (2 * np.pi ** 2 * (n + 2) + Lmax ** 2 / Lbo ** 2)
            / (4 * np.pi ** 2 + Lmax ** 2 / Lbo ** 2) ** (n / 2 + 1)
            - Lmin ** n
            / Lbo ** (n - 1)
            * (2 * np.pi ** 2 * (n + 2) + Lmin ** 2 / Lbo ** 2)
            / (4 * np.pi ** 2 + Lmin ** 2 / Lbo ** 2) ** (n / 2 + 1)
        )
        / (
            (Lmax / Lbo) ** (n - 1)
            * scipy.special.hyp2f1(
                (n - 1) / 2.0,
                (n + 4) / 2.0,
                (n + 1) / 2.0,
                -Lmax ** 2 / (4 * np.pi ** 2 * Lbo ** 2),
            )
            - (Lmin / Lbo) ** (n - 1)
            * scipy.special.hyp2f1(
                (n - 1) / 2.0,
                (n + 4) / 2.0,
                (n + 1) / 2.0,
                -Lmin ** 2 / (4 * np.pi ** 2 * Lbo ** 2),
            )
        )
    )

[adundovi]

I've opened a new PR, #287, which is not meant to be merged, but to serve as tracker for this model also. Please rebase your branch against it, put the TD13Field.h inside of turbulentField, and inherit the abstract class and use the generic spectrum provided by that class. For the moment, don't provide the Lc feature as we would have to think which one actually to implement.

[antoniusf]

Hi,

thanks for all of your effort so far! I'm sorry if the merge left some confusing parts, I must have missed that.

I can take care of rebasing and adapting TD13 to turbulentField in the next few days, if that's desired. I'm not sure if there are still open questions I can help answer, but most of the spectrum related stuff should be covered by the new base class now anyways, right?

– Antonius

[adundovi]

Hi Antonius,

that's great! Thanks. Yes, once the spectrum is inherited from the base class my concerns, at least about the spectrum, should be resolved! I have prepared several test-particle simulation checks which should give the same results between the grid and the plane-wave method (some of them also have analytical results). Before performed tests by @tfitoussi (as far as I saw from the attached figures) only checked if the field is divergence-free and if it reproduces the energy spectrum, but I'm curious what happens to the random phases since they can carry information about the turbulence structure that shouldn't be present in synthetic fields like these. But if there are some correlations between the phases in a way how you generate plane-waves, these may mimic a structure and affect test-particle simulations which would look more like done on an MHD than a synthetic field... we'll see, but I hope not :-)

One more thing, about the naming. Please, rename your class to PlaneWaveTurbulence but keep the reference to the TD13 paper in the description, of course. Plane-wave is something more recognisable by the community than TD13, and it should be seen as a complementary method to GridTurbulence that I now introduced in #287 . Therefore, they will be placed on an equal footing in the code... and that was my main intention.

[antoniusf]

Okay, will do! One more question: You said to rebase the current TD13 branch onto your PR. Would a merge do as well? The TD13 does not have a straightforward git history, since it contains modifications by both me and @tfitoussi, so I feel like a rebase might be somewhat involved.

Also on the topic of the git history: It was bugging me a bit that the history of our work on TD13 wasn't really reflected in git, so I've gone in and cleaned this up a bit. There's now a nice simple merge commit that cleanly merges both of our TD13 histories. The code is the exact same as the one that Thomas and I have produced in the past few weeks, it's just that this is now also reflected in the git history. The branch is here. Let me know if you all agree with this decision or if you would just like me to continue with the current TD13 branch. (Again, the files are the same, only the history is cleaner.)

[adundovi]

Hi @antoniusf, this will work too. Actually, I have already merged your td13-clean branch with the general turbulence PR #287, so you can continue from there. Once TD13 is "converted" to the TurbulenceField base class, and after I clean a bit other files and add some tests for all of these, we will proceed with the merge of #287 into the main branch.

[antoniusf]

Awesome, thank you!

[antoniusf]

One more question: which convention would you like to be used for the spectral index? As it is used in the TD13 paper, with 5/3 for a Kolmogorov-like spectrum? Or as it is used in initTurbulence, with –11/3? (Brief explanation for those wondering, since it took me way too long to understand this: TD13 uses a –5/3 spectrum, since the outer wavemodes have larger volumes to fill (proportional to k²). This means that when the individual amplitudes are distributed as k^-5/3, the effective density (I guess?) evaluates to k^-5/3 · k^-2 = k^–11/3, same as initTurbulence.) Or is this something I should not concern myself right now, and just assume that everything works out?

[antoniusf]

Please let me know if this covers your expecations or if there's something you'd still like me to do.

[adundovi]

I'm for 5/3, and I've already put it in the new base class (called sindex). Firstly, it is the original Komlogorov result widely recognized. Secondly, the additional k^-2 in -11/3 comes from the correction for the 3D geometry (to cancel out k^2 which comes when the integration over wave-numbers is performed: dk^3 is transformed to spherical coordinates k^2 dk), while in the slab or 2D geometry it is different (~1 and ~k, respectively), which can lead to confusion if you don't specify details. Here is the proper definition of the omnidirectional energy spectrum and it is universal for all geometries:

[adundovi]

Please let me know if this covers your expecations or if there's something you'd still like me to do.

It is fine as far as I can see. Great work! I'll see now how it works in practice, and improve the other components of this big PR (TurbulentField and GridTurbulence) so that it can be merged with the master. BTW, q = 4 should be the default value, not 0.

[antoniusf]

Oh, okay, thank you for the explanation! I wasn't aware that 5/3 was the standard value and 11/3 the corrected one. This now makes sense to me, thanks!

I must have missed the q = 4 thing, sorry. I'm assuming that you're going to fix this during your cleanup, or would it be helpful to you if I commited the fix?

[adundovi]

Hi @antoniusf, I've merged your changes in the PR. So far, looks good. I've removed L_c from your class and implemented the correct one inside of TurbulentField (luckily C++11 has the gamma function in the standard library), therefore, it works for all (except the Simple model which has its own L_c). I'm still undecided regarding the arguments' order in the constructors: should those regarding physics be first, or those technical better; should we provide default values at all or force users to check them before using classes... also, maybe there are too many of them :-/ I'll start now adding tests and checking functionality of both approaches, but that will take me several days. I'll add also more in-code comments. Then I'll ask others for the code review so that we can merge this PR.

[antoniusf]

Hi @adundovi,

that sounds great, thank you!

The order of arguments is also something I wasn't able to solve to my satisfaction, especially considering default values and such. Speaking from recent experience with Rust library conventions, it might be worth to consider the builder pattern? e.g.

PlaneWaveTurbulenceBuilder()
      .Brms(1*nG)
      .Lmin(...)
      ...
      .build()

(the build method then checks that all necessesary values are present, fills in the defaults, and constructs the actual object.)

Alternatively, maybe there's a middle road where we separate only the physical and technical aspects? Perhaps the spectrum could be separated out into its own class, taking the energySpectrum method with it?

TurbulenceSpectrum spectrum(1*nG, lmin, lmax, 5/3.); // use defaults for q and l_bendover
PlaneWaveTurbulence field(spectrum, 100, 0); // 100 modes, automatic seed

PlaneWaveTurbulence could then do something like Gk = spectrum.get_amplitude(k). Anyway, feel free to let me know what you think!

I'm looking forward to the review and the merge. Thanks again for taking on this issue and for all the work you're doing on it.

[adundovi]

Although I like the Rust way, I don't think that people of this community would generally recognise such a pattern (save that idea for later when we'll rewrite CRPropa in Rust ;-) ). I'm also not sure how it would work with SWIG which is quite limited when used in non-default C++ design patterns. I already started drafting something similar what you suggested as the alternative, yet, we have to check if this will impact the performance. If it will not have an impact, that is the way to go. Regarding separating technical from physical arguments, I'm not sure whether L_min and L_max are physical or technical parameters. In the ideal case of homogeneous turbulence, k_min = 0 and k_max = infty, therefore, L_min and L_max are more like approximated finite limits, i.e., technical parameters :-)

But before we continue further with UX decisions, I still have to convince myself if PlaneWaveTurbulence really works :-D

[antoniusf]

Okay, that makes sense to me. Awesome!

Huh, I see what you mean with regards to L_min and L_max. I'm curious to see what you will come up with.

I really do hope that integrating PlaneWaveTurbulence didn't break anything, but I think you were also referring to the actual physics. But please just keep in mind that I have not tested the integrated version of the code against the older versions.

[adundovi]

By convincing I meant both as an implementation and as an approach :-) So far I've been "specializing" the grid approach and I understand well what results produced by it makes physical sense and what are numerical artefacts. With the plane-wave approach, I don't know yet. I also don't understand how the number of modes affects the resonance of the scattering. TD13 is not clear on that, only suggesting that 8 (eight) modes are somehow fine (!) with the slab turbulence which makes no much sense to me. GC99 is even less clear. In TD13, figure 6 shows D_\perp (or "lambda") which is sub-diffusive in the slab case (well-known result) and so it is hard to argue how well the method works in the sub-diffusive regime as no theory covers it, thus saying anything about the number of modes, in this case, seems like guessing. Fig 7 is a better, however, would still like to see how it saturates in different regimes. It seems I have to try it myself.

The first test, comparison with the isotropic turbulence with no background field works really good - better than the grid approach as I can see.

[TobiasWinchen]

The build pattern should probably work in python and swigged objects, but unless we introduce it on a large scale in CRPropa we should not use it here: We should avoid having to many different interface patterns in different code parts to keep complexity low. In particular, as python kwargs solve the problem already for most of CRPropa steering already, so only the C++ interface would benefit here.

[adundovi]

The build pattern should probably work in python and swigged objects, but unless we introduce it on a large scale in CRPropa we should not use it here: We should avoid having to many different interface patterns in different code parts to keep complexity low. In particular, as python kwargs solve the problem already for most of CRPropa steering already, so only the C++ interface would benefit here.

Are you referring to the Rust-like pattern or providing the spectrum parameters through a new object like @antoniusf suggested as alternative? Do you think that we should put all 8+ parameters directly in the turbulent field constructor? If not, what approach you would use?

[TobiasWinchen]

I don't like the builder object, but you could use a static method as in https://en.wikibooks.org/wiki/More_C%2B%2B_Idioms/Named_Parameter which is nice for this purpose. We do not use this in CRPropa yet, but we could. However, python kwargs are preferable and easily automatically created having the parameters in the constructor, so the pattern helps only in constructing C++ objects and not in the order of parameters. Note that for the documentation you have to think about the order of parameters anyway.

[adundovi]

@TobiasWinchen The issue I have is that I would like to define the order of parameters in the constructor which follow the importance of parameters in physical sense, i.e., when you define a turbulent magnetic field, first you are interested in Brms, following by the spectrum of the inertial range, the L_c (which is directly related to the bendover scale)... finally, you reach "technical" parameters: the size of the box, and the grid (or number of modes), Lmin, Lmax, random seed, etc. However, for some parameters default values can be prescribed (like s=5/3), while for others they cannot. Neither Python nor C++ allow putting "prescribed" parameters before "free" parameters which is dictating me the order different than what I want. If I don't put any default value, I have all the freedom, yet too many parameters has to be specified (even the "common" ones) which can be a hurdle for users.

[TobiasWinchen]

OK - For the case here introducing a TurbulenceSpectrum struct with appropriate constructor is probably the easiest way to go, for the build pattern I would like to see a complete draft first to get a better picture, and we should think on how this would look in other parts of the code, so this would be more work to make it nice with probably only little benefit (KISS). Also, for python you will want kwargs so this will be more work, whereas the struct will allow

field = PlaneWaveTurbulence(TurbulenceSpectrum(1*nG, lmin, lmax, 5/3.), modes = 100, seed = 0)

automatically.

[adundovi]

OK, the test with a regular background field also went very well - matches results of the grid method. However, the performance of the evaluation dropped considerably - I'm not sure why since the PlaneWaveTurbulence instance stayed the same.

I'll clean up the interface of all classes now and request a code review...

[antoniusf]

Huh, the drop in performance seems weird, I could take a look at that if you want? By about what factor did it get slower?

[adundovi]

My performance issue is not something that blocks the merge - don't worry. :-) The implementation works great so far. Actually, looks superior to the grid method, although I still don't understand what is the proper (?) number of modes (Nm) that one has to adopt for her/his problem.

Regarding the performance issue. I turned all optimizations on, but the difference is the same with or without them just to make clear :-)

I have something like this:

BField = MagneticFieldList()
BField.addField(
        PlaneWaveTurbulence(Brms, s, q, bendover, lMin, lMax, Nm, seed)
    )
if B0 != 0:
    BField.addField(
         UniformMagneticField(Vector3d(0,0,B0))
    )

When B0 is 0 the simulation is incredible fast, but when I increase B0, it becomes slower and slower... even the order of magnitude slower when B0 > Brms. I use PropagationBP. I don't remember having the same difference with the grid method. How can it be?

[antoniusf]

Hey, sorry for the delayed response. I'm relieved that this doesn't seem to be the optimization's fault, I was a bit worried about that. Do you think that this might have something to do with the physics? For example, that it somehow has to do more steps per particle as B increases? (I still have some stepcount instrumentation lying around from when I was doing perf tests, I might be able to dig that back up...)

I'm very happy to hear that it's not worse than the grid, though I also really don't have an idea about the number of required wavemodes. (Oh, and a small suggestion: if you would like, you could try using “their” instead “her/his”, as a form of the singular they. It's a bit shorter and personally I think it's easier to read, plus there might be some other advantages. But you really don't have to, I know people can get a bit annoyed when this comes up, I just wanted to mention this in case you hadn't heard of it yet.)