Need more accuracy in matter power spectrum

Use of Healpix, camb, CLASS, cosmomc, compilers, etc.
Post Reply
Dournac Fabien
Posts: 8
Joined: May 18 2019
Affiliation: IRAP
Contact:

Need more accuracy in matter power spectrum

Post by Dournac Fabien » May 18 2019

Hello,

I would like to get more digits with float or doubles values about the power matter spectrum I am generating. I am using CAMB-1.0.4 and the Fortran version.

someone suggested me to look for in file source results.f90 but I didn't find how to increase this accuracy.

Basiclally, here is what I want :

k(h/Mpc) Pk/s8^2(Mpc/h)^3
5.3594794736e-06 1.8529569626e+01
5.6332442000e-06 1.9437295914e+01
5.9209928622e-06 2.0389484405e+01
6.2234403231e-06 2.1388326645e+01
6.5413364609e-06 2.2436098099e+01
6.8754711720e-06 2.3535198212e+01
7.2266739153e-06 2.4688137054e+01
7.5958159869e-06 2.5897554398e+01
7.9838137026e-06 2.7166225433e+01
8.3916311269e-06 2.8497039795e+01
8.8202796178e-06 2.9893053055e+01
9.2708232842e-06 3.1357446670e+01
9.7443817140e-06 3.2893573761e+01

For the moment, I can't have as many digits (different from zero) and my current result is :

k(h/Mpc) Pk/s8^2(Mpc/h)^3
5.2781500000e-06 1.9477400000e+01
5.5479700000e-06 2.0432300000e+01
5.8315700000e-06 2.1434000000e+01
6.1296700000e-06 2.2484700000e+01
6.4430100000e-06 2.3587000000e+01
6.7723700000e-06 2.4743400000e+01
7.1185600000e-06 2.5956400000e+01
7.4824500000e-06 2.7228900000e+01
7.8649500000e-06 2.8563800000e+01
8.2669900000e-06 2.9964100000e+01

The values are slightly different but the problem is not about his difference, this is about the multiple zeros digits that I whish to replace real digits number.

If someone could give me any tracks or suggestions to solve this issue ...

Regards, Fabien Dournac

ps : I have also tried the "high_accuracy_default=T" option and "accuracy_boost = 3" but this doesn't seem to get more precision

Antony Lewis
Posts: 1522
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: Need more accuracy in power matter spectrum

Post by Antony Lewis » May 20 2019

The easiest way to avoid any rounding in text files is to use python version (or call CAMB programmatically).

Otherwise you can change the Fortran output code in results.f90 (however the matter transfer function is only calculated to single precision - though that should be far more precision than the numerical error).

Dournac Fabien
Posts: 8
Joined: May 18 2019
Affiliation: IRAP
Contact:

Re: Need more accuracy in power matter spectrum

Post by Dournac Fabien » June 03 2019

Thanks Antony,

unfortunately, the source results.f90 is relatively long and I don't know where to modifiy the "single" to "double precision" for matter power spectrum values : have I got to change this part of the file at line 255 ? :

Code: Select all

    interface
    FUNCTION state_function(obj, a)
    use precision
    import
    class(CAMBdata) :: obj
    real(dl), intent(in) :: a 
    real(dl) :: state_function
    END FUNCTION  state_function
    end interface
If yes, how could I modify it ?

Any help is welcome, regards

Fabien Dournac

Antony Lewis
Posts: 1522
Joined: September 23 2004
Affiliation: University of Sussex
Contact:

Re: Need more accuracy in power matter spectrum

Post by Antony Lewis » June 04 2019

No.See Transfer_SaveMatterPower and trace back types of referenced objects from there.

Dournac Fabien
Posts: 8
Joined: May 18 2019
Affiliation: IRAP
Contact:

Re: Need more accuracy in power matter spectrum

Post by Dournac Fabien » June 06 2019

I took a look at

Code: Select all

Transfer_SaveMatterPower
function and tried the following modifications to get double precision or higher precision than original :

1 try )

Code: Select all

subroutine Transfer_SaveMatterPower(MTrans, State,FileNames, all21cm)
use constants

!Export files of total  matter power spectra in h^{-1} Mpc units, against k/h.
Type(MatterTransferData), intent(in) :: MTrans
Type(CAMBdata) :: State
character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
character(LEN=name_tag_len) :: columns(3)
integer itf, i, unit
integer points
! Added : way of declaring double precision
integer, parameter :: wp = selected_real_kind(15,307)
real(wp), dimension(:,:), allocatable :: outpower
real(wp minkh,dlnkh
Type(MatterPowerData) :: PK_data
integer ncol
logical, intent(in), optional :: all21cm
logical all21
!JD 08/13 Changes in here to PK arrays and variables
integer itf_PK
2 try)

