From b94db26db76498da3733304fb287b1b99e16de77 Mon Sep 17 00:00:00 2001 From: Jobilthomas Date: Fri, 21 Apr 2023 20:55:50 +0530 Subject: [PATCH 1/4] primordial_inflation_find_phi_pivot bug fix --- Makefile | 206 --------------- source/primordial.c | 602 ++++++++++++++++---------------------------- 2 files changed, 211 insertions(+), 597 deletions(-) delete mode 100644 Makefile diff --git a/Makefile b/Makefile deleted file mode 100644 index 9e4a8a779..000000000 --- a/Makefile +++ /dev/null @@ -1,206 +0,0 @@ -#Some Makefile for CLASS. -#Julien Lesgourgues, 28.11.2011 -#Nils Schöneberg, Matteo Lucca, 27.02.2019 - -MDIR := $(shell pwd) -WRKDIR = $(MDIR)/build - -.base: - if ! [ -e $(WRKDIR) ]; then mkdir $(WRKDIR) ; mkdir $(WRKDIR)/lib; fi; - touch build/.base - -vpath %.c source:tools:main:test -vpath %.o build -vpath .base build - -######################################################## -###### LINES TO ADAPT TO YOUR PLATEFORM ################ -######################################################## - -# your C compiler: -CC = gcc -#CC = icc -#CC = pgcc - -# your tool for creating static libraries: -AR = ar rv - -# Your python interpreter. -# In order to use Python 3, you can manually -# substitute python3 to python in the line below, or you can simply -# add a compilation option on the terminal command line: -# "PYTHON=python3 make all" (Thanks to Marius Millea for python3 compatibility) -PYTHON ?= python - -# your optimization flag -OPTFLAG = -O3 -#OPTFLAG = -Ofast -ffast-math #-march=native -#OPTFLAG = -fast - -# your openmp flag (comment for compiling without openmp) -OMPFLAG = -fopenmp -#OMPFLAG = -mp -mp=nonuma -mp=allcores -g -#OMPFLAG = -openmp - -# all other compilation flags -CCFLAG = -g -fPIC -LDFLAG = -g -fPIC - -# leave blank to compile without HyRec, or put path to HyRec directory -# (with no slash at the end: e.g. "external/RecfastCLASS") -HYREC = external/HyRec2020 -RECFAST = external/RecfastCLASS -HEATING = external/heating - -######################################################## -###### IN PRINCIPLE THE REST SHOULD BE LEFT UNCHANGED ## -######################################################## - -# pass current working directory to the code -CCFLAG += -D__CLASSDIR__='"$(MDIR)"' - -# where to find include files *.h -INCLUDES = -I../include -HEADERFILES = $(wildcard ./include/*.h) - -# automatically add external programs if needed. First, initialize to blank. -EXTERNAL = - -vpath %.c $(RECFAST) -#CCFLAG += -DRECFAST -INCLUDES += -I../$(RECFAST) -EXTERNAL += wrap_recfast.o -HEADERFILES += $(wildcard ./$(RECFAST)/*.h) - -vpath %.c $(HEATING) -#CCFLAG += -DHEATING -INCLUDES += -I../$(HEATING) -EXTERNAL += injection.o noninjection.o -HEADERFILES += $(wildcard ./$(HEATING)/*.h) - -# update flags for including HyRec -ifneq ($(HYREC),) -vpath %.c $(HYREC) -CCFLAG += -DHYREC -#LDFLAGS += -DHYREC -INCLUDES += -I../$(HYREC) -EXTERNAL += hyrectools.o helium.o hydrogen.o history.o wrap_hyrec.o energy_injection.o -HEADERFILES += $(wildcard ./$(HYREC)/*.h) -endif - -%.o: %.c .base $(HEADERFILES) - cd $(WRKDIR);$(CC) $(OPTFLAG) $(OMPFLAG) $(CCFLAG) $(INCLUDES) -c ../$< -o $*.o - -TOOLS = growTable.o dei_rkck.o sparse.o evolver_rkck.o evolver_ndf15.o arrays.o parser.o quadrature.o hyperspherical.o common.o trigonometric_integrals.o - -SOURCE = input.o background.o thermodynamics.o perturbations.o primordial.o fourier.o transfer.o harmonic.o lensing.o distortions.o - -INPUT = input.o - -PRECISION = precision.o - -BACKGROUND = background.o - -THERMO = thermodynamics.o - -PERTURBATIONS = perturbations.o - -TRANSFER = transfer.o - -PRIMORDIAL = primordial.o - -HARMONIC = harmonic.o - -FOURIER = fourier.o - -LENSING = lensing.o - -DISTORTIONS = distortions.o - -OUTPUT = output.o - -CLASS = class.o - -TEST_LOOPS = test_loops.o - -TEST_LOOPS_OMP = test_loops_omp.o - -TEST_HARMONIC = test_harmonic.o - -TEST_TRANSFER = test_transfer.o - -TEST_FOURIER = test_fourier.o - -TEST_PERTURBATIONS = test_perturbations.o - -TEST_THERMODYNAMICS = test_thermodynamics.o - -TEST_BACKGROUND = test_background.o - -TEST_HYPERSPHERICAL = test_hyperspherical.o - -C_TOOLS = $(addprefix tools/, $(addsuffix .c,$(basename $(TOOLS)))) -C_SOURCE = $(addprefix source/, $(addsuffix .c,$(basename $(SOURCE) $(OUTPUT)))) -C_TEST = $(addprefix test/, $(addsuffix .c,$(basename $(TEST_DEGENERACY) $(TEST_LOOPS) $(TEST_TRANSFER) $(TEST_FOURIER) $(TEST_PERTURBATIONS) $(TEST_THERMODYNAMICS)))) -C_MAIN = $(addprefix main/, $(addsuffix .c,$(basename $(CLASS)))) -C_ALL = $(C_MAIN) $(C_TOOLS) $(C_SOURCE) -H_ALL = $(addprefix include/, common.h svnversion.h $(addsuffix .h, $(basename $(notdir $(C_ALL))))) -PRE_ALL = cl_ref.pre clt_permille.pre -INI_ALL = explanatory.ini lcdm.ini -MISC_FILES = Makefile CPU psd_FD_single.dat myselection.dat myevolution.dat README bbn/sBBN.dat external_Pk/* cpp -PYTHON_FILES = python/classy.pyx python/setup.py python/cclassy.pxd python/test_class.py - -all: class libclass.a classy - -libclass.a: $(TOOLS) $(SOURCE) $(EXTERNAL) - $(AR) $@ $(addprefix build/, $(TOOLS) $(SOURCE) $(EXTERNAL)) - -class: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(CLASS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o class $(addprefix build/,$(notdir $^)) -lm - -test_loops: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(TEST_LOOPS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_loops_omp: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(TEST_LOOPS_OMP) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_harmonic: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_HARMONIC) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_transfer: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_TRANSFER) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_fourier: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_FOURIER) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_perturbations: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_PERTURBATIONS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_thermodynamics: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_THERMODYNAMICS) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_background: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_BACKGROUND) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm - -test_hyperspherical: $(TOOLS) $(TEST_HYPERSPHERICAL) - $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o test_hyperspherical $(addprefix build/,$(notdir $^)) -lm - - -tar: $(C_ALL) $(C_TEST) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES) - tar czvf class.tar.gz $(C_ALL) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES) - -classy: libclass.a python/classy.pyx python/cclassy.pxd -ifdef OMPFLAG - cp python/setup.py python/autosetup.py -else - grep -v "lgomp" python/setup.py > python/autosetup.py -endif - cd python; export CC=$(CC); $(PYTHON) autosetup.py install || $(PYTHON) autosetup.py install --user - rm python/autosetup.py - -clean: .base - rm -rf $(WRKDIR); - rm -f libclass.a - rm -f $(MDIR)/python/classy.c - rm -rf $(MDIR)/python/build - rm -f python/autosetup.py diff --git a/source/primordial.c b/source/primordial.c index 181ad29ca..0838656ce 100644 --- a/source/primordial.c +++ b/source/primordial.c @@ -2570,14 +2570,13 @@ int primordial_inflation_find_phi_pivot( ppm->error_message, ppm->error_message); - /** - case in which epsilon>1: hence we must find the value phi_stop < - phi_end where inflation ends up naturally */ + // assume that inflation ends up naturally - if (epsilon > 1.) { + /** - --> find latest value of the field such that epsilon = primordial_inflation_small_epsilon (default: 0.1) */ - // assume that inflation ends up naturally - /** - --> find latest value of the field such that epsilon = primordial_inflation_small_epsilon (default: 0.1) */ + if (epsilon > ppr->primordial_inflation_small_epsilon) { + /** - --> bracketing right-hand value is phi_end (but the potential will not be evaluated exactly there, only closeby */ phi_right = ppm->phi_end; @@ -2591,457 +2590,278 @@ int primordial_inflation_find_phi_pivot( ppm->error_message); } while (epsilon > ppr->primordial_inflation_small_epsilon); phi_left = ppm->phi_end-dphi; + } + else{ + /** - --> bracketing left-hand value is phi_end (but the potential will not be evaluated exactly there, only closeby */ + phi_left = ppm->phi_end; - /** - --> find value such that epsilon = primordial_inflation_small_epsilon by bisection */ + /** - --> bracketing right-hand value is found by iterating with logarithmic step until epsilon < primordial_inflation_small_epsilon */ + dphi = ppr->primordial_inflation_end_dphi; do { - phi_mid = 0.5*(phi_left+phi_right); - class_call(primordial_inflation_get_epsilon(ppm,phi_mid,&epsilon), + dphi *= ppr->primordial_inflation_end_logstep; + class_call(primordial_inflation_get_epsilon(ppm,ppm->phi_end+dphi,&epsilon), ppm->error_message, ppm->error_message); - if (epsilon < ppr->primordial_inflation_small_epsilon) phi_left=phi_mid; - else phi_right=phi_mid; - } while (fabs(epsilon-ppr->primordial_inflation_small_epsilon) > ppr->primordial_inflation_small_epsilon_tol); - - /** - --> value found and stored as phi_small_epsilon */ - phi_small_epsilon = phi_mid; - - /** - --> find inflationary attractor in phi_small_epsilon (should exist since epsilon<<1 there) */ - class_call(primordial_inflation_find_attractor(ppm, - ppr, - phi_small_epsilon, - ppr->primordial_inflation_attractor_precision_initial, - y, - dy, - &H_small_epsilon, - &dphidt_small_epsilon), - ppm->error_message, - ppm->error_message); - - /** - --> compute amount of inflation between this phi_small_epsilon and the end of inflation */ - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= phi_small_epsilon; - y[ppm->index_in_dphi]=y[ppm->index_in_a]*dphidt_small_epsilon; - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _end_inflation_, - 0., - _FALSE_, - forward, - conformal), - ppm->error_message, - ppm->error_message); - - // we have used here conformal time, so aH = dy[a]/y[a] - aH_ratio_after_small_epsilon = dy[ppm->index_in_a]/y[ppm->index_in_a]/H_small_epsilon; - a_ratio_after_small_epsilon = y[ppm->index_in_a]; - - switch (ppm->phi_pivot_method) { + } while (epsilon < ppr->primordial_inflation_small_epsilon); + phi_right = ppm->phi_end+dphi; + } + /** - --> find value such that epsilon = primordial_inflation_small_epsilon by bisection */ + do { + phi_mid = 0.5*(phi_left+phi_right); + class_call(primordial_inflation_get_epsilon(ppm,phi_mid,&epsilon), + ppm->error_message, + ppm->error_message); + if (epsilon < ppr->primordial_inflation_small_epsilon) phi_left=phi_mid; + else phi_right=phi_mid; + } while (fabs(epsilon-ppr->primordial_inflation_small_epsilon) > ppr->primordial_inflation_small_epsilon_tol); + + /** - --> value found and stored as phi_small_epsilon */ + phi_small_epsilon = phi_mid; + + /** - --> find inflationary attractor in phi_small_epsilon (should exist since epsilon<<1 there) */ + class_call(primordial_inflation_find_attractor(ppm, + ppr, + phi_small_epsilon, + ppr->primordial_inflation_attractor_precision_initial, + y, + dy, + &H_small_epsilon, + &dphidt_small_epsilon), + ppm->error_message, + ppm->error_message); - case ln_aH_ratio_auto: + /** - --> compute amount of inflation between this phi_small_epsilon and the end of inflation */ + y[ppm->index_in_a]=1.; + y[ppm->index_in_phi]= phi_small_epsilon; + y[ppm->index_in_dphi]=y[ppm->index_in_a]*dphidt_small_epsilon; - /* get the target value of ln_aH_ratio */ + class_call(primordial_inflation_evolve_background(ppm, + ppr, + y, + dy, + _end_inflation_, + 0., + _FALSE_, + forward, + conformal), + ppm->error_message, + ppm->error_message); - rho_end = 2./8./_PI_*pow(dy[ppm->index_in_a]/y[ppm->index_in_a],2); - rho_end = 8*_PI_/3.*rho_end/(_G_*_h_P_/pow(_c_,3))*pow(_Mpc_over_m_,2); - h = 0.7; - H0 = h * 1.e5 / _c_; - rho_c0 = pow(H0,2); + // we have used here conformal time, so aH = dy[a]/y[a] + aH_ratio_after_small_epsilon = dy[ppm->index_in_a]/y[ppm->index_in_a]/H_small_epsilon; + a_ratio_after_small_epsilon = y[ppm->index_in_a]; - sigma_B = 2. * pow(_PI_,5) * pow(_k_B_,4) / 15. / pow(_h_P_,3) / pow(_c_,2); - Omega_g0 = (4.*sigma_B/_c_*pow(2.726,4.)) / (3.*_c_*_c_*1.e10*h*h/_Mpc_over_m_/_Mpc_over_m_/8./_PI_/_G_); - Omega_r0 = 3.044*7./8.*pow(4./11.,4./3.)*Omega_g0; + switch (ppm->phi_pivot_method) { - target = log(H0/0.05*pow(Omega_r0,0.5)*pow(2./100.,1./12.)*pow(rho_end/rho_c0,0.25)); + case ln_aH_ratio_auto: - //fprintf(stderr,"auto: log(aH_end/aH_*)=%e\n",target); - break; + /* get the target value of ln_aH_ratio */ - case ln_aH_ratio: + rho_end = 2./8./_PI_*pow(dy[ppm->index_in_a]/y[ppm->index_in_a],2); + rho_end = 8*_PI_/3.*rho_end/(_G_*_h_P_/pow(_c_,3))*pow(_Mpc_over_m_,2); + h = 0.7; + H0 = h * 1.e5 / _c_; + rho_c0 = pow(H0,2); - target = ppm->phi_pivot_target; - //fprintf(stderr,"fixed: log(aH_end/aH_*)=%e\n",target); - break; + sigma_B = 2. * pow(_PI_,5) * pow(_k_B_,4) / 15. / pow(_h_P_,3) / pow(_c_,2); + Omega_g0 = (4.*sigma_B/_c_*pow(2.726,4.)) / (3.*_c_*_c_*1.e10*h*h/_Mpc_over_m_/_Mpc_over_m_/8./_PI_/_G_); + Omega_r0 = 3.044*7./8.*pow(4./11.,4./3.)*Omega_g0; - case N_star: + target = log(H0/0.05*pow(Omega_r0,0.5)*pow(2./100.,1./12.)*pow(rho_end/rho_c0,0.25)); - target = ppm->phi_pivot_target; - //fprintf(stderr,"fixed: log(a_end/a_*)=%e\n",target); - break; - } + //fprintf(stderr,"auto: log(aH_end/aH_*)=%e\n",target); + break; - /** - --> by starting from phi_small_epsilon and integrating an approximate - solution backward in time, try to estimate roughly a value close - to phi_pivot but a bit smaller. This is done by trying to reach - an amount of inflation equal to the requested one, minus the - amount after phi_small_epsilon, and plus - primordial_inflation_extra_efolds efolds (default: two). Note - that it is not aggressive to require two extra e-folds of - inflation before the pivot, since the calculation of the spectrum - in the observable range will require even more. */ - - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= phi_small_epsilon; - - switch (ppm->phi_pivot_method) { - - case ln_aH_ratio_auto: - case ln_aH_ratio: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _aH_, - H_small_epsilon/exp(target+ppr->primordial_inflation_extra_efolds)*aH_ratio_after_small_epsilon, - _TRUE_, - backward, - conformal), - ppm->error_message, - ppm->error_message); - break; + case ln_aH_ratio: - case N_star: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _a_, - 1./exp(target+ppr->primordial_inflation_extra_efolds)*a_ratio_after_small_epsilon, - _TRUE_, - backward, - conformal), - ppm->error_message, - ppm->error_message); - break; - } + target = ppm->phi_pivot_target; + //fprintf(stderr,"fixed: log(aH_end/aH_*)=%e\n",target); + break; - /* we now have a value phi_try believed to be close to and slightly smaller than phi_pivot */ + case N_star: - phi_try = y[ppm->index_in_phi]; + target = ppm->phi_pivot_target; + //fprintf(stderr,"fixed: log(a_end/a_*)=%e\n",target); + break; + } - /** - --> find attractor in phi_try */ + /** - --> by starting from phi_small_epsilon and integrating an approximate + solution backward in time, try to estimate roughly a value close + to phi_pivot but a bit smaller. This is done by trying to reach + an amount of inflation equal to the requested one, minus the + amount after phi_small_epsilon, and plus + primordial_inflation_extra_efolds efolds (default: two). Note + that it is not aggressive to require two extra e-folds of + inflation before the pivot, since the calculation of the spectrum + in the observable range will require even more. */ - class_call(primordial_inflation_find_attractor(ppm, - ppr, - phi_try, - ppr->primordial_inflation_attractor_precision_initial, - y, - dy, - &H_try, - &dphidt_try), - ppm->error_message, - ppm->error_message); + y[ppm->index_in_a]=1.; + y[ppm->index_in_phi]= phi_small_epsilon; - /** - --> check the total amount of inflation between phi_try and the end of inflation */ + switch (ppm->phi_pivot_method) { - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= phi_try; - y[ppm->index_in_dphi]= dphidt_try; + case ln_aH_ratio_auto: + case ln_aH_ratio: class_call(primordial_inflation_evolve_background(ppm, ppr, y, dy, - _end_inflation_, - 0., - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - - switch (ppm->phi_pivot_method) { - - case ln_aH_ratio_auto: - case ln_aH_ratio: + _aH_, + H_small_epsilon/exp(target+ppr->primordial_inflation_extra_efolds)*aH_ratio_after_small_epsilon, + _TRUE_, + backward, + conformal), + ppm->error_message, + ppm->error_message); + break; - // aH_ratio (we have used here proper time, so aH = dy[a]) - ratio_try = dy[ppm->index_in_a]/H_try; - break; + case N_star: - case N_star: + class_call(primordial_inflation_evolve_background(ppm, + ppr, + y, + dy, + _a_, + 1./exp(target+ppr->primordial_inflation_extra_efolds)*a_ratio_after_small_epsilon, + _TRUE_, + backward, + conformal), + ppm->error_message, + ppm->error_message); + break; + } - // a_ratio - ratio_try = y[ppm->index_in_a]; - break; - } + /* we now have a value phi_try believed to be close to and slightly smaller than phi_pivot */ - class_test(log(ratio_try) < target, - ppm->error_message, - "phi_try not small enough, log(aH_stop/aH_try) or log(a_stop/a_try) (depending on what you asked) is equal to %e instead of requested %e; must write here a loop to deal automatically with this situation (by decreasing phi_try iteratively), or must increase precision parameter primordial_inflation_extra_efolds", - log(ratio_try), - target); + phi_try = y[ppm->index_in_phi]; - phi_stop = y[1]; + /** - --> find attractor in phi_try */ - if (ppm->primordial_verbose > 1) - printf(" (inflation stops in phi_stop = %e)\n",phi_stop); - - /** - --> go back to phi_try, and now find phi_pivot such that the amount - of inflation between phi_pivot and the end of inflation is - exactly the one requested. */ - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= phi_try; - y[ppm->index_in_dphi]= dphidt_try; - - switch (ppm->phi_pivot_method) { - - case ln_aH_ratio_auto: - case ln_aH_ratio: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _aH_, - H_try*ratio_try/exp(target), - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - break; + class_call(primordial_inflation_find_attractor(ppm, + ppr, + phi_try, + ppr->primordial_inflation_attractor_precision_initial, + y, + dy, + &H_try, + &dphidt_try), + ppm->error_message, + ppm->error_message); - case N_star: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _a_, - ratio_try/exp(target), - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - break; - } + /** - --> check the total amount of inflation between phi_try and the end of inflation */ - ppm->phi_pivot = y[1]; + y[ppm->index_in_a]=1.; + y[ppm->index_in_phi]= phi_try; + y[ppm->index_in_dphi]= dphidt_try; - if (ppm->primordial_verbose > 1) { + class_call(primordial_inflation_evolve_background(ppm, + ppr, + y, + dy, + _end_inflation_, + 0., + _FALSE_, + forward, + proper), + ppm->error_message, + ppm->error_message); - printf(" (reached phi_pivot=%e)\n",ppm->phi_pivot); + switch (ppm->phi_pivot_method) { - /* - --> In verbose mode, check that phi_pivot is correct. Done by - restarting from phi_pivot and going again till the end of - inflation. */ + case ln_aH_ratio_auto: + case ln_aH_ratio: - aH_pivot = dy[0]; - a_pivot = y[0]; - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _end_inflation_, - 0., - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - printf(" (from phi_pivot till the end, ln(aH_2/aH_1) = %e, ln(a_2/a_1) = %e)\n",log(dy[0]/aH_pivot),log(y[0]/a_pivot)); - } + // aH_ratio (we have used here proper time, so aH = dy[a]) + ratio_try = dy[ppm->index_in_a]/H_try; + break; + case N_star: + // a_ratio + ratio_try = y[ppm->index_in_a]; + break; } - /** - case in which epsilon<1: */ - else { - /** - --> find inflationary attractor in phi_small_epsilon (should exist since epsilon<1 there) */ - class_call(primordial_inflation_find_attractor(ppm, - ppr, - ppm->phi_end, - ppr->primordial_inflation_attractor_precision_initial, - y, - dy, - &H_small_epsilon, - &dphidt_small_epsilon), - ppm->error_message, - ppm->error_message); + class_test(log(ratio_try) < target, + ppm->error_message, + "phi_try not small enough, log(aH_stop/aH_try) or log(a_stop/a_try) (depending on what you asked) is equal to %e instead of requested %e; must write here a loop to deal automatically with this situation (by decreasing phi_try iteratively), or must increase precision parameter primordial_inflation_extra_efolds", + log(ratio_try), + target); - /** - --> by starting from phi_end and integrating an approximate - solution backward in time, try to estimate roughly a value close - to phi_pivot but a bit smaller. This is done by trying to reach - an amount of inflation equal to the requested one, minus the - amount after phi_small_epsilon, and plus - primordial_inflation_extra_efolds efolds (default: two). Note - that it is not aggressive to require two extra e-folds of - inflation before the pivot, since the calculation of the spectrum - in the observable range will require even more. */ - - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= ppm->phi_end; - - switch (ppm->phi_pivot_method) { - - case ln_aH_ratio_auto: - case ln_aH_ratio: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _aH_, - H_small_epsilon/exp(target+ppr->primordial_inflation_extra_efolds)*aH_ratio_after_small_epsilon, - _TRUE_, - backward, - conformal), - ppm->error_message, - ppm->error_message); - break; + phi_stop = y[1]; - case N_star: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _a_, - 1./exp(target+ppr->primordial_inflation_extra_efolds)*a_ratio_after_small_epsilon, - _TRUE_, - backward, - conformal), - ppm->error_message, - ppm->error_message); - break; - } - - /** - --> we now have a value phi_try believed to be close to and slightly smaller than phi_pivot */ + if (ppm->primordial_verbose > 1) + printf(" (inflation stops in phi_stop = %e)\n",phi_stop); - phi_try = y[ppm->index_in_phi]; + /** - --> go back to phi_try, and now find phi_pivot such that the amount + of inflation between phi_pivot and the end of inflation is + exactly the one requested. */ + y[ppm->index_in_a]=1.; + y[ppm->index_in_phi]= phi_try; + y[ppm->index_in_dphi]= dphidt_try; - /** - --> find attractor in phi_try */ + switch (ppm->phi_pivot_method) { - class_call(primordial_inflation_find_attractor(ppm, - ppr, - phi_try, - ppr->primordial_inflation_attractor_precision_initial, - y, - dy, - &H_try, - &dphidt_try), - ppm->error_message, - ppm->error_message); - - /** - --> check the total amount of inflation between phi_try and the end of inflation */ - - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= phi_try; - y[ppm->index_in_dphi]= dphidt_try; + case ln_aH_ratio_auto: + case ln_aH_ratio: class_call(primordial_inflation_evolve_background(ppm, ppr, y, dy, - _phi_, - ppm->phi_end, + _aH_, + H_try*ratio_try/exp(target), _FALSE_, forward, proper), - ppm->error_message, - ppm->error_message); - - switch (ppm->phi_pivot_method) { - - case ln_aH_ratio_auto: - case ln_aH_ratio: - - // aH_ratio (we have used here proper time, so aH = dy[a]) - ratio_try = dy[ppm->index_in_a]/H_try; - break; - - case N_star: - - // a_ratio - ratio_try = y[ppm->index_in_a]; - break; - } - - class_test(log(ratio_try) < target, - ppm->error_message, - "phi_try not small enough, log(aH_stop/aH_try) or log(a_stop/a_try) (depending on what you asked) is equal to %e instead of requested %e; must write here a loop to deal automatically with this situation (by decreasing phi_try iteratively), or must increase precision parameter primordial_inflation_extra_efolds", - log(ratio_try), - target); + ppm->error_message, + ppm->error_message); + break; - phi_stop = y[1]; + case N_star: - if (ppm->primordial_verbose > 1) - printf(" (inflation stops in phi_stop = %e)\n",phi_stop); - - /** - --> go back to phi_try, and now find phi_pivot such that the amount - of inflation between phi_pivot and the end of inflation is - exactly the one requested. */ - y[ppm->index_in_a]=1.; - y[ppm->index_in_phi]= phi_try; - y[ppm->index_in_dphi]= dphidt_try; - - switch (ppm->phi_pivot_method) { - - case ln_aH_ratio_auto: - case ln_aH_ratio: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _aH_, - H_try*ratio_try/exp(target), - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - break; + class_call(primordial_inflation_evolve_background(ppm, + ppr, + y, + dy, + _a_, + ratio_try/exp(target), + _FALSE_, + forward, + proper), + ppm->error_message, + ppm->error_message); + break; + } - case N_star: - - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _a_, - ratio_try/exp(target), - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - break; - } + ppm->phi_pivot = y[1]; - ppm->phi_pivot = y[1]; + if (ppm->primordial_verbose > 1) { - if (ppm->primordial_verbose > 1) { + printf(" (reached phi_pivot=%e)\n",ppm->phi_pivot); - printf(" (reached phi_pivot=%e)\n",ppm->phi_pivot); + /* - --> In verbose mode, check that phi_pivot is correct. Done by + restarting from phi_pivot and going again till the end of + inflation. */ - /** - --> In verbose mode, check that phi_pivot is correct. Done by - restarting from phi_pivot and going again till the end of - inflation. */ + aH_pivot = dy[0]; + a_pivot = y[0]; + class_call(primordial_inflation_evolve_background(ppm, + ppr, + y, + dy, + _end_inflation_, + 0., + _FALSE_, + forward, + proper), + ppm->error_message, + ppm->error_message); + printf(" (from phi_pivot till the end, ln(aH_2/aH_1) = %e, ln(a_2/a_1) = %e)\n",log(dy[0]/aH_pivot),log(y[0]/a_pivot)); + } - aH_pivot = dy[0]; - a_pivot = y[0]; - class_call(primordial_inflation_evolve_background(ppm, - ppr, - y, - dy, - _phi_, - ppm->phi_end, - _FALSE_, - forward, - proper), - ppm->error_message, - ppm->error_message); - printf(" (from phi_pivot till the end, ln(aH_2/aH_1) = %e, ln(a_2/a_1) = %e)\n",log(dy[0]/aH_pivot),log(y[0]/a_pivot)); - } - } return _SUCCESS_; } From 551306cdfc346822be96f8e32e4af14f26364f72 Mon Sep 17 00:00:00 2001 From: Jobil Thomas <63110554+Jobilthomas@users.noreply.github.com> Date: Fri, 21 Apr 2023 21:13:15 +0530 Subject: [PATCH 2/4] Makefile --- Makefile | 206 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 206 insertions(+) create mode 100644 Makefile diff --git a/Makefile b/Makefile new file mode 100644 index 000000000..9e4a8a779 --- /dev/null +++ b/Makefile @@ -0,0 +1,206 @@ +#Some Makefile for CLASS. +#Julien Lesgourgues, 28.11.2011 +#Nils Schöneberg, Matteo Lucca, 27.02.2019 + +MDIR := $(shell pwd) +WRKDIR = $(MDIR)/build + +.base: + if ! [ -e $(WRKDIR) ]; then mkdir $(WRKDIR) ; mkdir $(WRKDIR)/lib; fi; + touch build/.base + +vpath %.c source:tools:main:test +vpath %.o build +vpath .base build + +######################################################## +###### LINES TO ADAPT TO YOUR PLATEFORM ################ +######################################################## + +# your C compiler: +CC = gcc +#CC = icc +#CC = pgcc + +# your tool for creating static libraries: +AR = ar rv + +# Your python interpreter. +# In order to use Python 3, you can manually +# substitute python3 to python in the line below, or you can simply +# add a compilation option on the terminal command line: +# "PYTHON=python3 make all" (Thanks to Marius Millea for python3 compatibility) +PYTHON ?= python + +# your optimization flag +OPTFLAG = -O3 +#OPTFLAG = -Ofast -ffast-math #-march=native +#OPTFLAG = -fast + +# your openmp flag (comment for compiling without openmp) +OMPFLAG = -fopenmp +#OMPFLAG = -mp -mp=nonuma -mp=allcores -g +#OMPFLAG = -openmp + +# all other compilation flags +CCFLAG = -g -fPIC +LDFLAG = -g -fPIC + +# leave blank to compile without HyRec, or put path to HyRec directory +# (with no slash at the end: e.g. "external/RecfastCLASS") +HYREC = external/HyRec2020 +RECFAST = external/RecfastCLASS +HEATING = external/heating + +######################################################## +###### IN PRINCIPLE THE REST SHOULD BE LEFT UNCHANGED ## +######################################################## + +# pass current working directory to the code +CCFLAG += -D__CLASSDIR__='"$(MDIR)"' + +# where to find include files *.h +INCLUDES = -I../include +HEADERFILES = $(wildcard ./include/*.h) + +# automatically add external programs if needed. First, initialize to blank. +EXTERNAL = + +vpath %.c $(RECFAST) +#CCFLAG += -DRECFAST +INCLUDES += -I../$(RECFAST) +EXTERNAL += wrap_recfast.o +HEADERFILES += $(wildcard ./$(RECFAST)/*.h) + +vpath %.c $(HEATING) +#CCFLAG += -DHEATING +INCLUDES += -I../$(HEATING) +EXTERNAL += injection.o noninjection.o +HEADERFILES += $(wildcard ./$(HEATING)/*.h) + +# update flags for including HyRec +ifneq ($(HYREC),) +vpath %.c $(HYREC) +CCFLAG += -DHYREC +#LDFLAGS += -DHYREC +INCLUDES += -I../$(HYREC) +EXTERNAL += hyrectools.o helium.o hydrogen.o history.o wrap_hyrec.o energy_injection.o +HEADERFILES += $(wildcard ./$(HYREC)/*.h) +endif + +%.o: %.c .base $(HEADERFILES) + cd $(WRKDIR);$(CC) $(OPTFLAG) $(OMPFLAG) $(CCFLAG) $(INCLUDES) -c ../$< -o $*.o + +TOOLS = growTable.o dei_rkck.o sparse.o evolver_rkck.o evolver_ndf15.o arrays.o parser.o quadrature.o hyperspherical.o common.o trigonometric_integrals.o + +SOURCE = input.o background.o thermodynamics.o perturbations.o primordial.o fourier.o transfer.o harmonic.o lensing.o distortions.o + +INPUT = input.o + +PRECISION = precision.o + +BACKGROUND = background.o + +THERMO = thermodynamics.o + +PERTURBATIONS = perturbations.o + +TRANSFER = transfer.o + +PRIMORDIAL = primordial.o + +HARMONIC = harmonic.o + +FOURIER = fourier.o + +LENSING = lensing.o + +DISTORTIONS = distortions.o + +OUTPUT = output.o + +CLASS = class.o + +TEST_LOOPS = test_loops.o + +TEST_LOOPS_OMP = test_loops_omp.o + +TEST_HARMONIC = test_harmonic.o + +TEST_TRANSFER = test_transfer.o + +TEST_FOURIER = test_fourier.o + +TEST_PERTURBATIONS = test_perturbations.o + +TEST_THERMODYNAMICS = test_thermodynamics.o + +TEST_BACKGROUND = test_background.o + +TEST_HYPERSPHERICAL = test_hyperspherical.o + +C_TOOLS = $(addprefix tools/, $(addsuffix .c,$(basename $(TOOLS)))) +C_SOURCE = $(addprefix source/, $(addsuffix .c,$(basename $(SOURCE) $(OUTPUT)))) +C_TEST = $(addprefix test/, $(addsuffix .c,$(basename $(TEST_DEGENERACY) $(TEST_LOOPS) $(TEST_TRANSFER) $(TEST_FOURIER) $(TEST_PERTURBATIONS) $(TEST_THERMODYNAMICS)))) +C_MAIN = $(addprefix main/, $(addsuffix .c,$(basename $(CLASS)))) +C_ALL = $(C_MAIN) $(C_TOOLS) $(C_SOURCE) +H_ALL = $(addprefix include/, common.h svnversion.h $(addsuffix .h, $(basename $(notdir $(C_ALL))))) +PRE_ALL = cl_ref.pre clt_permille.pre +INI_ALL = explanatory.ini lcdm.ini +MISC_FILES = Makefile CPU psd_FD_single.dat myselection.dat myevolution.dat README bbn/sBBN.dat external_Pk/* cpp +PYTHON_FILES = python/classy.pyx python/setup.py python/cclassy.pxd python/test_class.py + +all: class libclass.a classy + +libclass.a: $(TOOLS) $(SOURCE) $(EXTERNAL) + $(AR) $@ $(addprefix build/, $(TOOLS) $(SOURCE) $(EXTERNAL)) + +class: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(CLASS) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o class $(addprefix build/,$(notdir $^)) -lm + +test_loops: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(TEST_LOOPS) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_loops_omp: $(TOOLS) $(SOURCE) $(EXTERNAL) $(OUTPUT) $(TEST_LOOPS_OMP) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_harmonic: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_HARMONIC) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_transfer: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_TRANSFER) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_fourier: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_FOURIER) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_perturbations: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_PERTURBATIONS) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_thermodynamics: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_THERMODYNAMICS) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_background: $(TOOLS) $(SOURCE) $(EXTERNAL) $(TEST_BACKGROUND) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o $@ $(addprefix build/,$(notdir $^)) -lm + +test_hyperspherical: $(TOOLS) $(TEST_HYPERSPHERICAL) + $(CC) $(OPTFLAG) $(OMPFLAG) $(LDFLAG) -o test_hyperspherical $(addprefix build/,$(notdir $^)) -lm + + +tar: $(C_ALL) $(C_TEST) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES) + tar czvf class.tar.gz $(C_ALL) $(H_ALL) $(PRE_ALL) $(INI_ALL) $(MISC_FILES) $(HYREC) $(PYTHON_FILES) + +classy: libclass.a python/classy.pyx python/cclassy.pxd +ifdef OMPFLAG + cp python/setup.py python/autosetup.py +else + grep -v "lgomp" python/setup.py > python/autosetup.py +endif + cd python; export CC=$(CC); $(PYTHON) autosetup.py install || $(PYTHON) autosetup.py install --user + rm python/autosetup.py + +clean: .base + rm -rf $(WRKDIR); + rm -f libclass.a + rm -f $(MDIR)/python/classy.c + rm -rf $(MDIR)/python/build + rm -f python/autosetup.py From 87d3f73605142b0579810ad40219291a501d443f Mon Sep 17 00:00:00 2001 From: Jobilthomas Date: Sun, 30 Apr 2023 01:01:25 +0530 Subject: [PATCH 3/4] positive potential requirements --- include/primordial.h | 9 +++++++ source/input.c | 36 ++++++++++++++++++++++++++++ source/primordial.c | 56 ++++++++++++++++++++++++++++++++++++-------- 3 files changed, 91 insertions(+), 10 deletions(-) diff --git a/include/primordial.h b/include/primordial.h index 7c4fe9466..7e407e08d 100644 --- a/include/primordial.h +++ b/include/primordial.h @@ -69,6 +69,12 @@ enum inflation_module_behavior { analytical }; +/** enum specifying potential derivative dV/dphi sign */ + +enum potential_derivative_sign { + positive, + negative +}; /** * Structure containing everything about primordial spectra that other modules need to know. * @@ -181,6 +187,9 @@ struct primordial { enum phi_pivot_methods phi_pivot_method; /**< flag for method used to define and find the pivot scale */ double phi_pivot_target; /**< For each of the above methods, critical value to be reached between pivot and end of inflation (N_star, [aH]ratio, etc.) */ + /* parameter describing potential derivative dV/dphi sign */ + enum potential_derivative_sign potential_derivative; + /* behavior of the inflation module */ enum inflation_module_behavior behavior; /**< Specifies if the inflation module computes the primordial spectrum numerically (default) or analytically*/ diff --git a/source/input.c b/source/input.c index 6f3e3c85b..a12057a4c 100644 --- a/source/input.c +++ b/source/input.c @@ -4195,6 +4195,24 @@ int input_read_parameters_primordial(struct file_content * pfc, else if ((ppm->primordial_spec_type == inflation_V) || (ppm->primordial_spec_type == inflation_H)) { + + /* Sign of derivative of potential, dV/dphi */ + class_call(parser_read_string(pfc,"potential_derivative",&string1,&flag1,errmsg), + errmsg, + errmsg); + /* Complete set of parameters */ + if (flag1 == _TRUE_) { + if (strcmp(string1,"positive") == 0){ + ppm->potential_derivative = positive; + } + else if (strcmp(string1,"negative") == 0){ + ppm->potential = negative; + } + else{ + class_stop(errmsg,"You specified 'potential_derivative' as '%s'. It has to be one of {'positive','negative'}.",string1); + } + } + /** 1.c) For type 'inflation_V' */ if (ppm->primordial_spec_type == inflation_V) { @@ -4344,6 +4362,22 @@ int input_read_parameters_primordial(struct file_content * pfc, class_stop(errmsg,"You specified 'full_potential' as '%s'. It has to be one of {'polynomial','higgs_inflation'}.",string1); } } + /* Sign of derivative of potential, dV/dphi */ + class_call(parser_read_string(pfc,"potential_derivative",&string1,&flag1,errmsg), + errmsg, + errmsg); + /* Complete set of parameters */ + if (flag1 == _TRUE_) { + if (strcmp(string1,"positive") == 0){ + ppm->potential_derivative = positive; + } + else if (strcmp(string1,"negative") == 0){ + ppm->potential = negative; + } + else{ + class_stop(errmsg,"You specified 'potential_derivative' as '%s'. It has to be one of {'positive','negative'}.",string1); + } + } /** 1.e.3) Parameters of the potential */ /* Read */ @@ -5940,6 +5974,8 @@ int input_default_params(struct background *pba, ppm->phi_end=0.; /** 1.e.2) Shape of the potential */ ppm->potential=polynomial; + /* Sign of derivative of potential, dV/dphi */ + ppm->potential_derivative=negative; /** 1.e.4) Increase of scale factor or (aH) between Hubble crossing at pivot scale and end of inflation */ ppm->phi_pivot_method = N_star; diff --git a/source/primordial.c b/source/primordial.c index 0838656ce..fb08e7b81 100644 --- a/source/primordial.c +++ b/source/primordial.c @@ -2096,6 +2096,7 @@ int primordial_inflation_evolve_background( double quantity=0.; double V,dV,ddV; double sign_dtau=0.; + double sign_der=1.; pipaw.ppm = ppm; @@ -2167,6 +2168,18 @@ int primordial_inflation_evolve_background( break; case _phi_: // next (approximate) value of phi after next step + switch (ppm->potential_derivative){ + + case positive: + sign_der = -1.; + break; + + case negative: + sign_der = 1.; + break; + + } + quantity = y[ppm->index_in_phi]+dy[ppm->index_in_phi]*dtau; break; case _end_inflation_: @@ -2189,7 +2202,7 @@ int primordial_inflation_evolve_background( /* loop over time steps, checking that there will be no overshooting */ - while (sign_dtau*(quantity - stop) < 0.) { + while (sign_der*sign_dtau*(quantity - stop) < 0.) { /* check that V(phi) or H(phi) do not take forbidden values (negative or positive derivative) */ @@ -2419,12 +2432,22 @@ int primordial_inflation_check_potential( ppm->error_message, "This potential becomes negative at phi=%g, before the end of observable inflation. It cannot be treated by this code", phi); + switch (ppm->potential_derivative){ - class_test(*dV >= 0., + case positive: + class_test(*dV <= 0., ppm->error_message, - "All the code is written for the case dV/dphi<0. Here, in phi=%g, we have dV/dphi=%g. This potential cannot be treated by this code", + "The derivative of potential, dV/dphi<0. Here, in phi=%g, we have dV/dphi=%g.", phi,*dV); + break; + case negative: + class_test(*dV >= 0., + ppm->error_message, + "The derivative of potential, dV/dphi>0. Here, in phi=%g, we have dV/dphi=%g.", + phi,*dV); + break; + } return _SUCCESS_; } @@ -2563,6 +2586,7 @@ int primordial_inflation_find_phi_pivot( double sigma_B; double Omega_g0; double Omega_r0; + double sign_der=1.; /** - check whether in vicinity of phi_end, inflation is still ongoing */ @@ -2570,40 +2594,52 @@ int primordial_inflation_find_phi_pivot( ppm->error_message, ppm->error_message); + switch (ppm->potential_derivative){ + + case positive: + sign_der = -1.; + break; + + case negative: + sign_der = 1.; + break; + + } + // assume that inflation ends up naturally /** - --> find latest value of the field such that epsilon = primordial_inflation_small_epsilon (default: 0.1) */ - + if (epsilon > ppr->primordial_inflation_small_epsilon) { - /** - --> bracketing right-hand value is phi_end (but the potential will not be evaluated exactly there, only closeby */ + /** - --> bracketing right-hand value is phi_end (but the potential will not be evaluated exactly there, only closeby) */ phi_right = ppm->phi_end; /** - --> bracketing left-hand value is found by iterating with logarithmic step until epsilon < primordial_inflation_small_epsilon */ dphi = ppr->primordial_inflation_end_dphi; do { dphi *= ppr->primordial_inflation_end_logstep; - class_call(primordial_inflation_get_epsilon(ppm,ppm->phi_end-dphi,&epsilon), + class_call(primordial_inflation_get_epsilon(ppm,ppm->phi_end-sign_der*dphi,&epsilon), ppm->error_message, ppm->error_message); } while (epsilon > ppr->primordial_inflation_small_epsilon); - phi_left = ppm->phi_end-dphi; + phi_left = ppm->phi_end-sign_der*dphi; } else{ - /** - --> bracketing left-hand value is phi_end (but the potential will not be evaluated exactly there, only closeby */ + /** - --> bracketing left-hand value is phi_end (but the potential will not be evaluated exactly there, only closeby) */ phi_left = ppm->phi_end; /** - --> bracketing right-hand value is found by iterating with logarithmic step until epsilon < primordial_inflation_small_epsilon */ dphi = ppr->primordial_inflation_end_dphi; do { dphi *= ppr->primordial_inflation_end_logstep; - class_call(primordial_inflation_get_epsilon(ppm,ppm->phi_end+dphi,&epsilon), + class_call(primordial_inflation_get_epsilon(ppm,ppm->phi_end+sign_der*dphi,&epsilon), ppm->error_message, ppm->error_message); } while (epsilon < ppr->primordial_inflation_small_epsilon); - phi_right = ppm->phi_end+dphi; + phi_right = ppm->phi_end+sign_der*dphi; } /** - --> find value such that epsilon = primordial_inflation_small_epsilon by bisection */ do { From 6c654e6c9576e34dd376f45836fa35b66b99e45d Mon Sep 17 00:00:00 2001 From: Jobil Thomas <63110554+Jobilthomas@users.noreply.github.com> Date: Sat, 29 Jul 2023 20:41:41 +0530 Subject: [PATCH 4/4] input.c - bugfix bug fix. changed ppm->potential to ppm->potential_derivative in line 4209 of input.c --- source/input.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/input.c b/source/input.c index a12057a4c..ef4dc0696 100644 --- a/source/input.c +++ b/source/input.c @@ -4206,7 +4206,7 @@ int input_read_parameters_primordial(struct file_content * pfc, ppm->potential_derivative = positive; } else if (strcmp(string1,"negative") == 0){ - ppm->potential = negative; + ppm->potential_derivative = negative; } else{ class_stop(errmsg,"You specified 'potential_derivative' as '%s'. It has to be one of {'positive','negative'}.",string1);