LSS Simulation with CAMB/axionCAMB Transfer Function

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
+

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

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()