Code: Select all

subroutine Transfer_SaveMatterPower(MTrans, State,FileNames, all21cm)
use constants

!Export files of total  matter power spectra in h^{-1} Mpc units, against k/h.
Type(MatterTransferData), intent(in) :: MTrans
Type(CAMBdata) :: State
character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
character(LEN=name_tag_len) :: columns(3)
integer itf, i, unit
integer points
! Added : way of declaring double precision
double precision, dimension(:,:), allocatable :: outpower
double precision wp minkh,dlnkh
Type(MatterPowerData) :: PK_data
integer ncol
logical, intent(in), optional :: all21cm
logical all21
!JD 08/13 Changes in here to PK arrays and variables
integer itf_PK
But in both cases, I get this compilation error :

Code: Select all

../results.f90:3702:102:

                 call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf,  minkh,dlnkh, points)
                                                                                                      1
Error: Type mismatch in argument ‘outpower’ at (1); passed REAL(8) to REAL(4)
compilation terminated due to -fmax-errors=1.
make[1]: *** [results.o] Error 1
make: *** [camb] Error 2
Any idea about this issue and how to fix it ? thanks

Dournac Fabien
Posts: 8
Joined: May 18 2019
Affiliation: IRAP
Contact:

Re: Need more accuracy in power matter spectrum

Post by Dournac Fabien » June 06 2019

Compilation problem has been fixed. Now, I am facing another issue :

I have modified the source results.f90 in the following way :

Code: Select all

    module results
    use constants, only : const_pi, const_twopi
    use MiscUtils
    use RangeUtils
    use StringUtils
    use MathUtils
    use config
    use model
    implicit none 
    public

    integer, parameter :: wp = selected_real_kind(15,307)
    real(wp), dimension(:,:), allocatable :: outpower
    Type TBackgroundOutputs
        real(wp), allocatable :: H(:), DA(:), rs_by_D_v(:)
    end Type TBackgroundOutputs
and I replace into it all occurences of real(dp) by real(wp) to get double precision.

I have also replaced the original format E15.6 to E15.10 into Transfer_SaveMatterPower routine, like this :

Code: Select all

    subroutine Transfer_SaveMatterPower(MTrans, State,FileNames, all21cm)
    use constants
    !Export files of total  matter power spectra in h^{-1} Mpc units, against k/h.
    Type(MatterTransferData), intent(in) :: MTrans
    Type(CAMBdata) :: State
    character(LEN=Ini_max_string_len), intent(IN) :: FileNames(*)
    character(LEN=name_tag_len) :: columns(3)
    integer itf, i, unit
    integer points
    real, dimension(:,:), allocatable :: outpower
    real minkh,dlnkh
    Type(MatterPowerData) :: PK_data
    integer ncol
    logical, intent(in), optional :: all21cm
    logical all21
    !JD 08/13 Changes in here to PK arrays and variables
    integer itf_PK

    all21 = DefaultFalse(all21cm)
    if (all21) then
        ncol = 3
    else
        ncol = 1
    end if

    do itf=1, State%CP%Transfer%PK_num_redshifts
        if (FileNames(itf) /= '') then
            if (.not. transfer_interp_matterpower ) then
                itf_PK = State%PK_redshifts_index(itf)

                points = MTrans%num_q_trans
                allocate(outpower(points,ncol))

                !Sources
                if (all21) then
                    call Transfer_Get21cmPowerData(MTrans, State, PK_data, itf_PK)
                else
                    call Transfer_GetMatterPowerData(State, MTrans, PK_data, itf_PK)
                    !JD 08/13 for nonlinear lensing of CMB + LSS compatibility
                    !Changed (CP%NonLinear/=NonLinear_None) to CP%NonLinear/=NonLinear_none .and. CP%NonLinear/=NonLinear_Lens)
                    if(State%CP%NonLinear/=NonLinear_none .and. State%CP%NonLinear/=NonLinear_Lens) then
                        call State%CP%NonLinearModel%GetNonLinRatios(State, PK_data)
                        PK_data%matpower = PK_data%matpower +  2*log(PK_data%nonlin_ratio)
                        call MatterPowerdata_getsplines(PK_data)
                    end if
                end if

                outpower(:,1) = exp(PK_data%matpower(:,1))
                !Sources
                if (all21) then
                    outpower(:,3) = exp(PK_data%vvpower(:,1))
                    outpower(:,2) = exp(PK_data%vdpower(:,1))

                    outpower(:,1) = outpower(:,1)/1d10*const_pi*const_twopi/MTrans%TransferData(Transfer_kh,:,1)**3
                    outpower(:,2) = outpower(:,2)/1d10*const_pi*const_twopi/MTrans%TransferData(Transfer_kh,:,1)**3
                    outpower(:,3) = outpower(:,3)/1d10*const_pi*const_twopi/MTrans%TransferData(Transfer_kh,:,1)**3
                end if

                call MatterPowerdata_Free(PK_Data)
                columns = ['P   ', 'P_vd','P_vv']
                unit = open_file_header(FileNames(itf), 'k/h', columns(:ncol), 15)
                do i=1,points
                    write (unit, '(*(E15.10))') MTrans%TransferData(Transfer_kh,i,1),outpower(i,:)
                end do
                close(unit)
            else
                if (all21) stop 'Transfer_SaveMatterPower: if output all assume not interpolated'
                minkh = 1e-4
                dlnkh = 0.02
                points = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/dlnkh+1
                !             dlnkh = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/(points-0.999)
                allocate(outpower(points,1))
                call Transfer_GetMatterPowerS(State, MTrans, outpower(1,1), itf,  minkh,dlnkh, points)

                columns(1) = 'P'
                unit = open_file_header(FileNames(itf), 'k/h', columns(:1), 15)

                do i=1,points
                    write (unit, '(*(E15.10))') minkh*exp((i-1)*dlnkh),outpower(i,1)
                end do
                close(unit)
            end if

            deallocate(outpower)
        end if
    end do

    end subroutine Transfer_SaveMatterPower
