diff --git a/camb/_compilers.py b/camb/_compilers.py index 61a45cdd..1db45f4b 100644 --- a/camb/_compilers.py +++ b/camb/_compilers.py @@ -4,8 +4,6 @@ import struct import subprocess -from camb.baseconfig import CAMBError - is_windows = platform.system() == "Windows" is_32_bit = struct.calcsize("P") == 4 @@ -86,7 +84,7 @@ def check_gfortran(version=gfortran_min, msg=False, retry=False): version_str = str(subprocess.check_output("gfortran -dumpmachine", shell=True, env=compiler_environ)) ok = gfortran_bits in version_str if not ok and msg: - raise CAMBError( + raise RuntimeError( f"You need ifort or gfortran {version} or higher to compile (found: {gfortran_version}).\n" "See https://camb.readthedocs.io/en/latest/fortran_compilers.htm\n" "or install from pypi using pip ('pip install camb') to just use pre-built binary wheels." diff --git a/camb/baseconfig.py b/camb/baseconfig.py index 7e0a7e1b..a3bbb2a4 100644 --- a/camb/baseconfig.py +++ b/camb/baseconfig.py @@ -200,7 +200,7 @@ def AllocatableObject(cls=None): class _AllocatableArray(FortranAllocatable): # member corresponding to allocatable :: d(:) member in fortran - _fields_ = ("allocatable", ctypes.c_void_p * (_f_allocatable_array_size // ctypes.sizeof(ctypes.c_void_p))) + _fields_ = (("allocatable", ctypes.c_void_p * (_f_allocatable_array_size // ctypes.sizeof(ctypes.c_void_p))),) def get_allocatable(self): size = self._get_allocatable_1D_array(byref(self), byref(_reuse_pointer)) diff --git a/camb/initialpower.py b/camb/initialpower.py index 51746f84..7ff27cdc 100644 --- a/camb/initialpower.py +++ b/camb/initialpower.py @@ -190,8 +190,11 @@ def set_params( # set from inflationary consistency if ntrun: raise CAMBError("ntrun set but using inflation consistency (nt=None)") - if self.tensor_parameterization != tensor_param_rpivot: - raise CAMBError("tensor parameterization not tensor_param_rpivot with inflation consistency") + if self.tensor_parameterization != "tensor_param_rpivot": + raise CAMBError( + "tensor parameterization not tensor_param_rpivot with inflation consistency: " + f"{self.tensor_parameterization}" + ) self.nt = -r / 8.0 * (2.0 - ns - r / 8.0) self.ntrun = r / 8.0 * (r / 8.0 + ns - 1) else: diff --git a/camb/nonlinear.py b/camb/nonlinear.py index 0b1a6b83..695edbce 100644 --- a/camb/nonlinear.py +++ b/camb/nonlinear.py @@ -1,6 +1,9 @@ -from ctypes import c_double, c_int +from ctypes import POINTER, byref, c_double, c_int -from .baseconfig import F2003Class, fortran_class +import numpy as np +from numpy.ctypeslib import ndpointer + +from .baseconfig import F2003Class, fortran_class, numpy_1d class NonLinearModel(F2003Class): @@ -107,3 +110,57 @@ class SecondOrderPK(NonLinearModel): def set_params(self): pass + + +@fortran_class +class ExternalNonLinearRatio(NonLinearModel): + """ + Non-linear model that applies a user-supplied ratio sqrt(P_NL/P_L) + from an external source (e.g. CCL, axionHMcode). + + Use :meth:`set_ratio` to provide the ratio grid, then assign this + as ``params.NonLinearModel`` before calling :func:`camb.get_results`. + Can be called after only time transfers computed to then get + lensed C_l with a consistent non-linear model. If linear P_k sampling + points in k_h are used, no interpolation from provided grid is required. + """ + + _fortran_class_module_ = "ExternalNonLinearRatio" + _fortran_class_name_ = "TExternalNonLinearRatio" + + _methods_ = ( + ( + "SetRatio", + [ + POINTER(c_int), + POINTER(c_int), + numpy_1d, + numpy_1d, + ndpointer(c_double, flags="F_CONTIGUOUS", ndim=2), + ], + ), + ("ClearRatio", []), + ) + + def set_ratio(self, k_h, z, ratio): + """ + Set the non-linear ratio grid sqrt(P_NL/P_L). + + :param k_h: 1D array of k values in h/Mpc units (ascending) + :param z: 1D array of redshift values (ascending) + :param ratio: 2D array of sqrt(P_NL/P_L), shape (len(z), len(k_h)), + matching the convention of CAMB's get_matter_power_spectrum + """ + k_h = np.ascontiguousarray(k_h, dtype=np.float64) + z = np.ascontiguousarray(z, dtype=np.float64) + if ratio.shape != (len(z), len(k_h)): + raise ValueError(f"ratio shape {ratio.shape} must be (len(z), len(k_h)) = ({len(z)}, {len(k_h)})") + # Fortran expects (nk, nz) column-major; C-order (nz, nk) has the same memory layout + ratio_f = np.asfortranarray(ratio.T, dtype=np.float64) + self.f_SetRatio(byref(c_int(len(k_h))), byref(c_int(len(z))), k_h, z, ratio_f) + + def clear_ratio(self): + """ + Clear the stored ratio, releasing interpolation data. + """ + self.f_ClearRatio() diff --git a/fortran/ExternalNonLinearRatio.f90 b/fortran/ExternalNonLinearRatio.f90 new file mode 100644 index 00000000..2efff8c9 --- /dev/null +++ b/fortran/ExternalNonLinearRatio.f90 @@ -0,0 +1,100 @@ + + module ExternalNonLinearRatio + ! Non-linear model that applies a user-supplied ratio sqrt(P_NL/P_L) + ! from an external source (e.g. CCL, axionHMcode). + ! The ratio is stored as a 2D interpolation grid over (k/h, z). + use results + use transfer + use Interpolation, only : TInterpGrid2D + implicit none + private + + type, extends(TNonLinearModel) :: TExternalNonLinearRatio + Type(TInterpGrid2D) :: Ratio + logical :: ratio_set = .false. + contains + procedure :: GetNonLinRatios => TExternalNonLinearRatio_GetNonLinRatios + procedure :: GetNonLinRatios_All => TExternalNonLinearRatio_GetNonLinRatios_All + procedure :: SetRatio => TExternalNonLinearRatio_SetRatio + procedure :: ClearRatio => TExternalNonLinearRatio_ClearRatio + procedure, nopass :: SelfPointer => TExternalNonLinearRatio_SelfPointer + end type TExternalNonLinearRatio + + public TExternalNonLinearRatio + contains + + subroutine TExternalNonLinearRatio_SelfPointer(cptr, P) + use iso_c_binding + Type(c_ptr) :: cptr + Type(TExternalNonLinearRatio), pointer :: PType + class(TPythonInterfacedClass), pointer :: P + + call c_f_pointer(cptr, PType) + P => PType + + end subroutine TExternalNonLinearRatio_SelfPointer + + subroutine TExternalNonLinearRatio_SetRatio(this, nk, nz, k_h_arr, z_arr, ratio_arr) + ! Set the non-linear ratio grid sqrt(P_NL/P_L) + class(TExternalNonLinearRatio) :: this + integer, intent(in) :: nk, nz + real(dl), intent(in) :: k_h_arr(nk) + real(dl), intent(in) :: z_arr(nz) + real(dl), intent(in) :: ratio_arr(nk, nz) + + call this%Ratio%Init(k_h_arr, z_arr, ratio_arr) + this%ratio_set = .true. + + end subroutine TExternalNonLinearRatio_SetRatio + + subroutine TExternalNonLinearRatio_ClearRatio(this) + ! Clear the stored ratio, releasing interpolation data + class(TExternalNonLinearRatio) :: this + + call this%Ratio%Clear() + this%ratio_set = .false. + + end subroutine TExternalNonLinearRatio_ClearRatio + + subroutine TExternalNonLinearRatio_GetNonLinRatios(this, State, CAMB_Pk) + ! Fill CAMB_Pk%nonlin_ratio from the stored interpolation grid. + ! Clamps k and z to the grid range to prevent unreliable polynomial extrapolation. + class(TExternalNonLinearRatio) :: this + class(TCAMBdata) :: State + type(MatterPowerData), target :: CAMB_Pk + integer :: ik, iz + real(dl) :: kh, z, kh_clamped, z_clamped + real(dl) :: k_min, k_max, z_min, z_max + + if (.not. this%ratio_set) & + error stop 'ExternalNonLinearRatio: ratio not set. Call SetRatio first.' + + ! Get bounds of the ratio grid to clamp values + k_min = this%Ratio%x(1) + k_max = this%Ratio%x(this%Ratio%nx) + z_min = this%Ratio%y(1) + z_max = this%Ratio%y(this%Ratio%ny) + + CAMB_Pk%nonlin_ratio = 1.0_dl + do iz = 1, CAMB_Pk%num_z + z = CAMB_Pk%redshifts(iz) + z_clamped = max(z_min, min(z_max, z)) + do ik = 1, CAMB_Pk%num_k + kh = exp(CAMB_Pk%log_kh(ik)) + kh_clamped = max(k_min, min(k_max, kh)) + CAMB_Pk%nonlin_ratio(ik, iz) = this%Ratio%Value(kh_clamped, z_clamped) + end do + end do + + end subroutine TExternalNonLinearRatio_GetNonLinRatios + + subroutine TExternalNonLinearRatio_GetNonLinRatios_All(this, State, CAMB_Pk) + class(TExternalNonLinearRatio) :: this + class(TCAMBdata) :: State + type(MatterPowerData), target :: CAMB_Pk + + error stop 'ExternalNonLinearRatio: GetNonLinRatios_All not supported (no velocity corrections)' + + end subroutine TExternalNonLinearRatio_GetNonLinRatios_All + + end module ExternalNonLinearRatio diff --git a/fortran/Makefile_main b/fortran/Makefile_main index 64d1f066..745c6477 100644 --- a/fortran/Makefile_main +++ b/fortran/Makefile_main @@ -5,7 +5,7 @@ POWERSPECTRUM_FILES ?= InitialPower REIONIZATION_FILES ?= reionization RECOMBINATION_FILES ?= recfast -NONLINEAR_FILES ?= halofit SecondOrderPK +NONLINEAR_FILES ?= halofit SecondOrderPK ExternalNonLinearRatio DARKENERGY_FILES ?= DarkEnergyFluid DarkEnergyPPF PowellMinimize DarkEnergyQuintessence BISPECTRUM ?= SeparableBispectrum