CRPropa / CRPropa3

CRPropa is a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles. https://crpropa.github.io/CRPropa3/
https://crpropa.desy.de
GNU General Public License v3.0
68 stars 67 forks source link

Uniform interface for the turbulence generation #388

Closed LeanderSchlegel closed 2 years ago

LeanderSchlegel commented 2 years ago

Dear all,

with this pullrequest we want to provide a fix to the issue of the non-uniform spectral index handling, as described in #387. Analyzing the problem, we made two changes:

First, to solve the problem for the PlaneWaveTurbulence method using SimpleTurbulenceSpectrum that was shown in #387, we effectively removed the "+1" in the exponent in the denominator of the formula of the power spectrum, to be accordant to the Tautz and Dosch 2013 paper, as suggested in the discussion of the issue.

Second, we found that for a uniform handling, it would be good to also change the argument that is given to energySpectrum() in the general case of TurbulenceSpectrum. Instead of giving it k^ = k * l_bendover, we put the multiplication with the bendoverscale inside of the function to have the same argument (just k) as in SimpleTurbulenceSpectrum.

We tested the combinations of the turbulence models (SimpleTurbulenceSpectrum and TurbulenceSpectrum) with the available methods (PlaneWaveTurbulence and gridbased GridTurbulence and SimpleGridTurbulence) and found them working as expected, producing the right spectral behaviour, as shown in the plots below, where the spectrum for s=5/3. is shown reweighted with 11/3. On the left the PlaneWaveMethod is shown and on the right the gridbased methods.

fixed

While testing the FastWaves optimization for the PlaneWaveTurbulence method, we found a minor typo preventing it from being built, that we fixed within this PR.

Thank you in advance, @JulienDoerner , @LeanderSchlegel

antoniusf commented 2 years ago

Hiya,

just a short update from me: I've been unable to work the past few weeks due to wrist problems, but I will hopefully get around to reviewing this soon.

A quick note of confusion from me: Why is my name undersigned on the PR comment? I've not had any part in this PR up until now.

Cheers, Antonius

LeanderSchlegel commented 2 years ago

Hi,

no problem, I hope it got better! Thanks for reviewing it. We wanted to credit your analysis tool and suggestion of the missing "+1" in the issue:)

Cheers, Leander

antoniusf commented 2 years ago

I see, thanks for explaining!

antoniusf commented 2 years ago

Hey, just so we're on the same page on this (because in the mess of the past two months I don't think I've ever communicated this publicly) – the current PR doesn't implement the exact fix suggested in #387 (i.e. fixing the spectrum generation in TurbulenceSpectrum::energySpectrum), but instead fixes the generated spectrum in PlaneWaveTurbulence, is that correct? If I see this correctly, the resulting spectrum only mostly follows the formula from Tautz & Dosch, since it now outputs

which for k >> 1 (or rather k >> l_bendover) is basically the same as TD13's formula

but since they're not exactly equal (and the approximate equality is only valid for large k), I thought it was worth mentioning this.

Also, @LeanderSchlegel, would you mind removing my name from the PR comment? I don't think it's bad that you want to credit me for some tooling I've written a few years ago, but to me it reads as if I've been involved in the PR somehow and I feel like that might be confusing to others too.

I'm sorry for the delays my injury (and general life circumstances) has caused with respect to my response to this issue; but I hope it's still valuable to you.

JulienDoerner commented 2 years ago

hey,

I corrected the difference @antoniusf reported. Now the correction term is 1 + k^2. Therefore it is necessary to introduce a bendoverscale in the simpleTurbulenceSpectrum. In this case I set it to 1000 times the maximal turbulence scale to ensure to be in the limit k^2 >> 1. I checked all tests mentioned before and there is no change in the expected behavior.

Additionally I updated the documentation as @reichherzerp requested.

adundovi commented 2 years ago

Hi @JulienDoerner, thank you for the hard work you put in to resolve this issue I unintentionally caused by merging these two approaches of generating turbulence. I'll checkout the PR now locally and run a few tests myself. However, putting a number, 1000, is a bit confusing - I would either like to avoid it or define a constant equals to 1000 which is named according to some idea from where this 1000 comes from.

JulienDoerner commented 2 years ago

Hey @adundovi, I introduced a variable for the scaling factor of 1000 and documented the meaning.

antoniusf commented 2 years ago

I haven't tested it, but everything looks good to me too. Thank you so much for your work!