From 1754b9eef5bb819d138ffb04ea906dff0771c53b Mon Sep 17 00:00:00 2001 From: Antony Lewis Date: Mon, 10 Jul 2023 08:46:03 +0100 Subject: [PATCH] fix nasty bug getting Cl from time sources --- fortran/camb.f90 | 2 +- fortran/cmbmain.f90 | 41 ++++++++++++++++++++++++----------------- 2 files changed, 25 insertions(+), 18 deletions(-) diff --git a/fortran/camb.f90 b/fortran/camb.f90 index 97d08a2f..aefb0b1b 100644 --- a/fortran/camb.f90 +++ b/fortran/camb.f90 @@ -32,7 +32,7 @@ subroutine CAMB_TransfersToPowers(CData) Cdata%OnlyTransfer = .true. !prevent ClTransferToCl Cdata%CP%WantTensors = .false. CData%CP%WantVectors = .false. - call TimeSourcesToCl + call TimeSourcesToCl(CData%ClData%CTransScal) Cdata%CP%WantTensors = want_tensors CData%CP%WantVectors = want_vectors Cdata%OnlyTransfer = .false. diff --git a/fortran/cmbmain.f90 b/fortran/cmbmain.f90 index c209bc82..53c6c73e 100644 --- a/fortran/cmbmain.f90 +++ b/fortran/cmbmain.f90 @@ -105,7 +105,6 @@ module CAMBmain integer :: max_bessels_l_index = 1000000 real(dl) :: max_bessels_etak = 1000000*2 - Type(ClTransferData), pointer :: ThisCT => null() Type(TTimeSources) , pointer :: ThisSources => null() Type(TTimeSources), allocatable, target :: TempSources @@ -126,6 +125,7 @@ subroutine cmbmain type(EvolutionVars) EV Type(TTimer) :: Timer real(dl) starttime + Type(ClTransferData), pointer :: ThisCT WantLateTime = CP%DoLensing .or. State%num_redshiftwindows > 0 .or. CP%CustomSources%num_custom_sources>0 @@ -216,7 +216,7 @@ subroutine cmbmain ! integrating the sources over time and over k. if (CP%WantCls .and. (.not. CP%WantScalars .or. .not. State%HasScalarTimeSources)) then - call TimeSourcesToCl + call TimeSourcesToCl(ThisCT) if (CP%WantScalars) then deallocate(State%ScalarTimeSources) @@ -231,9 +231,10 @@ subroutine cmbmain end subroutine cmbmain - subroutine TimeSourcesToCl - Type(TTimer) :: Timer + subroutine TimeSourcesToCl(ThisCT) + Type(ClTransferData) :: ThisCT integer q_ix + Type(TTimer) :: Timer if (CP%WantScalars) ThisSources => State%ScalarTimeSources @@ -255,21 +256,21 @@ subroutine TimeSourcesToCl max_bessels_l_index = ThisCT%ls%nl max_bessels_etak = maximum_qeta - if (CP%WantScalars) call GetLimberTransfers + if (CP%WantScalars) call GetLimberTransfers(ThisCT) ThisCT%max_index_nonlimber = max_bessels_l_index if (State%flat) call InitSpherBessels(ThisCT%ls, CP, max_bessels_l_index,max_bessels_etak ) !This is only slow if not called before with same (or higher) Max_l, Max_eta_k !Preferably stick to Max_l being a multiple of 50 - call SetkValuesForInt + call SetkValuesForInt(ThisCT) if (DebugMsgs .and. Feedbacklevel > 0) call WriteFormat('Set %d integration k values',ThisCT%q%npoints) !Begin k-loop and integrate Sources*Bessels over time !$OMP PARALLEL DO DEFAULT(SHARED), SCHEDULE(STATIC,4) do q_ix=1,ThisCT%q%npoints - call SourceToTransfers(q_ix) + call SourceToTransfers(ThisCT, q_ix) end do !q loop !$OMP END PARALLEL DO @@ -398,7 +399,8 @@ subroutine CalcLimberScalCls(CTrans) end subroutine CalcLimberScalCls - subroutine GetLimberTransfers + subroutine GetLimberTransfers(ThisCT) + Type(ClTransferData), target :: ThisCT integer ell, ell_needed integer i, s_ix, s_ix_lens type(TRedWin), pointer :: W @@ -517,7 +519,8 @@ subroutine GetLimberTransfers end subroutine GetLimberTransfers - subroutine SourceToTransfers(q_ix) + subroutine SourceToTransfers(ThisCT, q_ix) + type(ClTransferData), target :: ThisCT integer q_ix type(IntegrationVars) :: IV @@ -532,7 +535,7 @@ subroutine SourceToTransfers(q_ix) call InterpolateSources(IV) - call DoSourceIntegration(IV) + call DoSourceIntegration(IV, ThisCT) if (.not.State%flat) deallocate(IV%ddSource_q) deallocate(IV%Source_q) @@ -1231,7 +1234,8 @@ subroutine InitSourceInterpolation end subroutine InitSourceInterpolation - subroutine SetkValuesForInt + subroutine SetkValuesForInt(ThisCT) + Type(ClTransferData) :: ThisCT integer no real(dl) dk,dk0,dlnk1, dk2, max_k_dk, k_max_log, k_max_0 integer lognum @@ -1392,10 +1396,11 @@ subroutine IntegrationVars_Init(IV) end subroutine IntegrationVars_Init - subroutine DoSourceIntegration(IV) !for particular wave number q + subroutine DoSourceIntegration(IV, ThisCT) !for particular wave number q + type(IntegrationVars) IV + Type(ClTransferData) :: ThisCT integer j,ll,llmax real(dl) nu - type(IntegrationVars) IV real(dl) :: sixpibynu nu=IV%q*State%curvature_radius @@ -1418,12 +1423,12 @@ subroutine DoSourceIntegration(IV) !for particular wave number q end if if (State%flat) then - call DoFlatIntegration(IV,llmax) + call DoFlatIntegration(IV,ThisCT, llmax) else do j=1,ThisCT%ls%nl ll=ThisCT%ls%l(j) if (ll>llmax) exit - call IntegrateSourcesBessels(IV,j,ll,nu) + call IntegrateSourcesBessels(IV,ThisCT,j,ll,nu) end do !j loop end if @@ -1449,9 +1454,10 @@ function UseLimber(l) end function UseLimber !flat source integration - subroutine DoFlatIntegration(IV, llmax) + subroutine DoFlatIntegration(IV, ThisCT, llmax) implicit none type(IntegrationVars) IV + Type(ClTransferData) :: ThisCT integer llmax integer j logical DoInt @@ -1641,9 +1647,10 @@ end subroutine DoFlatIntegration !non-flat source integration - subroutine IntegrateSourcesBessels(IV,j,l,nu) + subroutine IntegrateSourcesBessels(IV,ThisCT,j,l,nu) use SpherBessels type(IntegrationVars) IV + Type(ClTransferData) :: ThisCT logical DoInt integer l,j, nstart,nDissipative,ntop,nbot,nrange,nnow real(dl) nu,ChiDissipative,ChiStart,tDissipative,y1,y2,y1dis,y2dis