diff --git a/README b/README new file mode 100644 index 0000000..4dbc6db --- /dev/null +++ b/README @@ -0,0 +1,94 @@ +********************************************************************** +This is an addon module for CAMB/cosmoMC which computes the +primordial power spectra by an exact numerical integration of some +standard one field inflationary models. It implements reheating +through the "reheating parameter" R allowing robust inflationary +parameter estimations and inference on the reheating energy scale. The +underlying perturbation code actually deals with N fields +minimally-coupled and/or non-minimally coupled to gravity. Works for +flat FLRW only. +********************************************************************** + +-These patches modify the cosmomc or camb codes and therefore may have +bugs that are not present in the originals. The authors of the +original versions should not be bothered with problems resulting in +the use of these patches. Feel free to insult me, or at least contact +me, in case of bugs! + +-Apply the patch in the cosmomc/ directory: + +patch -p1 < purpose_MMYY.patch + +-See astro-ph/0509727, astro-ph/0605367, astro-ph/0703486, arXiv:1004.5525 + +-Please note that the genuine Makefiles have been heavily changed and +may require edition to work on your system. They assume the existence +of a directory "cosmomc/WMAP" which points to the WMAP team likelihood +code directory. The WMAP data are also assumed to reside in the +"cosmomc/data" directory. + +********************************************************************** + +Basic controls are included in the cosmomc "params.ini" file and +consist on choosing an inflationary model and specifying the prior +range on its potential parameters. The CAMB control file, of the same +name, includes a detailed interface to the inflationary module. + +The inflation code itself is in the sub-directory "fieldinf" and a +short description of the relevant files is given in the following. + +*infbgmodel.f90 + +Contains the encoded inflationary potentials together with their +short-name. If you need to add a model, this starts there. + +*infsrmodel.F90 + +Contains routines to solve the slow-roll approximated equations for +each of the above mentioned model. This is used only to provide a +guess (optional) of the initial field values producing a numerically +optimal *total* number of e-folds of inflation (long enough to perform +computations on the attractor, but not too long to spare computing +time). Additional "hard" theoretical priors on the allowed field +values are coded in this file. + +*infbg.f90 + +Integrates the field and metric background from given values of the +potential parameters and initial conditions. The function +"set_infbg_ini" is precisely calling the routines of "infsrmodel.f90", +but only if you set all initial field value to 0. Otherwise, the +background is computed from the input initial field value. A lot of +tunable options can be set in the function "bg_field_evol" to specify +how inflation ends, to explore parametric oscillations... + +*infpert.f90 + +Integrates the perturbations in Fourier space over the previously +computed background. For each perturbation mode, you get the tensor +and scalar power spectra. + +*inftorad.f90 + +Contains all routines needed to relate the inflationary era to the +radiation era by means of the reheating parameter R. Extra assumptions +on reheating can be fixed there. + + +*power_inf.f90 + +This is the interface between CAMB and the inflationary code. It feeds +CAMB with each perturbation mode. It contains various optimisations, +tests for hard priors. A spline of the power spectra is switched on +by the boolean "useSplineDefault". The range and number of points are +set in the routine "SetDefPowerParams". + +********************************************************************** + +All source files contain the switches "display" and "dump_files" which +can be used for debugging. They output data at various stages of the +computation. Mind that they slow down a lot the execution and you +don't want to use them for production runs. + + + diff --git a/camb/Makefile b/camb/Makefile index 5c09a67..01ca343 100644 --- a/camb/Makefile +++ b/camb/Makefile @@ -1,83 +1,134 @@ -#CAMB Makefile +# >>> DESIGNED FOR GMAKE <<< +#CAMB System unified Makefile -#Set FISHER=Y to compile bispectrum fisher matrix code -FISHER= +ext=$(shell uname | cut -c1-3) -#Edit for your compiler -#Note there are many old ifc versions, some of which behave oddly +ifeq ($(ext),IRI) +F90C=f90 +FFLAGS = -Ofast=ip35 -n32 +LFLAGS= +endif + +ifeq ($(ext),Lin) +F90C=gfortran +FFLAGS= -O -fopenmp +LFLAGS= +endif + +ifeq ($(ext),OSF) +F90C=f90 +FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1 +LFLAGS= +endif +ifeq ($(ext),Sun) +F90C=f90 +FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none +LFLAGS= +endif -#Intel , -openmp toggles mutli-processor: -#note version 10.0 gives wrong result for lensed when compiled with -openmp [fixed in 10.1] -F90C = ifort -FFLAGS = -openmp -fast -W0 -WB -fpp2 -vec_report0 -ifneq ($(FISHER),) -FFLAGS += -mkl +ifeq ($(ext),AIX) +F90C=xlf90_r +FFLAGS= -O4 -q64 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90 +LFLAGS= endif -#Gfortran compiler: -#The options here work in v4.5, delete from RHS in earlier versions (15% slower) -#if pre v4.3 add -D__GFORTRAN__ -#With v4.6+ try -Ofast -march=native -fopenmp -#On my machine v4.5 is about 20% slower than ifort -#F90C = gfortran -#FFLAGS = -O3 -fopenmp -ffast-math -march=native -funroll-loops - - -#Old Intel ifc, add -openmp for multi-processor (some have bugs): -#F90C = ifc -#FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2 -#some systems can can also add e.g. -tpp7 -xW - -#G95 compiler -#F90C = g95 -#FFLAGS = -O2 - -#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems. -#F90C = f90 -#FFLAGS = -Ofast -mp - -#Digital/Compaq fortran, -omp toggles multi-processor -#F90C = f90 -#FFLAGS = -omp -O4 -arch host -math_library fast -tune host -fpe1 - -#Absoft ProFortran, single processor: -#F90C = f95 -#FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" - -#NAGF95, single processor: -#F90C = f95 -#FFLAGS = -DNAGF95 -O3 - -#PGF90 -#F90C = pgf90 -#FFLAGS = -O2 -DESCAPEBACKSLASH -Mpreprocess - -#Sun V880 -#F90C = mpf90 -#FFLAGS = -O4 -openmp -ftrap=%none -dalign -DMPI - -#Sun parallel enterprise: -#F90C = f95 -#FFLAGS = -O2 -xarch=native64 -openmp -ftrap=%none -#try removing -openmp if get bus errors. -03, -04 etc are dodgy. - -#IBM XL Fortran, multi-processor (run gmake) -#F90C = xlf90_r -#FFLAGS = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3 - -#Settings for building camb_fits -#Location of FITSIO and name of library -FITSDIR = /home/cpac/cpac-tools/lib -FITSLIB = cfitsio + +#Files containing evolution equations initial power spectrum module +EQUATIONS = equations +POWERSPECTRUM = power_inf +REIONIZATION = reionization +RECOMBINATION = recfast +BISPECTRUM = SeparableBispectrum +DENABLE_FISHER = + +#inf module +INFDIR = ../fieldinf +INFOBJ = infprec.o binfspline.o hyper2F1.o specialinf.o infinout.o inftools.o infbgmodel.o infsrmodel.o infbg.o infbgspline.o cosmopar.o inftorad.o infpert.o infpowspline.o + +#Module doing non-linear scaling +NONLINEAR = halofit + +#Driver program +DRIVER = inidriver.F90 +#DRIVER = sigma8.f90 +#DRIVER = tester.f90 + #Location of HEALPIX for building camb_fits -HEALPIXDIR = /home/cpac/cpac-tools/healpix +HEALPIXDIR = /usr -ifneq ($(FISHER),) +CAMBLIB = libcamb.a +INFLIB = libinf.a + +ifneq ($(DENABLE_FISHER),) FFLAGS += -DFISHER +LFLAGS += -llapack -lblas EXTCAMBFILES = Matrix_utils.o else EXTCAMBFILES = endif -include ./Makefile_main + +#Shouldn't need to change anything else... + +F90FLAGS = $(FFLAGS) -DPP5 +HEALPIXLD = -L$(HEALPIXDIR)/lib -lhealpix +FC = $(F90C) + +CAMBOBJ = constants.o utils.o subroutines.o inifile.o $(POWERSPECTRUM).o\ + $(RECOMBINATION).o $(REIONIZATION).o modules.o \ + bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o \ + cmbmain.o camb.o + + +default: camb.$(ext) + +all: camb.$(ext) $(CAMBLIB) $(INFLIB) + + +subroutines.o: constants.o utils.o +$(POWERSPECTRUM): subroutines.o inifile.o +$(RECOMBINATION).o: subroutines.o inifile.o +$(REIONIZATION).o: constants.o inifile.o +modules.o: $(REIONIZATION).o $(POWERSPECTRUM).o $(RECOMBINATION).o +bessels.o: modules.o +$(EQUATIONS): bessels.o +$(NONLINEAR).o: modules.o +lensing.o: bessels.o +$(BISPECTRUM).o: lensing.o modules.o +cmbmain.o: lensing.o $(NONLINEAR).o $(EQUATIONS).o $(BISPECTRUM).o +camb.o: cmbmain.o +Matrix_utils.o: utils.o + +camb.$(ext): $(INFOBJ) $(CAMBOBJ) $(DRIVER) + $(F90C) $(F90FLAGS) $(INFOBJ) $(CAMBOBJ) $(DRIVER) $(LFLAGS) -o $@ + +$(CAMBLIB): $(CAMBOBJ) + ar -r $@ $? + +$(INFLIB): $(INFOBJ) + ar -r $@ $? + +camb_fits.$(ext): writefits.f90 $(CAMBOBJ) $(DRIVER) + $(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) $(LFLAGS) -DWRITE_FITS -o $@ + +%.o: %.f90 + $(F90C) $(F90FLAGS) -c $< + +%.o: $(INFDIR)/%.f90 + $(F90C) $(F90FLAGS) -c $< + +%.o: $(INFDIR)/%.F90 + $(F90C) $(F90FLAGS) -c $< + +utils.o: + $(F90C) $(F90FLAGS) -c utils.F90 + +$(BISPECTRUM).o: + $(F90C) $(F90FLAGS) -c $(BISPECTRUM).F90 + +Matrix_utils.o: + $(F90C) $(F90FLAGS) -c Matrix_utils.F90 +clean: + -rm -f *.o *.a *.d core *.mod + diff --git a/camb/cmbmain.f90 b/camb/cmbmain.f90 index 90ce721..82a8aee 100644 --- a/camb/cmbmain.f90 +++ b/camb/cmbmain.f90 @@ -1991,13 +1991,17 @@ contains integer, intent(in) :: numks, pix real(dl) pows(numks), ks(numks) integer i - +!fields +!$omp parallel do & +!$omp default(shared) & +!$omp private(i) & +!$omp schedule(dynamic) do i = 1, numks !!change to vec... pows(i) = ScalarPower(ks(i) ,pix) - if (global_error_flag/=0) exit end do - +!$omp end parallel do +!end fields end subroutine GetInitPowerArrayVec @@ -2005,12 +2009,16 @@ contains integer, intent(in) :: numks, pix real(dl) pows(numks), ks(numks) integer i - +!fields +!$omp parallel do & +!$omp default(shared) & +!$omp private(i) & +!$omp schedule(dynamic) do i = 1, numks pows(i) = TensorPower(ks(i) ,pix) - if (global_error_flag/=0) exit end do - +!$omp end parallel do +!end fields end subroutine GetInitPowerArrayTens @@ -2026,7 +2034,11 @@ contains do pix=1,CP%InitPower%nn - +!field: calling ScalarPower is time consumming when it comes from inflation code. Better to parallelise +!$omp parallel do & +!$omp default(shared) & +!$omp private(q_ix) & +!$omp schedule(dynamic,1) do q_ix = 1, CTrans%q%npoints if (CP%flat) then @@ -2038,10 +2050,11 @@ contains end if pows(q_ix) = ScalarPower(ks(q_ix) ,pix) - if (global_error_flag/=0) return end do - +!$omp end parallel do + if (global_error_flag/=0) return +!end field !$OMP PARAllEl DO DEFAUlT(SHARED),SCHEDUlE(STATIC,4) & !$OMP & PRIVATE(j,q_ix,dlnk,apowers,ctnorm,dbletmp) diff --git a/camb/inidriver.F90 b/camb/inidriver.F90 index b355761..db0c05a 100644 --- a/camb/inidriver.F90 +++ b/camb/inidriver.F90 @@ -332,6 +332,10 @@ call CAMB_cleanup +!fields + call FreePowers(P%InitPower) +!end fields + end program driver diff --git a/camb/params.ini b/camb/params.ini index cf5b821..32158e8 100644 --- a/camb/params.ini +++ b/camb/params.ini @@ -1,12 +1,12 @@ #Parameters for CAMB #output_root is prefixed to output file names -output_root = test +output_root = fields #What to do get_scalar_cls = T get_vector_cls = F -get_tensor_cls = F +get_tensor_cls = T get_transfer = F #if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT} @@ -65,19 +65,6 @@ nu_mass_degeneracies = 0 #Fraction of total omega_nu h^2 accounted for by each eigenstate, eg. 0.5 0.5 nu_mass_fractions = 1 -#Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}. -initial_power_num = 1 -pivot_scalar = 0.05 -pivot_tensor = 0.05 -scalar_amp(1) = 2.1e-9 -scalar_spectral_index(1) = 0.96 -scalar_nrun(1) = 0 -tensor_spectral_index(1) = 0 -#ratio is that of the initial tens/scal power spectrum amplitudes -initial_ratio(1) = 1 -#note vector modes use the scalar settings above - - #Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5 reionization = T @@ -98,6 +85,63 @@ RECFAST_fudge_He = 0.86 RECFAST_Heswitch = 6 RECFAST_Hswitch = T +#inflationary model parameters (see infbgmodel.f90) +#param(1)=M; param(2)=p; param(3)=nu; param(4)=mu, param(5)=q +#largef: V = M^4 phi^p p > 0 +#smallf: V = M^4 [1 - (phi/mu)^p] p > 2 +#hybrid: V = M^4 [1 + (phi/mu)^p] +#runmas: V = M^4 {1 + nu [1/p + ln(phi/mu)]phi^p} +#kklmmt: V = M^4 [1 - q (phi/mu)^-p]^q p > 0; q=+-1 + +#number of read parameters +inf_param_number = 7 + +#should be consistent with the above defs +inf_model_name = largef +inf_param(1) = 2e-3 +inf_param(2) = 2 +inf_param(3) = 0 +inf_param(4) = 0 +inf_param(5) = 1 + + +#Hard prior that reject models with field values greater that +#field_bound. May be used for keeping models out of the +#self-reproducing regime, or constrain the branes to be inside the +#throat. In unit of kappa +inf_check_bound = F +inf_param(6) = 0 + + +#stops artificially inflation for this matter field value. Unit of +#kappa for large field and running mass; unit of mu for small fields +#and kklmmt; unit of a maximal field value that gives 120 efolds of +#hybrid inflation +inf_check_stop = F +inf_param(7) = 0 + +#initial fields value. Unit of kappa for large field and running; unit +#of mu for small field, hybrid and kklmmt. For inf_matter = 0, uses +#a slow-roll approximation to guess the relevant values +inf_conform(1) = 1 +inf_matter(1) = 0 + +#This is ln(aend/areh) - 1/4 ln(rhoreh kappa^4) + 1/2 ln(rhoend kappa^4) +inf_ln_reheat = 0 + + +#speed up A LOT the computations (only for featureless power spectra) +power_spectra_spline = T + +#if false then every k required are computed from their quantum birth +#during the inflationary era. Otherwise, only a few of them are +#computed and the others interpolated. The following parameters set +#the range and number (ln=log base e) +spline_lnkmpc_min = -14 +spline_lnkmpc_max = 0 +spline_lnkmpc_num = 12 + + #Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, # neutrino density iso =4, neutrino velocity iso = 5) initial_condition = 1 diff --git a/camb/readme.html b/camb/readme.html index f5ada26..cff517b 100644 --- a/camb/readme.html +++ b/camb/readme.html @@ -112,7 +112,7 @@ where all quantities are in the synchronous gauge and evaluated at the requested
halofit.f90
-Implements the NonLinear module, to calculate non linear scalings of the matter power spectrum as a function of redshift. Uses HALOFIT (astro-ph/0207664, code thanks to Robert Smith. Note this is only reliable at the several percent level for standard ΛCDM models with power law initial power spectra. This module can be replaced to use a different non-linear fitting method. +Implements the NonLinear module, to calculate non linear scalings of the matter power spectrum as a function of redshift. Uses HALOFIT (astro-ph/0207664, code thanks to Robert Smith, with tweaks from arXiv:1109.4416. Note this is only reliable at the several percent level for standard ΛCDM models with power law initial power spectra. This module can be replaced to use a different non-linear fitting method.
@@ -563,7 +563,9 @@ Implements the NonLinear module, to calculate non linear scalings of the matter
-Scalar errors should rarely exceed 0.3% for min(2500, L well into the damping tail) at default accuracy setting, and 0.1% for 500<L<2000 with high_accuracy_default=T. Matter power spectrum errors are usually dominated by interpolation in the acoustic oscillations, with about 0.2% accuracy with high_accuracy_default (but much better rms accuracy). +Scalar numerical errors should rarely exceed 0.3% for min(2500, L well into the damping tail) at default accuracy setting, and 0.1% for 500<L<2000 with high_accuracy_default=T. Matter power spectrum errors are usually dominated by interpolation in the acoustic oscillations, with about 0.2% accuracy with high_accuracy_default (but much better rms accuracy). For a detailed study of numerical accuracy as of January 2012 see arXiv:1201.3654. + + See also comparison with CMBFAST. Accuracy of course assumes the model is correct, and is dependent on RECFAST being the correct ionization history. Lensed C_l TT, TE and EE are accurate at the same level (to within the approximation that the lensing potential is linear, or the accuracy of the the HALOFIT non-linear model).Extreme models (e.g. scale > 4, h>1) may give errors of 5% or more. @@ -594,10 +596,15 @@ See also comparison with CMBFAST. Accuracy of course assu
REFERENCES
-Some notes and relevant Maple derivations are given here (see also the Appendix of astro-ph/0406096). The CAMB notes outline the equations and approximations used, and relation to standard synchronous-gauge and Newtonian-gauge variables. +Some notes and relevant Maple derivations are given here (see also the Appendix of astro-ph/0406096). The CAMB notes outline the equations and approximations used, and relation to standard synchronous-gauge and Newtonian-gauge variables; see also arXiv:1201.3654. There is a BibTex file of references (including CosmoMC).
+CMB power spectrum parameter degeneracies in the era of precision cosmology
+Cullan Howlett, Antony Lewis, Alex Hall, Anthony Challinor arXiv:1201.3654. +Efficient computation of CMB anisotropies in closed FRW Models
Antony Lewis, Anthony Challinor and Anthony Lasenby astro-ph/9911177 Ap. J. 538:473-476, 2000. @@ -647,7 +654,8 @@ Antony Lewis, astro-ph/0403583 HALOFITStable clustering, the halo model and nonlinear cosmological power spectra
-Smith, R. E. and others, astro-ph/0207664 +Smith, R. E. and others, astro-ph/0207664. +RECOMBINATION
@@ -696,13 +704,25 @@ The Cosmic Linear Anisotropy Solving System (CLASS) II: Blas, Diego and Lesgourgues, Julien and Tram, Thomas. arXiv:1104.2933
+Massive Neutrinos ++CMB power spectrum parameter degeneracies in the era of precision cosmology +
+Cullan Howlett, Antony Lewis, Alex Hall, Anthony Challinor. +arXiv:1201.3654 ++
Evolution of cosmological dark matter perturbations
+Antony Lewis and Anthony Challinor astro-ph/0203507 +Phys. Rev. D66, 023531 (2002) + +Synchronous gauge theory and non-flat models
Complete treatment of CMB anisotropies in a FRW universe
Wayne Hu, Uros Seljak and Matias Zaldarriaga. Phys. Rev. D57:6, 3290-3301, 1998. astro-ph/9709066. -
WKB approx to hyperspherical Bessel functions @@ -719,10 +739,16 @@ Blas, Diego and Lesgourgues, Julien and Tram, Thomas. astro-ph/9603033 Ap.J. 469:2 437-444, 1996
+ +Integral solution for the microwave background + anisotropies in nonflat universes
+ Matias Zaldarriaga, Uros Seljak, Edmund Bertschinger. + ApJ. 494:491-501, 1998. astro-ph/9704265. + +CMBFAST for spatially closed universes
Uros Seljak and Matias Zaldariaga, astro-ph/9911219 --See also the references on the CMBFAST home page. +