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)
230 stars 285 forks source link

Are there any changes required to the python wrapper when changing perturbations? #418

Closed kabeleh closed 3 years ago

kabeleh commented 3 years ago

I am trying to run a modified version of CLASS 2.9.4 together with MontePython. Just calling CLASS works fine, but MontePython fails. Now I am wondering, if there might be some changes required to the python wrapper? I've seen similar issues already around (#390 , MP#183, MP#198) yet no viable solution. I am on python 3.9.3 and I am using CLASS 2.9.4 with modified input.c, background.c, and perturbations.c (with corresponding changes in background.h and perturbations.h, as well as classy.c, classy.pyx, and cclassy.pyxd). All changes are related to the scalar field or dark matter. Everything runs fine with the default perturbations, but not with the modified perturbations. CLASS itself runs with the modified perturbations, but when I try using it in MontePython, it throws between 3 and 10 times the following error message (with different interval: [X:Y] values)

Error in Class: perturb_init(L:841) :error in perturb_solve(ppr, pba, pth, ppt, index_md, index_ic, index_k, pppw[thread]);
=>perturb_solve(L:3096) :error in generic_evolver(perturb_derivs, interval_limit[index_interval], interval_limit[index_interval+1], ppw->pv->y, ppw->pv->used_in_sources, ppw->pv->pt_size, &ppaw, ppr->tol_perturb_integration, ppr->smallest_allowed_variation, perturb_timescale, ppr->perturb_integration_stepsize, ppt->tau_sampling, tau_actual_size, perturb_sources, perhaps_print_variables, ppt->error_message);
=>evolver_ndf15(L:462) :condition (absh <= hmin) is true; Step size too small: step:2.22045e-16, minimum:2.22045e-16, in interval: [1645.87:14213.7]

before it aborts with

Traceback (most recent call last):
  File "/home/MontePython/montepython_public-3.4/montepython/MontePython.py", line 40, in <module>
    sys.exit(run())
  File "/home/MontePython/montepython_public-3.4/montepython/run.py", line 45, in run
    sampler.run(cosmo, data, command_line)
  File "/home//MontePython/montepython_public-3.4/montepython/sampler.py", line 46, in run
    mcmc.chain(cosmo, data, command_line)
  File "/home/MontePython/montepython_public-3.4/montepython/mcmc.py", line 787, in chain
    newloglike = sampler.compute_lkl(cosmo, data)
  File "/home/MontePython/montepython_public-3.4/montepython/sampler.py", line 776, in compute_lkl
    value = likelihood.loglkl(cosmo, data)
  File "/home/MontePython/montepython_public-3.4/montepython/likelihood_class.py", line 952, in loglkl
    cl = self.get_cl(cosmo)
  File "/home/MontePython/montepython_public-3.4/montepython/likelihood_class.py", line 191, in get_cl
    cl = cosmo.lensed_cl(int(l_max))
  File "classy.pyx", line 567, in classy.Class.lensed_cl
classy.CosmoSevereError: 

Error in Class: No lensed Cl computed

I played around with the parameter ranges in MontePython and I also tried fixing the scalar field parameters in MontePython to values that work with CLASS when just calling CLASS through an *.ini.file, i.e. setting those to constant ones in the MontePython parameter file. I also tried including data.cosmo_arguments['lensing'] = 'yes' and data.cosmo_arguments['output'] = 'tCl,pCl,lCl,mPk' in MontePython without success.

Is there anything else I could try or something I should check? Do I need to make some specific changes to the python wrapper, when changing something in perturbations.c?

Thank you very much in advance. Any help is highly appreciated.

kabeleh commented 3 years ago

I realised that the results for the matter power specturm P(k) and cl(lensed) are either very high, 0.00000 or -nan. CLASS itself did not complain about this, but MontePython aborted. So it is rather an issue with my own code and I apologise for opening this thread. I didn't realise this earlier, because I wasn't aware of the fact that MontePython automatically askes for 'output': 'mPk tCl lCl pCl ', even if it is not specified by the user.

kabeleh commented 3 years ago

I can now happily close this issue. It was not actually an issue with my code per se. But the perturbed equations of motion were in a way that weren't handled well by the implicit solver. Setting evolver=0 in the *.ini-file solved my problem. The spectra look good now and MontePython is running without further hiccups (so far).