Once compilation passed, at the execution, I have this type of output data in files generated :

Code: Select all

.5278154731E-05.1454488945E+02
.5547966339E-05.1525794315E+02
.5831570434E-05.1600595665E+02
.6129671874E-05.1679063797E+02
.6443011898E-05.1761379051E+02
.6772369488E-05.1847729492E+02
.7118563190E-05.1938313484E+02
.7482453839E-05.2033338356E+02
.7864946383E-05.2133021545E+02
.8266991244E-05.2237591171E+02
.8689587958E-05.2347287369E+02
.9133786989E-05.2462361717E+02
.9600692465E-05.2583077240E+02
.1009146672E-04.2709711075E+02
.1060732757E-04.2842552185E+02
.1114955921E-04.2981906509E+02
.1171950862E-04.3128091812E+02
.1231859278E-04.3281444168E+02
.1294830145E-04.3442313766E+02
.1361019986E-04.3611070633E+02
.1430593420E-04.3788099289E+02
.1503723252E-04.3973807144E+02
.1580591379E-04.4168618774E+02
.1661389069E-04.4372981262E+02
.1746316775E-04.4587360382E+02
It seems that the 2 columns of matter power spectrum have not enough space to be separated by an empty space character.

1) Why the width of the 2 columns is fixed (and where in the source ?) ?

If I set in results.f90 the original format (E15.6), everything is ok since I get into test_matterpower*.dat files :

Code: Select all

   0.527815E-05   0.797896E+01
   0.554797E-05   0.837012E+01
   0.583157E-05   0.878046E+01
   0.612967E-05   0.921092E+01
   0.644301E-05   0.966248E+01
   0.677237E-05   0.101362E+02
   0.711856E-05   0.106331E+02
   0.748245E-05   0.111544E+02
   0.786495E-05   0.117012E+02
   0.826699E-05   0.122749E+02
   0.868959E-05   0.128766E+02
   0.913379E-05   0.135079E+02
   0.960069E-05   0.141701E+02
   0.100915E-04   0.148648E+02
   0.106073E-04   0.155935E+02
   0.111496E-04   0.163580E+02
   0.117195E-04   0.171599E+02
   0.123186E-04   0.180012E+02
   0.129483E-04   0.188837E+02
2) So, as conclusion, how to keep the empty space characters between the 2 columns into test_matterpower*.dat output files, with the format E15.10 (original format is E15.6) ?

Thanks, Fabien Dournac

Dournac Fabien
Posts: 8
Joined: May 18 2019
Affiliation: IRAP
Contact:

Re: Need more accuracy in power matter spectrum

Post by Dournac Fabien » June 06 2019

ok, my Fortran memories are old, 'E15.10,1x' did the trick !

Post Reply