BEoRN coupled to 21cmFast¶
[1]:
import numpy as np
import matplotlib.pyplot as plt
[2]:
import beorn
import tools21cm as t2c
Define parameters¶
beorn contains a parameter object beorn.par() conatining the default parameter values. This object is essential for most of the beorn simulation functions/classes.
[3]:
param = beorn.par()
Cosmology¶
[4]:
# cosmology
param.cosmo.Om = 0.315
param.cosmo.Ob = 0.045
param.cosmo.Ol = 0.685
param.cosmo.h = 0.673
param.cosmo.ps = '../../examples/PCDM_Planck.dat'
param.cosmo.corr_fct = '../../examples/corr_fct.dat'
Source parameters¶
[5]:
# lyman-alpha photons
param.source.N_al = 9690 #1500
param.source.alS_lyal = 0.0
# ionising photons
param.source.Nion = 3000 #5000
# X-ray photons
param.source.E_min_xray = 500
param.source.E_max_xray = 10000
param.source.E_min_sed_xray = 200
param.source.E_max_sed_xray = 10000
param.source.alS_xray = 1.5
param.source.cX = 3.4e40
# photon escape fraction
param.source.f0_esc = 0.2
param.source.pl_esc = 0.5
# stellar-to-halo
param.source.f_st = 0.14
param.source.g1 = 0.49
param.source.g2 = -0.61
param.source.g3 = 4
param.source.g4 = -4
param.source.Mp = 1.6e11 * param.cosmo.h
param.source.Mt = 1e9
# Minimum star forming halo
param.source.M_min= 1e8
Simulation parameters¶
[6]:
# sim
param.sim.M_i_min = 3e1 * np.exp(0.79*(25-40))
param.sim.M_i_max = 8e7 * np.exp(0.79*(25-40))
param.sim.model_name = 'py21cmfast_test'
param.sim.cores = 1 # nbr of cores to use
param.sim.binn = 40 # nbr of halo mass bin
# Box size and Number of pixels
param.sim.Lbox = 100
param.sim.Ncell = 128
param.sim.halo_catalogs = None
param.sim.thresh_pixel = 20*(param.sim.Ncell/128)**3
param.sim.dens_fields = None
param.sim.dens_field_type = '21cmfast'
param.sim.data_dir = '../../examples/py21cmfast_test'
param.sim.store_grids = True #'replace' # True # False
# bubble expansion solver
param.solver.z_ini = 40
param.solver.z_end = 6
param.solver.Nz = 100
Simulate matter evolution¶
beorn can either use density snapshots and halo catalogues from any N-body simulation or run a simulation code that it is coupled to. Here an example of running 21cmFast to model matter evolution is shown.
[7]:
redshift_list = [9,8,7]
matter_field = beorn.simulate_matter_21cmfast(param, redshift_list=redshift_list, IC=None)
You have not specified a data_dir to store the 21cmFAST data. By default it is ./21cmFAST_data.
Simulating matter evolution with 21cmFast...
random seed: 12345
We strongly advice to specify a path for param.sim.dens_field to write the densities.
It should be of the form path+dir+name. For instance param.sim.dens_field = ./dir_density/dens_21cmFast_z
We strongly advice to specify a path for param.sim.halo_catalogs to write the halo catalogues.
It should be of the form path+dir+name. For instance param.sim.dens_field = ./dir_density/dens_21cmFast_z
You have not specified a data_dir to store the 21cmFAST data. By default it is ./21cmFAST_data.
The following parameters to FlagOptions are not supported: ['USE_VELS_AUX']
The following parameters to FlagOptions are not supported: ['USE_VELS_AUX']
The following parameters to FlagOptions are not supported: ['USE_VELS_AUX']
param.sim.Lbox is in Mpc/h. Halo catalogs catalogs have masses in Msol/h and positions in Mpc/h.
param.sim.Lbox is in Mpc/h. Halo catalogs catalogs have masses in Msol/h and positions in Mpc/h.
The following parameters to FlagOptions are not supported: ['USE_VELS_AUX']
The following parameters to FlagOptions are not supported: ['USE_VELS_AUX']
param.sim.Lbox is in Mpc/h. Halo catalogs catalogs have masses in Msol/h and positions in Mpc/h.
...done | Runtime = 00:00:48
[8]:
print('No. of density and halo fields =', len(matter_field['dens']), len(matter_field['halo_list']))
No. of density and halo fields = 3 3
Simulating reionization and cosmic dawn¶
The results from 21cmFast is defined in the parameter object.
[9]:
param.sim.dens_fields = matter_field['dens']
param.sim.halo_catalogs = matter_field['halo_list']
param.solver.Nz = 100
[10]:
beorn.run.compute_profiles(param)
Computing Temperature (Tk), Lyman-α and ionisation fraction (xHII) profiles...
param.solver.Nz is given as an integer. We define z values in linspace from 40 to 6
param.solver.fXh is set to constant. We will assume f_X,h = 2e-4**0.225
... Profiles stored in dir ./profiles.
It took 00:00:40 to compute the profiles.
[12]:
param.solver.Nz = [ke for ke in matter_field['dens'].keys()]
beorn.run.paint_boxes(param,
RSD=False,
ion=True,
temp=True,
dTb=True,
lyal=True,
check_exists=False,
cross_corr=False,
)
Painting profiles on a grid with 128 pixels per dim. Box size is 100 cMpc/h.
param.solver.Nz is given as a list.
Core nbr 0 is taking care of z = 9
----- Painting 3D map for z = 9 -------
param.sim.dens_field_type should be either 21cmFAST or pkdgrav.
There are 74862 halos at z= 9
Looping over halo mass bins and painting profiles on 3D grid ....
Quick calculation from the profiles predicts xHII = 0.0721
30616 halos in mass bin 23 . It took 00:00:01 to paint the profiles.
23307 halos in mass bin 24 . It took 00:00:01 to paint the profiles.
14657 halos in mass bin 25 . It took 00:00:02 to paint the profiles.
2445 halos in mass bin 26 . It took 00:00:03 to paint the profiles.
2586 halos in mass bin 27 . It took 00:00:04 to paint the profiles.
887 halos in mass bin 28 . It took 00:00:04 to paint the profiles.
278 halos in mass bin 29 . It took 00:00:05 to paint the profiles.
72 halos in mass bin 30 . It took 00:00:06 to paint the profiles.
14 halos in mass bin 31 . It took 00:00:07 to paint the profiles.
.... Done painting profiles.
Dealing with the overlap of ionised bubbles....
initial sum of ionized fraction : 146522.754
Universe not fully ionized : xHII is 0.0699
7702 connected regions.
there are 7221 connected regions with less than 20.0 pixels. They contain a fraction 0.109 of the total ionisation fraction.
final xion sum: 146522.754
.... Done. It took: 00:00:04 to redistribute excess photons from the overlapping regions.
--- Including Salpha fluctuations in dTb ---
--- Including xcoll fluctuations in dTb ---
----- Snapshot at z = 9 is done -------
Core nbr 0 is taking care of z = 8
----- Painting 3D map for z = 8 -------
param.sim.dens_field_type should be either 21cmFAST or pkdgrav.
There are 151215 halos at z= 8
Looping over halo mass bins and painting profiles on 3D grid ....
Quick calculation from the profiles predicts xHII = 0.2056
61138 halos in mass bin 22 . It took 00:00:01 to paint the profiles.
46064 halos in mass bin 23 . It took 00:00:01 to paint the profiles.
30930 halos in mass bin 24 . It took 00:00:02 to paint the profiles.
8216 halos in mass bin 25 . It took 00:00:03 to paint the profiles.
3141 halos in mass bin 26 . It took 00:00:04 to paint the profiles.
757 halos in mass bin 27 . It took 00:00:05 to paint the profiles.
693 halos in mass bin 28 . It took 00:00:05 to paint the profiles.
217 halos in mass bin 29 . It took 00:00:06 to paint the profiles.
53 halos in mass bin 30 . It took 00:00:07 to paint the profiles.
6 halos in mass bin 31 . It took 00:00:08 to paint the profiles.
.... Done painting profiles.
Dealing with the overlap of ionised bubbles....
initial sum of ionized fraction : 405662.182
Universe not fully ionized : xHII is 0.1934
9415 connected regions.
there are 8758 connected regions with less than 20.0 pixels. They contain a fraction 0.0707 of the total ionisation fraction.
final xion sum: 405662.182
.... Done. It took: 00:00:07 to redistribute excess photons from the overlapping regions.
--- Including Salpha fluctuations in dTb ---
--- Including xcoll fluctuations in dTb ---
----- Snapshot at z = 8 is done -------
Core nbr 0 is taking care of z = 7
----- Painting 3D map for z = 7 -------
param.sim.dens_field_type should be either 21cmFAST or pkdgrav.
There are 330146 halos at z= 7
Looping over halo mass bins and painting profiles on 3D grid ....
Quick calculation from the profiles predicts xHII = 0.4822
114201 halos in mass bin 20 . It took 00:00:01 to paint the profiles.
52168 halos in mass bin 21 . It took 00:00:01 to paint the profiles.
111683 halos in mass bin 22 . It took 00:00:02 to paint the profiles.
27258 halos in mass bin 23 . It took 00:00:03 to paint the profiles.
15047 halos in mass bin 24 . It took 00:00:03 to paint the profiles.
5903 halos in mass bin 25 . It took 00:00:04 to paint the profiles.
2586 halos in mass bin 26 . It took 00:00:05 to paint the profiles.
884 halos in mass bin 27 . It took 00:00:06 to paint the profiles.
214 halos in mass bin 28 . It took 00:00:07 to paint the profiles.
152 halos in mass bin 29 . It took 00:00:07 to paint the profiles.
45 halos in mass bin 30 . It took 00:00:08 to paint the profiles.
5 halos in mass bin 31 . It took 00:00:09 to paint the profiles.
.... Done painting profiles.
Dealing with the overlap of ionised bubbles....
initial sum of ionized fraction : 973517.009
Universe not fully ionized : xHII is 0.4642
8666 connected regions.
there are 8178 connected regions with less than 20.0 pixels. They contain a fraction 0.0303 of the total ionisation fraction.
final xion sum: 973517.009
.... Done. It took: 00:00:05 to redistribute excess photons from the overlapping regions.
--- Including Salpha fluctuations in dTb ---
--- Including xcoll fluctuations in dTb ---
----- Snapshot at z = 7 is done -------
Finished painting the maps. They are stored in ./grid_output. It took in total: 00:00:43 to paint the grids.
[13]:
zz = 7
xHII_grid = beorn.load_grid(param,z=zz,type='bubbles')
beorn.plot_2d_map(xHII_grid, 100, 64, qty='xHII')
dTb_grid = beorn.load_grid(param,z=zz,type='dTb')
beorn.plot_2d_map(dTb_grid, 100, 64, scale='lin', qty='dTb [mK]')
Ncell is 128
Ncell is 128
