search

lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
223 stars 292 forks source link

CLASS evolver pk_ref.pre #453

Open moserb94 opened 2 years ago

moserb94 commented 2 years ago

Hi all,

this is mostly a conceptual question. In the CLASS II paper you discuss the different methods for solving the ODEs in CLASS and underline that the ndf15 stiff integrator is more robust and faster than the Runge-Kutta method and thus set as default. But then, in the pk_ref.pre precision setting file, you switch to evolver = 0 (which is the rk solver). Is there a reason for that? In my benchmarks, I actually see that evolver = 0 is an order of magnitude faster than evolver = 1 when all other pk_ref.pre parameters are kept fixed, but I find this counter-intuitive (and in contradiction with the CLASS II paper). I guess this depends on all the other precision parameters, but would like to hear your insight.

Thanks a lot in advance, Beatrice

ThomasTram commented 2 years ago

Hi Beatrice

That is a good question. NDF15 is an implicit solver, so it has to solve an algebraic system of equations at each time-step. This makes it intrinsically slower than an explicit Runge-Kutta method. (Lots of tricks are being employed to reduce the time required for solving this system of course, such as sparse LU decompositions.) However, the step size required by an explicit (Runge-Kutta) method can be severely limited by the stiffness of the problem. At early times, tight-coupling of the photon-baryon fluid leads to such a scenario. It is nice to have NDF15 as the default evolver because it is robust against time-scales, but if we say in the CLASS II paper that it is always faster, we are overselling it!

The Runge-Kutta integrator in CLASS is furthermore not highly optimised -- for a long time I wanted (someone) to add a Runge-Kutta integrator based on the Dormand-Prince pair with dense output functionality. I have suggested this as student projects a number of times, but so far no-one has been interested! :)

I have an implementation here, https://github.com/ThomasTram/LASAGNA_public/blob/cf679266d7c31a01716d157acd9be9fcb7c112e7/tools/evolver_rk45.c#L235 which should be almost plug-and-play....

Cheers, Thomas