diff --git a/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/README.txt b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/README.txt new file mode 100644 index 0000000000..5ee98ee450 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/README.txt @@ -0,0 +1,81 @@ +The procedures of this LSS simulation includes: +--------------------------------------------------------------------------------------------------------------- + Flow chart + if not axion transfer function if axion transfer function + | | + | "planck_2018.ini" | "planck_2018_axion.ini" + | | + v v + [(a) CAMB] [(a) axionCAMB] + [camb planck_2018.ini] [./camb planck_2018_axion.ini] + | | + | "planck_2018_transfer_out.dat" | "planck_2018_axion_transfer_out.dat" + | | + | | + -------> [(b) set_velocities_zero.py] <----- + [python3 ./set_velocities_zero.py F planck_2018_trnasfer_out.dat planck_2018_transfer_out_no_vel.dat (CAMB)] + [python3 ./set_velocities_zero.py T planck_2018_axion_trnasfer_out.dat planck_2018_axion_transfer_out_no_vel.dat (axionCAMB)] + | + | "planck_2018_transfer_out_no_vel.dat" (CAMB) + | "planck_2018_axion_transfer_out_no_vel.dat" (axionCAMB) + | + "ics_example.conf" v + ----------------------> [(c) MUSIC] + [./MUSIC ics_example.conf] + | + | "planck_2018_tf_no_vel.hdf5" (CAMB) + | "planck_2018_axion_tf_no_vel.hdf5" (axionCAMB) + | + "ics_example.conf" v + --------------> [(d) make_umic_from_hdf5.py] + [python3 ./make_umic_from_hdf5.py S(D) planck_2018_tf_no_vel.hdf (CAMB)] + [python3 ./make_umic_from_hdf5.py S(D) planck_2018_axion_tf_no_vel.hdf (axionCAMB)] + | + | "UM_IC" + | + v + [(e) GAMER] + +--------------------------------------------------------------------------------------------------------------- + +(a) After successful installation of CAMB/axionCAMB (see Notes for installation), use CAMB/axionCAMB to to produce transfer function. + The input file planck_2018.ini/planck_2018_axion.ini for CAMB/axionCAMB is attached for reference. + The expected output files are: + planck_2018_trnasfer_out.dat -> for CAMB + planck_2018_axion_trnasfer_out.dat -> for axionCAMB + +(b) Run the script: "python3 ./set_velocities_zero.py F(T) planck_2018_trnasfer_out.dat(planck_2018_axion_trnasfer_out.dat) planck_2018_transfer_out_no_vel.dat(planck_2018_axion_transfer_out_no_vel.dat)". + The script reads-in the output files in step(a), and produce the transfer function requireded by MUSIC. + The boolean F/T suggests whether axion-physics is included in the input transfer function files. + The expected output files are: + planck_2018_transfer_out_no_vel.dat -> for CAMB + planck_2018_axion_transfer_out_no_vel.dat -> for axionCAMB + +(c) Modifying the cosmology paramters to be consistent with CAMB/axionCAMB .ini files, as well as the input file parameters in MUSIC input file ics_example.conf: + transfer = camb_file + transfer_file = planck_2018_transfer_out_no_vel.dat/planck_2018_axion_transfer_out_no_vel.dat + format = generic + filename = planck_2018_tf_no_vel.hdf5/planck_2018_axion_tf_no_vel.hdf5) + Install (see Notes for installation) and run the MUSIC: "./MUSIC ics_example.conf". + +(d) Run the script "python3 ./make_umic_from_hdf5.py S(D) planck_2018_axion_tf_no_vel.hdf5" to produce UM_IC file for single(double) precision UM_IC for FDM simulation (or "python3 ./make_umic_from_hdf5.py S(D) planck_2018_tf_no_vel.hdf5" for CDM simulation). + +(e) Run the GAMER. + +# Notes + 1. CAMB will automatically compute \sigma_8 based on the given cosomology parameters. + If you find the given \sigma_8 is not the desired value, it can be tuned by changing the As (scalar_amp in the input file .ini). + Say if one wants to change the original \sigma_8 = 0.8119112 with As = 2.100549e-9 , to new value \sigma_8 = 0.818 , one just changes the new As to (0.818/0.8119112)^2*2.100549e-9 = 2.132173e-09 , then CAMB will give new \sigma_8 = 0.818. + This is due to power spectrum is proportional to As*f(k) . + 2. CAMB (Code for Anisotropies in the Microwave Background) + GitHub page: https://github.com/cmbant/CAMB + Install procedures: https://github.com/cmbant/CAMB#description-and-installation + 3. axionCAMB(modified version of CAMB to include axion physics) + GitHub page: https://github.com/dgrin1/axionCAMB + Install procedures: https://github.com/dgrin1/axionCAMB?tab=readme-ov-file#compiling-axioncamb + 4. MUSIC(MUlti-Scale Initial Conditions for cosmological simulations) + Bitbucket page: https://bitbucket.org/ohahn/music/src/master/ + Website (for user's guide pdf): https://www-n.oca.eu/ohahn/MUSIC/ + Install procedures: modified the path and parameters in MUSIC-repo-provided Makefile, then make + or, can refer to the section: "Building MUSIC" in the Bitbucket page + diff --git a/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/ics_example.conf b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/ics_example.conf new file mode 100644 index 0000000000..13e920951f --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/ics_example.conf @@ -0,0 +1,79 @@ +[setup] +boxlength = 1.4 +zstart = 100 +levelmin = 10 +levelmin_TF = 10 +levelmax = 10 +padding = 4 +overlap = 1 +ref_center = 0.5, 0.5, 0.5 +ref_extent = 0.2, 0.2, 0.2 +align_top = no +baryons = no +use_2LPT = yes +use_LLA = no +periodic_TF = yes + + +[cosmology] +Omega_m = 0.3158230904284232 +Omega_L = 0.6841769095715768 +w0 = -1.0 +wa = 0.0 +Omega_b = 0.04938682464547351 +H0 = 67.32117 +sigma_8 = 0.81191119 +nspec = 0.96605 +transfer = camb_file +#transfer_file = planck_2018_transfer_out_no_vel.dat +transfer_file = planck_2018_axion_transfer_out_no_vel.dat + + +[random] +seed[7] = 12345 +seed[8] = 23456 +seed[9] = 34567 +seed[10] = 45678 +seed[11] = 56789 +seed[12] = 67890 + + +[output] +##generic MUSIC data format (used for testing) +##requires HDF5 installation and HDF5 enabled in Makefile +format = generic +filename = planck_2018_axion_tf_no_vel.hdf5 +#filename = planck_2018_tf_no_vel.hdf5 + +##ENZO - also outputs the settings for the parameter file +##requires HDF5 installation and HDF5 enabled in Makefile +#format = enzo +#filename = ic.enzo + +##Gadget-2 (type=1: high-res particles, type=5: rest) +#format = gadget2 +#filename = ics_gadget.dat + +##Grafic2 compatible format for use with RAMSES +##option 'ramses_nml'=yes writes out a startup nml file +#format = grafic2 +#filename = ics_ramses +#ramses_nml = yes + +##TIPSY compatible with PKDgrav and Gasoline +#format = tipsy +#filename = ics_tipsy.dat + +## NYX compatible output format +##requires boxlib installation and boxlib enabled in Makefile +#format = nyx +#filename = init + +[poisson] +fft_fine = yes +accuracy = 1e-5 +pre_smooth = 3 +post_smooth = 3 +smoother = gs +laplace_order = 6 +grad_order = 6 diff --git a/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/make_umic_from_hdf5.py b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/make_umic_from_hdf5.py new file mode 100755 index 0000000000..8c68adf752 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/make_umic_from_hdf5.py @@ -0,0 +1,156 @@ +############################################################################################# +# This script is used for genertating the wave function needed for GAMER with ELBDM from # +# MUSIC, included its real and imaginary parts. Be sure the requirements below are met: # +# (1) Input configuration file (ics_example.conf here) for MUSIC is under the same folder. # +# (2) Output generic hdf5 file (generated by the above input file)is under the same folder. # +# (3) Adjust the psidm_mass in physical constant if needed. # +# # +# Use handle S/D to select output UM_IC as single/double precision and assign the hdf5 file # +# when running this script, e.g. ./make_umic_from_hdf5.py S(D) hdf5_filename # +# # +# The code does the follwing things: # +# (1) Scan through the input file to collect parameters # +# (2) Calculate the density from over-density given by MUSIC # +# (3) Calculate the phase from velocity by solving the poisson equation in spcetrum way: # +# a*psidm_mass*(div(velocity))/h_bar = del(phase) # +# div is evaluated via Richardson extrapolation with an adjustable order # +# del is the \nabla operator, i.e., del(phase) is the phase gradient # +# (4) Convert density and phase field to real and imaginary part of wave function # +# (5) Output as binary file "UM_IC" # +# # +# Unit(s): MUSIC v.s. here # +# (1) density: back ground density ; back ground density # +# (2) velocity: box length*H_0 ; m/s # +# # +# Physical meaning for MUSIC's velocity recorded in HDF5 is peculiar velocity (in physical # +# frame, not comoving frame), i.e., after multiplying box_length*km_to_m, we have: # +# (v_hdf5),i = (v_peculiar,physical),i = a*(dx_{comoving}),i/dt # +# where i=x, y, z # +############################################################################################# + +import os, h5py, subprocess, re, sys +import numpy as np + +# Calculate the gradient by Richardson extrapolation (with periodic boundary condition) +def GRAD(field, axis, h, order): + dim = list(field.shape) + dim.insert(0, order) + grad = np.zeros(tuple(dim)) + for o in range(order): + interval = 2**(order-1-o) + grad[o] = (np.roll(field, -interval, axis=axis) - np.roll(field, interval, axis=axis)) / (2*interval) + for o in range(1,order): + grad[o:] = (4.**o*grad[o:]-grad[o-1:-1]) / (4.**o-1.) + return grad[-1]/h +# Physical Constant +psidm_mass = 8e-23 #eV +h_bar = 1.0545718e-34 +eV_to_kg = 1.78266192162790e-36 +Mpc_to_m = 3.0856776e22 +km_to_m = 1.0e3 +##################################################################################################################### + +if len(sys.argv)!=3: + raise RuntimeError("Input Error: Please select output UMIC as single/double precision (S/D) and file name of hdf5!") +else: + input_file = "ics_example.conf" + input_hdf5 = sys.argv[2] + if sys.argv[1] == "S": + print("Output UMIC as single precision.") + flag_D = False + elif sys.argv[1] == "D": + print("Output UMIC as double precision.") + flag_D = True + else: + raise RuntimeError("Input Error: Please select output UMIC as single/double precision (S/D)!") + + filename_hdf5 = sys.argv[2] + if os.path.isfile("./%s"%filename_hdf5): + print("hdf5 file name extracted successfully! Filename is %s ."%filename_hdf5) + else: + raise RuntimeError("File %s cannot be found!"%filename_hdf5) + + + output_file = "UM_IC" + ghost_zone = 4 + + cmd = "sed -e '/^#/d' %s | grep levelmax | awk {'print $3'}"%input_file + level = subprocess.check_output(cmd, shell=True).decode('utf-8') + level = eval(re.sub("\n","", level)) + print("Maximum level is %d ."%level) + + cmd = "grep boxlength %s | awk '{printf $3}' "%input_file + box_length = subprocess.check_output(cmd, shell=True).decode('utf-8') + box_length = eval(re.sub("\n","", box_length)) + print("Box length is %.4f Mpc/h."%box_length) + + cmd = "grep H0 %s | awk '{printf $3}' "%input_file + H0 = subprocess.check_output(cmd, shell=True).decode('utf-8') + H0 = eval(re.sub("\n","", H0)) + print("H0 is %.4f km/s/Mpc."%H0) + + cmd = "grep zstart %s | awk '{printf $3}' "%input_file + z = subprocess.check_output(cmd, shell=True).decode('utf-8') + z = eval(re.sub("\n","", z)) + print("z_start is %.2f ."%z) + + factor = box_length*km_to_m + N = 2**level + h = 1./N + k_factor = 2.*np.pi/N + box_length *= Mpc_to_m/(H0/100.) # meter + a = 1./(1.+z) + + with h5py.File(input_hdf5,'r') as f: + density_hdf5 = f['level_%03d_DM_rho'%level][ghost_zone:-ghost_zone, ghost_zone:-ghost_zone, ghost_zone:-ghost_zone] + vx_hdf5 = f['level_%03d_DM_vx'%level ][ghost_zone:-ghost_zone, ghost_zone:-ghost_zone, ghost_zone:-ghost_zone] + vy_hdf5 = f['level_%03d_DM_vy'%level ][ghost_zone:-ghost_zone, ghost_zone:-ghost_zone, ghost_zone:-ghost_zone] + vz_hdf5 = f['level_%03d_DM_vz'%level ][ghost_zone:-ghost_zone, ghost_zone:-ghost_zone, ghost_zone:-ghost_zone] + f.close() + + + criteria = (density_hdf5<-1.) + print("Percentage of over-density smaller than -1: %.8f %%."%(100*criteria.sum()/2**(3*level)) ) + density_hdf5[criteria] = -1. + vx_hdf5[criteria] = 0. + vy_hdf5[criteria] = 0. + vz_hdf5[criteria] = 0. + + density_hdf5 = (density_hdf5+1.)**0.5 + vx_hdf5 *= factor + vy_hdf5 *= factor + vz_hdf5 *= factor + + # Calculate div(v) + vx_x = GRAD(vx_hdf5, axis = 0, h=h ,order=3) + vy_y = GRAD(vy_hdf5, axis = 1, h=h ,order=3) + vz_z = GRAD(vz_hdf5, axis = 2, h=h ,order=3) + v_div = vx_x + vy_y + vz_z + v_div *= a*psidm_mass*eV_to_kg/h_bar + + # Do forward DFT + v_div_k = np.fft.rfftn(v_div) + # Do inverse Laplacian + kx, ky, kz = np.arange(N), np.arange(N), np.arange(N//2+1.) + kxx, kyy, kzz = np.meshgrid(kx, ky, kz) + v_div_k /= 2.*(np.cos(k_factor*kxx)+np.cos(k_factor*kyy)+np.cos(k_factor*kzz)-3.) + v_div_k[0,0,0] = 0. + # Do inverse DFT + phase = np.fft.irfftn(v_div_k) + # Rescale to correct unit + phase *= box_length/N**2 + phase -= phase.min() + + # Convert to wave function + wf_real = density_hdf5*np.cos(phase) + wf_imag = density_hdf5*np.sin(phase) + + new_data_hdf5 = np.zeros((2, N, N, N)) + new_data_hdf5[0] = np.swapaxes(wf_real,0,2) + new_data_hdf5[1] = np.swapaxes(wf_imag,0,2) + if not flag_D: + new_data_hdf5 = new_data_hdf5.astype(np.float32) + print("Writing wave function to binary file...") + with open(output_file,"wb") as f: + new_data_hdf5.tofile(f) + f.close() diff --git a/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018.ini b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018.ini new file mode 100644 index 0000000000..fd7efef592 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018.ini @@ -0,0 +1,283 @@ +#Parameters for CAMB +#Here, best-fit parameters from the baseline Planck 2018 LCDM model +#using base_plikHM_TTTEEE_lowl_lowE_lensing +# +#Accuracy settings are default for Planck, i.e. not accurate for lensing at high L +#For high accuracy to high L see planck_2018_base_plikHM_TTTEEE_lowl_lowE_lensing_lmax10000.ini + +#output_root is prefixed to output file names +output_root = planck_2018 + +#What to do +get_scalar_cls = T +get_vector_cls = F +get_tensor_cls = F +get_transfer = T + +#if do_lensing then lens_potential_output_file contains the unlensed CMB and lensing potential power spectra +#and lensed CMB Cls (without tensors) are in lensed_output_file, total in lensed_total_output_file. +do_lensing = T + +# 0: linear, 1: non-linear matter power (HALOFIT), 2: non-linear CMB lensing (HALOFIT), +# 3: both non-linear matter power and CMB lensing (HALOFIT) +do_nonlinear = 2 + +#Maximum multipole and k*eta. +# Note that C_ls near l_max are inaccurate (about 5%), go to 50 more than you need +# Lensed power spectra are computed to l_max_scalar-100 +# To get accurate lensed BB need to have l_max_scalar>2000, k_eta_max_scalar > 10000 +# To get accurate lensing potential you also need k_eta_max_scalar > 10000 +# Otherwise k_eta_max_scalar=2*l_max_scalar usually suffices, or don't set to use default +l_max_scalar = 2700 +k_eta_max_scalar = 18000 + +# Tensor settings should be less than or equal to the above +l_max_tensor = 1500 +k_eta_max_tensor = 3000 + +#Main cosmological parameters +#Set physical densities in baryons, CDM and neutrinos + Omega_k +ombh2 = 0.0223828 +omnuh2 = 0.6451439E-03 +omch2 = 0.1201075 +omk = 0 +hubble = 67.32117 + +#Dark energy can be "fluid" (constant w), or "PPF" (allowing wa, and w=-1 crossing) +dark_energy_model = ppf + +#effective equation of state parameter for dark energy +w = -1 +#constant comoving sound speed of the dark energy (1=quintessence), for fluid model +cs2_lam = 1 + +#For PPF model can set these instead: +#wa = 0 +##if use_tabulated_w read (a,w) from the following user-supplied file instead of above +#use_tabulated_w = T +#wafile = wa.dat + +temp_cmb = 2.7255 +helium_fraction = 0.2454006 + +#for share_delta_neff = T, the fractional part of massless_neutrinos gives the change in the effective number +#(for QED + non-instantaneous decoupling) i.e. the increase in neutrino temperature, +#so Neff = massless_neutrinos + sum(massive_neutrinos) +#For full neutrino parameter details see https://cosmologist.info/notes/CAMB.pdf +massless_neutrinos = 2.046 + +#number of distinct mass eigenstates +nu_mass_eigenstates = 1 +#array of the integer number of physical neutrinos per eigenstate, e.g. massive_neutrinos = 2 1 +massive_neutrinos = 1 +#specify whether all neutrinos should have the same temperature, specified from fractional part of massless_neutrinos +share_delta_neff = T +#nu_mass_fractions specifies how Omeganu_h2 is shared between the eigenstates +#i.e. to indirectly specify the mass of each state; e.g. nu_mass_factions= 0.75 0.25 +nu_mass_fractions = 1 +#if share_delta_neff = F, specify explicitly the degeneracy for each state (e.g. for sterile with different temperature to active) +#(massless_neutrinos must be set to degeneracy for massless, i.e. massless_neutrinos does then not include Deleta_Neff from massive) +#if share_delta_neff=T then degeneracies is not given and set internally +#e.g. for massive_neutrinos = 2 1, this gives equal temperature to 4 neutrinos: nu_mass_degeneracies = 2.030 1.015, massless_neutrinos = 1.015 +nu_mass_degeneracies = + +#Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}. +initial_power_num = 1 +pivot_scalar = 0.05 +pivot_tensor = 0.05 +scalar_amp = 2.100549e-9 +scalar_spectral_index = 0.96605 +scalar_nrun = 0 +scalar_nrunrun = 0 +tensor_spectral_index = 0 +tensor_nrun = 0 +#Three parameterizations (1,2,3) for tensors, see https://cosmologist.info/notes/CAMB.pdf +tensor_parameterization = 1 +#ratio is that of the initial tens/scal power spectrum amplitudes, depending on parameterization +#for tensor_parameterization == 1, P_T = initial_ratio*scalar_amp*(k/pivot_tensor)^tensor_spectral_index +#for tensor_parameterization == 2, P_T = initial_ratio*P_s(pivot_tensor)*(k/pivot_tensor)^tensor_spectral_index +#Note that for general pivot scales and indices, tensor_parameterization==2 has P_T depending on n_s +initial_ratio = 1 +#tensor_amp is used instead if tensor_parameterization == 3, P_T = tensor_amp *(k/pivot_tensor)^tensor_spectral_index +#tensor_amp = 4e-10 + +#note vector modes use the scalar settings above + + +#Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5 +reionization = T + +re_use_optical_depth = T +re_optical_depth = 0.05430842 +#If re_use_optical_depth = F then use following, otherwise ignored +re_redshift = 11 +#width of reionization transition. CMBFAST model was similar to re_delta_redshift~0.5. +re_delta_redshift = 0.5 +#re_ionization_frac=-1 sets it to become fully ionized using Yhe to get helium contribution +#Otherwise x_e varies from 0 to re_ionization_frac +re_ionization_frac = -1 + +#Parameters for second reionization of helium +re_helium_redshift = 3.5 +re_helium_delta_redshift = 0.5 + +#Can compile with HyRec and/or CosmoRec support as well; Recfast is default +recombination_model = Recfast +#RECFAST 1.5.x recombination parameters; +RECFAST_fudge = 1.14 +RECFAST_fudge_He = 0.86 +RECFAST_Heswitch = 6 +RECFAST_Hswitch = T + +# compile with RECOMBINATION_FILES including cosmorec +# +# cosmorec_runmode== 0: with diffusion +# 1: without diffusion +# 2: RECFAST++ run (equivalent of the original RECFAST version) +# 3: RECFAST++ run with correction function of Calumba & Thomas, 2010 +# +# For 'cosmorec_accuracy' and 'cosmorec_fdm' see CosmoMC for explanation +#--------------------------------------------------------------------------------------- +#cosmorec_runmode = 0 +#cosmorec_accuracy = 0 +#cosmorec_fdm = 0 + +#Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, +# neutrino density iso =4, neutrino velocity iso = 5) +initial_condition = 1 +#If above is zero, use modes in the following (totally correlated) proportions +#Note: we assume all modes have the same initial power spectrum +initial_vector = -1 0 0 0 0 + +#For vector modes: 0 for regular (neutrino vorticity mode), 1 for magnetic +vector_mode = 0 + +#Normalization +COBE_normalize = F +##CMB_outputscale scales the output Culs +#To get MuK^2 set realistic initial amplitude (e.g. scalar_amp(1) = 2.3e-9 above) and +#otherwise for dimensionless transfer functions set scalar_amp(1)=1 and use +CMB_outputscale = 1 +#CMB_outputscale = 7.42835025e12 + +#Transfer function settings, transfer_kmax=0.5 is enough for sigma_8 +#transfer_k_per_logint=0 sets sensible non-even sampling; +#transfer_k_per_logint=5 samples fixed spacing in log-k +#transfer_interp_matterpower =T produces matter power in regular interpolated grid in log k; +# use transfer_interp_matterpower =F to output calculated values (e.g. for later interpolation) +transfer_high_precision = T +transfer_kmax = 1000 +transfer_k_per_logint = 5 +transfer_num_redshifts = 1 +transfer_interp_matterpower = T +transfer_redshift(1) = 0.0000 +transfer_filename(1) = transfer_out.dat +#Matter power spectrum output against k/h in units of h^{-3} Mpc^3 +transfer_matterpower(1) = matterpower.dat + +#which variable to use for defining the matter power spectrum and sigma8 +#main choices are 2: CDM, 7: CDM+baryon+neutrino, 8: CDM+baryon, 9: CDM+baryon+neutrino+de perts +transfer_power_var = 7 + +#Output files not produced if blank. make camb_fits to use the FITS setting. +scalar_output_file = scalCls.dat +vector_output_file = vecCls.dat +tensor_output_file = tensCls.dat +total_output_file = totCls.dat +lensed_output_file = lensedCls.dat +lensed_total_output_file = lensedtotCls.dat +lens_potential_output_file = lenspotentialCls.dat +FITS_filename = scalCls.fits + +#Bispectrum parameters if required; primordial is currently only local model (fnl=1) +#lensing is fairly quick, primordial takes several minutes on quad core +do_lensing_bispectrum = F +do_primordial_bispectrum = F + +#1 for just temperature, 2 with E +bispectrum_nfields = 1 +#set slice non-zero to output slice b_{bispectrum_slice_base_L L L+delta} +bispectrum_slice_base_L = 0 +bispectrum_ndelta = 3 +bispectrum_delta(1) = 0 +bispectrum_delta(2) = 2 +bispectrum_delta(3) = 4 +#bispectrum_do_fisher estimates errors and correlations between bispectra +#note you need to compile with LAPACK and FISHER defined to use get the Fisher info +bispectrum_do_fisher = F +#Noise is in muK^2, e.g. 2e-4 roughly for Planck temperature +bispectrum_fisher_noise = 0 +bispectrum_fisher_noise_pol = 0 +bispectrum_fisher_fwhm_arcmin = 7 +#Filename if you want to write full reduced bispectrum (at sampled values of l_1) +bispectrum_full_output_file = +bispectrum_full_output_sparse = F +#Export alpha_l(r), beta_l(r) for local non-Gaussianity +bispectrum_export_alpha_beta = F + +##Optional parameters to control the computation speed,accuracy and feedback + +#If feedback_level > 0 print out useful information computed about the model +feedback_level = 1 + +#whether to start output files with comment describing columns +output_file_headers = T + +#write out various derived parameters +derived_parameters = T + +# 1: curved correlation function, 2: flat correlation function, 3: inaccurate harmonic method +lensing_method = 1 +accurate_BB = F + + +#massive_nu_approx: 0 - integrate distribution function +# 1 - switch to series in velocity weight once non-relativistic +massive_nu_approx = 1 + +#Whether you are bothered about polarization. +accurate_polarization = T + +#Whether you are bothered about percent accuracy on EE from reionization +accurate_reionization = T + +#whether or not to include neutrinos in the tensor evolution equations +do_tensor_neutrinos = T + +#whether you care about accuracy of the neutrino transfers themselves +accurate_massive_neutrino_transfers = F + +#Whether to turn off small-scale late time radiation hierarchies (save time,v. accurate) +do_late_rad_truncation = T + +#Which version of Halofit approximation to use (default currently Takahashi): +#1. Original Smith et al. (2003; arXiv:astro-ph/0207664) HALOFIT +#2. Bird et al. (arXiv:1109.4416) updated HALOFIT +#3. Original plus fudge from http://www.roe.ac.uk/~jap/haloes/, +#4. Takahashi (2012; arXiv:1208.2701) HALOFIT update +#5. HMcode (Mead et al. 2016; arXiv 1602.02154) +#6. A standard (inaccurate) halo model power spectrum calcultion +#7. PKequal (Casarini et al. arXiv:0810.0190, arXiv:1601.07230) +#8. HMcode (Mead et al. 2015; arXiv 1505.07833) +halofit_version = 5 + +#Computation parameters +#if number_of_threads=0 assigned automatically +number_of_threads = 0 + +#The default target accuracy is 0.1% at L>600 (with boost parameter=1 below) +#Try accuracy_boost=2, l_accuracy_boost=2 if you want to check stability/even higher accuracy +#Note increasing accuracy_boost parameters is very inefficient if you want higher accuracy, + +#Increase accuracy_boost to decrease time steps, use more k values, etc. +#Decrease to speed up at cost of worse accuracy. Suggest 0.8 to 3. +accuracy_boost = 2.1 + +#Larger to keep more terms in the hierarchy evolution. +l_accuracy_boost = 2.1 + +#Increase to use more C_l values for interpolation. +#Increasing a bit will improve the polarization accuracy at l up to 200 - +#interpolation errors may be up to 3% +#Decrease to speed up non-flat models a bit +l_sample_boost = 128 diff --git a/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018_axion.ini b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018_axion.ini new file mode 100644 index 0000000000..c542a248e8 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018_axion.ini @@ -0,0 +1,243 @@ +#Parameters for CAMB + +#output_root is prefixed to output file names +output_root = planck_2018_axion + +#What to do +get_scalar_cls = T +get_vector_cls = F +get_tensor_cls = F +get_transfer = T + +#if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT} +#where p is the projected potential. Output lensed CMB Cls (without tensors) are in lensed_output_file below. +do_lensing = T + +# 0: linear, 1: non-linear matter power (HALOFIT), 2: non-linear CMB lensing (HALOFIT) +do_nonlinear = 2 +halofit_version = 1 + +#Maximum multipole and k*eta. +# Note that C_ls near l_max are inaccurate (about 5%), go to 50 more than you need +# Lensed power spectra are computed to l_max_scalar-100 +# To get accurate lensed BB need to have l_max_scalar>2000, k_eta_max_scalar > 10000 +# Otherwise k_eta_max_scalar=2*l_max_scalar usually suffices, or don't set to use default +l_max_scalar = 2700 +k_eta_max_scalar = 18000 + +# Tensor settings should be less than or equal to the above +l_max_tensor = 1500 +k_eta_max_tensor = 3000 + +#Main cosmological parameters, neutrino masses are assumed degenerate +# If use_phyical set phyiscal densities in baryons, CDM and neutrinos + Omega_k +use_physical = T +ombh2 = 0.0223828 +#omch2 = .1201075 +omch2 = 0.0000001 +omnuh2 = 0.6451439e-3 +omk = 0 +hubble = 67.32117 +#effective equation of state parameter for dark energy, assumed constant +w = -1.0 +#constant comoving sound speed of the dark energy (1=quintessence) +cs2_lam = 1 + +#Axion variables +# Mass is in units eV +use_axfrac=T +#now we will ignore omaxh2 if use_axfrac=T +axion_isocurvature=F +alpha_ax = 0.0 +Hinf=13.7 #Log Hinflation in GeV + +#Axion variables +# Mass is in units eV +omaxh2 = 0.1201075 +#omaxh2 = 0.000001 +m_ax = 8e-23 + +# if use_axfrac = T set parameters as here, ignored if use_axfrac=F +omdah2 = 0.1201075 +axfrac = 1.0 +#omdah2 = 0.119 +#axfrac = 0.01 + +#if use_physical = F set parameters as here +#omega_baryon = 0.0462 +#omega_cdm = 0.2538 +#omega_lambda = 0.7 +#omega_neutrino = 0 +#omega_axion = 0 + + +temp_cmb = 2.7255 +helium_fraction = 0.2454006 +# massless_neutrinos is the effective number (for QED + non-instantaneous decoupling) +# fractional part of the number is used to increase the neutrino temperature, e.g. +# 2.99 correponds to 2 neutrinos with a much higher temperature, 3.04 correponds to +# 3 neutrinos with a slightly higher temperature +massless_neutrinos = 2.0406 +massive_neutrinos = 1 +share_delta_neff= T +#Neutrino mass splittings +nu_mass_eigenstates = 1 +#nu_mass_degeneracies = 0 sets nu_mass_degeneracies = massive_neutrinos +#otherwise should be an array +#e.g. for 3 neutrinos with 2 non-degenerate eigenstates, nu_mass_degeneracies = 2 1 +#nu_mass_degeneracies = 0 #overriden when share_delta_neff= true, otherwise set +#Fraction of total omega_nu h^2 accounted for by each eigenstate, eg. 0.5 0.5 +nu_mass_fractions = 1 + +#Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}. +initial_power_num = 1 +pivot_scalar = 0.05 +pivot_tensor = 0.05 +scalar_amp(1) = 2.100549e-9 +scalar_spectral_index(1) = 0.96605 +scalar_nrun(1) = 0 +tensor_spectral_index(1) = 0 +#ratio is that of the initial tens/scal power spectrum amplitudes +initial_ratio(1) = 0 +tens_ratio = 0 +#DM: gave this twice now because of problem using ini_driver to read initial_ratio +#note vector modes use the scalar settings above + + +#Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5 +reionization = T + +re_use_optical_depth = T +re_optical_depth = 0.05430842 +#If re_use_optical_depth = F then use following, otherwise ignored +re_redshift = 11 +#width of reionization transition. CMBFAST model was similar to re_delta_redshift~0.5. +re_delta_redshift = 0.5 +#re_ionization_frac=-1 sets to become fully ionized using YHe to get helium contribution +#Otherwise x_e varies from 0 to re_ionization_frac +re_ionization_frac = -1 + + +#RECFAST 1.5 recombination parameters; +RECFAST_fudge = 1.14 +RECFAST_fudge_He = 0.86 +RECFAST_Heswitch = 6 +RECFAST_Hswitch = T + +#Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, +# neutrino density iso =4, neutrino velocity iso = 5, axion isocurvature = 6) +initial_condition = 1 +#If above is zero, use modes in the following (totally correlated) proportions +#Note: we assume all modes have the same initial power spectrum +initial_vector = 1 0 0 0 0 1 + +#For vector modes: 0 for regular (neutrino vorticity mode), 1 for magnetic +vector_mode = 0 + +#Normalization +COBE_normalize = F +##CMB_outputscale scales the output Cls +#To get MuK^2 set realistic initial amplitude (e.g. scalar_amp(1) = 2.3e-9 above) and +#otherwise for dimensionless transfer functions set scalar_amp(1)=1 and use +CMB_outputscale = 1 +#CMB_outputscale = 7.42835025e12 + +#Transfer function settings, transfer_kmax=0.5 is enough for sigma_8 +#transfer_k_per_logint=0 sets sensible non-even sampling; +#transfer_k_per_logint=5 samples fixed spacing in log-k +#transfer_interp_matterpower =T produces matter power in regular interpolated grid in log k; +# use transfer_interp_matterpower =F to output calculated values (e.g. for later interpolation) +transfer_high_precision = T +transfer_kmax = 1000 +transfer_k_per_logint = 5 +transfer_num_redshifts = 1 +transfer_interp_matterpower = T +transfer_redshift(1) = 100 +transfer_filename(1) = transfer_out.dat +#Matter power spectrum output against k/h in units of h^{-3} Mpc^3 +transfer_matterpower(1) = matterpower.dat + +#Output files not produced if blank. make camb_fits to use use the FITS setting. +scalar_output_file = scalCls.dat +vector_output_file = vecCls.dat +tensor_output_file = tensCls.dat +total_output_file = totCls.dat +lensed_output_file = lensedCls.dat +lensed_total_output_file =lensedtotCls.dat +lens_potential_output_file = lenspotentialCls.dat +FITS_filename = scalCls.fits + +#Bispectrum parameters if required; primordial is currently only local model (fnl=1) +#lensing is fairly quick, primordial takes several minutes on quad core +do_lensing_bispectrum = F +do_primordial_bispectrum = F + +#1 for just temperature, 2 with E +bispectrum_nfields = 1 +#set slice non-zero to output slice b_{bispectrum_slice_base_L L L+delta} +bispectrum_slice_base_L = 0 +bispectrum_ndelta=3 +bispectrum_delta(1)=0 +bispectrum_delta(2)=2 +bispectrum_delta(3)=4 +#bispectrum_do_fisher estimates errors and correlations between bispectra +#note you need to compile with LAPACK and FISHER defined to use get the Fisher info +bispectrum_do_fisher= F +#Noise is in muK^2, e.g. 2e-4 roughly for Planck temperature +bispectrum_fisher_noise=0 +bispectrum_fisher_noise_pol=0 +bispectrum_fisher_fwhm_arcmin=7 +#Filename if you want to write full reduced bispectrum (at sampled values of l_1) +bispectrum_full_output_file= +bispectrum_full_output_sparse=F + +##Optional parameters to control the computation speed,accuracy and feedback + +#If feedback_level > 0 print out useful information computed about the model +feedback_level = 1 + +# 1: curved correlation function, 2: flat correlation function, 3: inaccurate harmonic method +lensing_method = 1 +accurate_BB = F + + +#massive_nu_approx: 0 - integrate distribution function +# 1 - switch to series in velocity weight once non-relativistic +massive_nu_approx = 1 + +#Whether you are bothered about polarization. +accurate_polarization = T + +#Whether you are bothered about percent accuracy on EE from reionization +accurate_reionization = T + +#whether or not to include neutrinos in the tensor evolution equations +do_tensor_neutrinos = T + +#Whether to turn off small-scale late time radiation hierarchies (save time,v. accurate) +do_late_rad_truncation = T + +#Computation parameters +#if number_of_threads=0 assigned automatically +number_of_threads = 0 + +#Default scalar accuracy is about 0.3% (except lensed BB) if high_accuracy_default=F +#If high_accuracy_default=T the default taget accuracy is 0.1% at L>600 (with boost parameter=1 below) +#Try accuracy_boost=2, l_accuracy_boost=2 if you want to check stability/even higher accuracy +#Note increasing accuracy_boost parameters is very inefficient if you want higher accuracy, +#but high_accuracy_default is efficient + +high_accuracy_default=T + +#Increase accuracy_boost to decrease time steps, use more k values, etc. +#Decrease to speed up at cost of worse accuracy. Suggest 0.8 to 3. +accuracy_boost = 2.1 + +#Larger to keep more terms in the hierarchy evolution. +l_accuracy_boost = 2.1 + +#Increase to use more C_l values for interpolation. +#Increasing a bit will improve the polarization accuracy at l up to 200 - +#interpolation errors may be up to 3% +#Decrease to speed up non-flat models a bit +l_sample_boost = 128 diff --git a/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/set_velocities_zero.py b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/set_velocities_zero.py new file mode 100755 index 0000000000..21b9943c50 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/set_velocities_zero.py @@ -0,0 +1,58 @@ +################################################################################################## +# This script is used for setting the velocity-dependent transfer fucntion (TF) to 1e-8 for TF # +# file generated by CAMB(planck_2018_transfer_out.dat) and axsionCAMB(planck_2018_axion_transfer # +# _out.dat). Be sure to prepare the input transfer function files and specify their names!! The # +# output transfer function file will be produced accordingly (corresponding to CAMB/axionCAMB # +# input). # +################################################################################################## + +import sys +import numpy as np +import pandas as pd + +if len(sys.argv) != 4: + raise RuntimeError("Input Error: Please select whether-axion-physics-included (T/F), file name of the input transfer function, and file name of the modified output transfer function!") +else: + if sys.argv[1] not in ["T","F"]: + raise RuntimeError("Input Error: Please select T/F if the input transfer function includes/excludes axion-physics!") + input_file = sys.argv[2] + output_file = sys.argv[3] + +if sys.argv[1] == "F": + # no axion + df_trans = pd.read_csv(input_file, delimiter='\s+', index_col=False) + col = list(df_trans.columns[1:]) + col.append("None") + df_trans.columns = col + df_trans = df_trans.iloc[:,:-1] + df_trans['v_CDM'] = 1e-8*np.ones(df_trans['v_CDM'].shape) + df_trans['v_b'] = 1e-8*np.ones(df_trans['v_b'].shape) + df_trans['v_b-v_c'] = 1e-8*np.ones(df_trans['v_b-v_c'].shape) + + df_trans.to_csv(output_file, header=None, sep='\t', float_format = "%.6e", index=None) +#################################################################################################################### +elif sys.argv[1] == "T": + # with axion + df_trans_axion = pd.read_csv(input_file, delimiter='\s+', index_col=False, header=None) + + # switch the CDM column with axion column + temp = df_trans_axion.iloc[:,1].copy() + df_trans_axion.iloc[:,1] = df_trans_axion.iloc[:,6].copy() + df_trans_axion.iloc[:,6] = temp.copy() + + # switch the CDM column with total matter column + temp = df_trans_axion.iloc[:,6].copy() + df_trans_axion.iloc[:,6] = df_trans_axion.iloc[:,-1].copy() + df_trans_axion.iloc[:,-1] = temp.copy() + + # subtitute the data with negative TF by its absolute value + df_trans_axion = abs(df_trans_axion) + + # subtitute the data with negative TF by 1e-8 + for i in range(9,13): + df_trans_axion[i] = 1e-8*np.ones(len(df_trans_axion.iloc[:,0])) + + df_trans_axion.to_csv(output_file, header=None, sep='\t', float_format = "%.6e", index=None) +#################################################################################################################### +else: + raise RuntimeError("Input Error: Please select do axion or not (T/F)!")