From d26a599cf4e2933b24390ef091336e52b0b69451 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Tue, 7 May 2024 16:49:41 +0800 Subject: [PATCH 01/11] add CAMB/axionCAMB initial condition tool for LSS simulation --- .../README.txt | 13 + .../ics_example.conf | 79 +++++ .../make_umic_from_hdf5.py | 167 ++++++++++ .../planck_2018.ini | 284 ++++++++++++++++++ .../planck_2018_axion.ini | 243 +++++++++++++++ .../set_velocities_zero.py | 56 ++++ 6 files changed, 842 insertions(+) create mode 100644 tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/README.txt create mode 100644 tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/ics_example.conf create mode 100755 tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/make_umic_from_hdf5.py create mode 100644 tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018.ini create mode 100644 tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018_axion.ini create mode 100755 tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/set_velocities_zero.py 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..ab509dd78b --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/README.txt @@ -0,0 +1,13 @@ +The procedures of this LSS simulation includes: + +(a) Using CAMB/axionCAMB to produce transfer function (the input file planck_2018.ini/planck_2018_axion.ini for CAMB/axionCAMB is attached for reference) +# Note + 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 \simga_8 = 0.8119112 with As = 2.100549e-9 , to new value \simga_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 \simga_8 = 0.818. This is due to power spectrum is proportional to As*f(k) . + +(b) Run the script "./set_velocities_zero.py F(T)" to produce the CAMB(axionCAMB) transfer function: planck_2018_transfer_out_no_vel.dat(planck_2018_axion_transfer_out_no_vel.dat) needed by MUSIC. + +(c) Modifying the cosmology paramters consistent with CAMB/axionCAMB .ini files, as well as the input file parameters in ics_example.conf (transfer = camb_file; transfer_file = planck_2018_axion_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). Then run the MUSIC: "./MUSIC ics_example.conf". + +(d) Run the script "./make_umic_from_hdf5.py S(D) hdf5_filename" to produce UM_IC file for single/double precision GAMER. + +(e) Run the GAMER. 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..f1b0ce5efc --- /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_axion_transfer_out_no_vel.dat +#transfer_file = planck_2018_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..4b13144b4a --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/make_umic_from_hdf5.py @@ -0,0 +1,167 @@ +#!/usr/bin/env python3.7 + +############################################################################################# +# 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 phidm_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 quation in spcetrum way: # +# a*phidm_mass*(div(velocity))/h_bar = del(phase) # +# div is evaluated via Richardson extrapolation with an adjustable order # +# (4) Convert density and phase field to real and imaginary part of wave function # +# (5) Output as binary file "UM_IC" # +# # +# Unit(s): MUCIS v.s. here # +# (1) density: back ground density ; back ground density # +# (2) velocity: box length*H_0 ; m/s # +############################################################################################# + +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 +phidm_mass = 8e-23 #eV +h_bar i = 1.0545718e-34 +eV_to_kg = 1.78266192162790e-36 +Mpc_to_m = 3.0856776e22 +##################################################################################################################### + +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*1000. + 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() + +# Bug test################################### + #growing_factor = 5./3. + #density_hdf5 = (growing_factor*density_hdf5+1.) + #criteria = (density_hdf5<0.) + #print("Percentage of over-density smaller than 0.: %.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 +############################################# + +# Normal version############################# + 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*phidm_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 + phi_fft = np.fft.irfftn(v_div_k) + # Rescale to correct unit + phi_fft *= box_length/N**2 + phi_fft -= phi_fft.min() + + # Convert to wave function + wf_real = density_hdf5*np.cos(phi_fft) + wf_imag = density_hdf5*np.sin(phi_fft) + + 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..bbedb5b7b4 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/planck_2018.ini @@ -0,0 +1,284 @@ +#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 +#transfer_power_var = 2 + +#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..cd4ea41ae4 --- /dev/null +++ b/tool/inits/GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/set_velocities_zero.py @@ -0,0 +1,56 @@ +#!/usr/bin/env python3.7 + +################################################################################################## +# 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 name!! The # +# output file planck_2018_transfer_out_no_vel.dat/planck_2018_axion_transfer_out.dat will be # +# produced accordingly(corresponding to CAMB/axionCAMB input). # +################################################################################################## + +import sys +import numpy as np +import pandas as pd + +if len(sys.argv) != 2: + raise RuntimeError("Input Error: Please select do axion or not (T/F)!") + +if sys.argv[1] == "F": + # no axion + df_trans = pd.read_csv('planck_2018_transfer_out.dat', 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('planck_2018_transfer_out_no_vel.dat',header=None, sep='\t', float_format = "%.6e", index=None) +#################################################################################################################### +elif sys.argv[1] == "T": + # with axion + df_trans_axion = pd.read_csv('planck_2018_axion_transfer_out.dat', 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.iloc[:,:100].to_csv('planck_2018_axion_transfer_out_no_vel.dat', \ + header=None, sep='\t', float_format = "%.6e", index=None) +#################################################################################################################### +else: + raise RuntimeError("Input Error: Please select do axion or not (T/F)!") From 31ad83681fb85fe1547a0ac6d4a316cb1ed8f973 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Tue, 7 May 2024 18:22:11 +0800 Subject: [PATCH 02/11] minor --- .../make_umic_from_hdf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index 4b13144b4a..d259679e2c 100755 --- 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 @@ -40,7 +40,7 @@ def GRAD(field, axis, h, order): return grad[-1]/h # Physical Constant phidm_mass = 8e-23 #eV -h_bar i = 1.0545718e-34 +h_bar = 1.0545718e-34 eV_to_kg = 1.78266192162790e-36 Mpc_to_m = 3.0856776e22 ##################################################################################################################### From 3081181b6816a3a40f0a44aee4040ce5eb6c1184 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Tue, 9 Jul 2024 16:37:06 +0800 Subject: [PATCH 03/11] add explanations in the header for the physical meaning of MUSIC HDF5 velocity --- .../make_umic_from_hdf5.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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 index d259679e2c..d27cee6e87 100755 --- 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 @@ -19,9 +19,12 @@ # (4) Convert density and phase field to real and imaginary part of wave function # # (5) Output as binary file "UM_IC" # # # -# Unit(s): MUCIS v.s. here # +# Unit(s): MUSIC v.s. here # # (1) density: back ground density ; back ground density # # (2) velocity: box length*H_0 ; m/s # +# # +# Phyiscal meaning for MUSIC's velocity recorded in HDF5 is perculiar velocity (in physical # +# frame, not comoving frame) # ############################################################################################# import os, h5py, subprocess, re, sys From 35a916151fde778a72c8cb45acb8f0318d4ef747 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Tue, 9 Jul 2024 16:58:40 +0800 Subject: [PATCH 04/11] add formula to the explanation for physical meaning of MUSIC HDF5 velocity --- .../make_umic_from_hdf5.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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 index d27cee6e87..969a6e1764 100755 --- 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 @@ -24,7 +24,9 @@ # (2) velocity: box length*H_0 ; m/s # # # # Phyiscal meaning for MUSIC's velocity recorded in HDF5 is perculiar velocity (in physical # -# frame, not comoving frame) # +# frame, not comoving frame), i.e., after multiplying box_length*1000, 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 From bcb52ed85fed142648f9b90873110e27b6be49b0 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Tue, 9 Jul 2024 17:02:36 +0800 Subject: [PATCH 05/11] fix typo --- .../make_umic_from_hdf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index 969a6e1764..d880eb4913 100755 --- 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 @@ -25,7 +25,7 @@ # # # Phyiscal meaning for MUSIC's velocity recorded in HDF5 is perculiar velocity (in physical # # frame, not comoving frame), i.e., after multiplying box_length*1000, we have: # -# (v_hdf5),i = (v_peculiar,physical),i = a*(dx_{comoving}),i/dt # +# (v_hdf5),i = (v_perculiar,physical),i = a*(dx_{comoving}),i/dt # # where i=x, y, z # ############################################################################################# From 70b0cebe81f69daf43e9471a3876d28240752a01 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Tue, 9 Jul 2024 17:03:29 +0800 Subject: [PATCH 06/11] fix typo --- .../make_umic_from_hdf5.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 index d880eb4913..56a63eeee5 100755 --- 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 @@ -23,9 +23,9 @@ # (1) density: back ground density ; back ground density # # (2) velocity: box length*H_0 ; m/s # # # -# Phyiscal meaning for MUSIC's velocity recorded in HDF5 is perculiar velocity (in physical # +# Physical meaning for MUSIC's velocity recorded in HDF5 is peculiar velocity (in physical # # frame, not comoving frame), i.e., after multiplying box_length*1000, we have: # -# (v_hdf5),i = (v_perculiar,physical),i = a*(dx_{comoving}),i/dt # +# (v_hdf5),i = (v_peculiar,physical),i = a*(dx_{comoving}),i/dt # # where i=x, y, z # ############################################################################################# From 973f2a7ce3fd2483fa2e6210c87c8bb6dec5436f Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Wed, 15 Jan 2025 21:00:02 +0800 Subject: [PATCH 07/11] modify according to PR review --- .../README.txt | 80 ++++++++++++-- .../ics_example.conf | 100 +++++++++--------- .../make_umic_from_hdf5.py | 52 ++++----- .../set_velocities_zero.py | 24 +++-- 4 files changed, 155 insertions(+), 101 deletions(-) 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 index ab509dd78b..2371f8013f 100644 --- 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 @@ -1,13 +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 (CAMB)] + [python3 ./set_velocities_zero.py T planck_2018_axion_trnasfer_out.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) Using CAMB/axionCAMB to produce transfer function (the input file planck_2018.ini/planck_2018_axion.ini for CAMB/axionCAMB is attached for reference) -# Note - 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 \simga_8 = 0.8119112 with As = 2.100549e-9 , to new value \simga_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 \simga_8 = 0.818. This is due to power spectrum is proportional to As*f(k) . +--------------------------------------------------------------------------------------------------------------- -(b) Run the script "./set_velocities_zero.py F(T)" to produce the CAMB(axionCAMB) transfer function: planck_2018_transfer_out_no_vel.dat(planck_2018_axion_transfer_out_no_vel.dat) needed by MUSIC. +(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 -(c) Modifying the cosmology paramters consistent with CAMB/axionCAMB .ini files, as well as the input file parameters in ics_example.conf (transfer = camb_file; transfer_file = planck_2018_axion_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). Then run the MUSIC: "./MUSIC ics_example.conf". +(b) Run the script: "python3 ./set_velocities_zero.py F(T) planck_2018_trnasfer_out.dat(planck_2018_axion_trnasfer_out.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 -(d) Run the script "./make_umic_from_hdf5.py S(D) hdf5_filename" to produce UM_IC file for single/double precision GAMER. +(c) Modifying the cosmology paramters 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 \simga_8 = 0.8119112 with As = 2.100549e-9 , to new value \simga_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 \simga_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 index f1b0ce5efc..13e920951f 100644 --- 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 @@ -1,79 +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 +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_axion_transfer_out_no_vel.dat +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 +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 +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 +#format = enzo +#filename = ic.enzo ##Gadget-2 (type=1: high-res particles, type=5: rest) -#format = gadget2 -#filename = ics_gadget.dat +#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 +#format = grafic2 +#filename = ics_ramses +#ramses_nml = yes ##TIPSY compatible with PKDgrav and Gasoline -#format = tipsy -#filename = ics_tipsy.dat +#format = tipsy +#filename = ics_tipsy.dat ## NYX compatible output format ##requires boxlib installation and boxlib enabled in Makefile -#format = nyx -#filename = init +#format = nyx +#filename = init [poisson] -fft_fine = yes -accuracy = 1e-5 -pre_smooth = 3 -post_smooth = 3 -smoother = gs +fft_fine = yes +accuracy = 1e-5 +pre_smooth = 3 +post_smooth = 3 +smoother = gs laplace_order = 6 -grad_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 index 56a63eeee5..f661d99d17 100755 --- 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 @@ -1,11 +1,9 @@ -#!/usr/bin/env python3.7 - ############################################################################################# # 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 phidm_mass in physical constant if needed. # +# (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 # @@ -13,9 +11,10 @@ # 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 quation in spcetrum way: # -# a*phidm_mass*(div(velocity))/h_bar = del(phase) # +# (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" # # # @@ -24,7 +23,7 @@ # (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*1000, we have: # +# 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 # ############################################################################################# @@ -44,10 +43,11 @@ def GRAD(field, axis, h, order): grad[o:] = (4.**o*grad[o:]-grad[o-1:-1]) / (4.**o-1.) return grad[-1]/h # Physical Constant -phidm_mass = 8e-23 #eV +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: @@ -94,7 +94,7 @@ def GRAD(field, axis, h, order): z = eval(re.sub("\n","", z)) print("z_start is %.2f ."%z) - factor = box_length*1000. + factor = box_length*km_to_m. N = 2**level h = 1./N k_factor = 2.*np.pi/N @@ -103,27 +103,12 @@ def GRAD(field, axis, h, order): 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] + 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() -# Bug test################################### - #growing_factor = 5./3. - #density_hdf5 = (growing_factor*density_hdf5+1.) - #criteria = (density_hdf5<0.) - #print("Percentage of over-density smaller than 0.: %.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 -############################################# - -# Normal version############################# + criteria = (density_hdf5<-1.) print("Percentage of over-density smaller than -1: %.8f %%."%(100*criteria.sum()/2**(3*level)) ) density_hdf5[criteria] = -1. @@ -135,14 +120,13 @@ def GRAD(field, axis, h, order): 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*phidm_mass*eV_to_kg/h_bar + v_div *= a*psidm_mass*eV_to_kg/h_bar # Do forward DFT v_div_k = np.fft.rfftn(v_div) @@ -152,14 +136,14 @@ def GRAD(field, axis, h, order): 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 - phi_fft = np.fft.irfftn(v_div_k) + phase = np.fft.irfftn(v_div_k) # Rescale to correct unit - phi_fft *= box_length/N**2 - phi_fft -= phi_fft.min() + phase *= box_length/N**2 + phase -= phase.min() # Convert to wave function - wf_real = density_hdf5*np.cos(phi_fft) - wf_imag = density_hdf5*np.sin(phi_fft) + 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) 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 index cd4ea41ae4..052433e909 100755 --- 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 @@ -1,23 +1,26 @@ -#!/usr/bin/env python3.7 - ################################################################################################## # 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 name!! The # -# output file planck_2018_transfer_out_no_vel.dat/planck_2018_axion_transfer_out.dat will be # -# produced accordingly(corresponding to CAMB/axionCAMB input). # +# 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) != 2: - raise RuntimeError("Input Error: Please select do axion or not (T/F)!") +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('planck_2018_transfer_out.dat', delimiter='\s+', index_col=False) + 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 @@ -26,11 +29,11 @@ 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('planck_2018_transfer_out_no_vel.dat',header=None, sep='\t', float_format = "%.6e", index=None) + 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('planck_2018_axion_transfer_out.dat', delimiter='\s+', index_col=False, header=None) + 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() @@ -49,8 +52,7 @@ for i in range(9,13): df_trans_axion[i] = 1e-8*np.ones(len(df_trans_axion.iloc[:,0])) - df_trans_axion.iloc[:,:100].to_csv('planck_2018_axion_transfer_out_no_vel.dat', \ - header=None, sep='\t', float_format = "%.6e", index=None) + 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)!") From 4a63f1bb6c899b474a893bafd750d25a3f8d774d Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Wed, 15 Jan 2025 21:27:46 +0800 Subject: [PATCH 08/11] align equal sign --- .../planck_2018.ini | 151 +++++++++--------- 1 file changed, 75 insertions(+), 76 deletions(-) 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 index bbedb5b7b4..fd7efef592 100644 --- 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 @@ -6,21 +6,21 @@ #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 +output_root = planck_2018 #What to do -get_scalar_cls = T -get_vector_cls = F -get_tensor_cls = F -get_transfer = T +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 +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 +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 @@ -37,44 +37,44 @@ 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 +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 +dark_energy_model = ppf #effective equation of state parameter for dark energy -w = -1 +w = -1 #constant comoving sound speed of the dark energy (1=quintessence), for fluid model -cs2_lam = 1 +cs2_lam = 1 #For PPF model can set these instead: -#wa = 0 +#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 +#use_tabulated_w = T +#wafile = wa.dat -temp_cmb = 2.7255 -helium_fraction = 0.2454006 +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 +massless_neutrinos = 2.046 #number of distinct mass eigenstates -nu_mass_eigenstates = 1 +nu_mass_eigenstates = 1 #array of the integer number of physical neutrinos per eigenstate, e.g. massive_neutrinos = 2 1 -massive_neutrinos = 1 +massive_neutrinos = 1 #specify whether all neutrinos should have the same temperature, specified from fractional part of massless_neutrinos -share_delta_neff = T +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 +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 @@ -99,7 +99,7 @@ tensor_parameterization = 1 #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 +#tensor_amp = 4e-10 #note vector modes use the scalar settings above @@ -118,13 +118,13 @@ re_delta_redshift = 0.5 re_ionization_frac = -1 #Parameters for second reionization of helium -re_helium_redshift = 3.5 +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 = 1.14 RECFAST_fudge_He = 0.86 RECFAST_Heswitch = 6 RECFAST_Hswitch = T @@ -144,112 +144,111 @@ RECFAST_Hswitch = T #Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, # neutrino density iso =4, neutrino velocity iso = 5) -initial_condition = 1 +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 +initial_vector = -1 0 0 0 0 #For vector modes: 0 for regular (neutrino vorticity mode), 1 for magnetic -vector_mode = 0 +vector_mode = 0 #Normalization -COBE_normalize = F +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 +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_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 +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 +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 -#transfer_power_var = 2 +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 +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 +do_lensing_bispectrum = F +do_primordial_bispectrum = F #1 for just temperature, 2 with E -bispectrum_nfields = 1 +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_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 +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 +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 +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 +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 +feedback_level = 1 #whether to start output files with comment describing columns output_file_headers = T #write out various derived parameters -derived_parameters = T +derived_parameters = T # 1: curved correlation function, 2: flat correlation function, 3: inaccurate harmonic method -lensing_method = 1 -accurate_BB = F +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 +massive_nu_approx = 1 #Whether you are bothered about polarization. -accurate_polarization = T +accurate_polarization = T #Whether you are bothered about percent accuracy on EE from reionization -accurate_reionization = T +accurate_reionization = T #whether or not to include neutrinos in the tensor evolution equations -do_tensor_neutrinos = T +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 +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 @@ -260,7 +259,7 @@ do_late_rad_truncation = T #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 +halofit_version = 5 #Computation parameters #if number_of_threads=0 assigned automatically From 02e35640a110758ec6a455024623dec7fc09ecf2 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Mon, 3 Feb 2025 14:42:49 +0800 Subject: [PATCH 09/11] Minor: fix typo & add documnet for new parameters --- .../GEN_UM_IC_WITH_PHASE_WITH_TRANSFER_FUNC/README.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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 index 2371f8013f..5ee98ee450 100644 --- 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 @@ -13,8 +13,8 @@ The procedures of this LSS simulation includes: | | | | -------> [(b) set_velocities_zero.py] <----- - [python3 ./set_velocities_zero.py F planck_2018_trnasfer_out.dat (CAMB)] - [python3 ./set_velocities_zero.py T planck_2018_axion_trnasfer_out.dat (axionCAMB)] + [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) @@ -44,14 +44,14 @@ The procedures of this LSS simulation includes: 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)". +(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 consistent with CAMB/axionCAMB .ini files, as well as the input file parameters in MUSIC input file ics_example.conf: +(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 @@ -65,7 +65,7 @@ The procedures of this LSS simulation includes: # 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 \simga_8 = 0.8119112 with As = 2.100549e-9 , to new value \simga_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 \simga_8 = 0.818. + 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 From a20c06c89d0b60a6c13c46852bfe8f5aef3a6156 Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Mon, 3 Feb 2025 21:33:08 +0800 Subject: [PATCH 10/11] Minor: bug fix --- .../make_umic_from_hdf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 index f661d99d17..8c68adf752 100755 --- 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 @@ -94,7 +94,7 @@ def GRAD(field, axis, h, order): z = eval(re.sub("\n","", z)) print("z_start is %.2f ."%z) - factor = box_length*km_to_m. + factor = box_length*km_to_m N = 2**level h = 1./N k_factor = 2.*np.pi/N From 189742dabbc01822e7650649429673bb411885fd Mon Sep 17 00:00:00 2001 From: koarakawaii Date: Wed, 24 Sep 2025 15:59:20 +0800 Subject: [PATCH 11/11] Minor: fix grammar error --- .../set_velocities_zero.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 index 052433e909..21b9943c50 100755 --- 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 @@ -1,8 +1,8 @@ ################################################################################################## # 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 name!! The # -# output transfer function file will be produced accordingly(corresponding to CAMB/axionCAMB # +# _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). # ##################################################################################################