Need more accuracy in matter power spectrum
-
- Posts: 74
- Joined: May 18 2019
- Affiliation: IRAP
- Contact:
Need more accuracy in matter power spectrum
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
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
-
- Posts: 1941
- Joined: September 23 2004
- Affiliation: University of Sussex
- Contact:
Re: Need more accuracy in power matter spectrum
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).
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).
-
- Posts: 74
- Joined: May 18 2019
- Affiliation: IRAP
- Contact:
Re: Need more accuracy in power matter spectrum
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 ? :
If yes, how could I modify it ?
Any help is welcome, regards
Fabien Dournac
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
Any help is welcome, regards
Fabien Dournac
-
- Posts: 1941
- Joined: September 23 2004
- Affiliation: University of Sussex
- Contact:
Re: Need more accuracy in power matter spectrum
No.See Transfer_SaveMatterPower and trace back types of referenced objects from there.
-
- Posts: 74
- Joined: May 18 2019
- Affiliation: IRAP
- Contact:
Re: Need more accuracy in power matter spectrum
I took a look at function and tried the following modifications to get double precision or higher precision than original :
1 try )
2 try)
But in both cases, I get this compilation error :
Any idea about this issue and how to fix it ? thanks
Code: Select all
Transfer_SaveMatterPower
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
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
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
-
- Posts: 74
- Joined: May 18 2019
- Affiliation: IRAP
- Contact:
Re: Need more accuracy in power matter spectrum
Compilation problem has been fixed. Now, I am facing another issue :
I have modified the source results.f90 in the following way :
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 :
Once compilation passed, at the execution, I have this type of output data in files generated :
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 :
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
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
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
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
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
Thanks, Fabien Dournac
-
- Posts: 74
- Joined: May 18 2019
- Affiliation: IRAP
- Contact:
Re: Need more accuracy in power matter spectrum
ok, my Fortran memories are old, 'E15.10,1x' did the trick !