pyC\(^2\)Ray Simulation Class¶
This tutorial presents the Python class that manages the pyC\(^2\)Ray simulations for your custom EoR simulation.
[2]:
import pyc2ray as pc2r
import numpy as np, yaml
import matplotlib.pyplot as plt
import astropy.units as u
import astropy.constants as cst
The C2Ray Python class is a fundamental tool of the pyC\(2\)Ray simulation. This object group includes the basic required functions to set up a simulation (e.g., cosmology, time evolution, I/O, raytracing, chemistry, etc.) and access and manage the parameters in the parameter file.
Look at the tutorial on the parameter file for an overview.
pyC\(^2\)Ray provides the basic class, named C2Ray, that is later inherited by the other existing and more extensive class
[3]:
sim = pc2r.C2Ray(paramfile='parameters.yml')
Number of GPUS 1 _________ ____
____ __ __/ ____/__ \ / __ \____ ___ __
/ __ \/ / / / / __/ // /_/ / __ `/ / / /
/ /_/ / /_/ / /___ / __// _, _/ /_/ / /_/ /
/ .___/\__, /\____//____/_/ |_|\__,_/\__, /
/_/ /____/ /____/
GPU Device ID 0: "NVIDIA RTX A1000 6GB Laptop GPU" with compute capability 8.6
Successfully allocated 536.871 Mb of device memory for grid of size N = 256, with source batch size 1
Welcome! Mesh size is N = 256.
Simulation Box size (comoving Mpc): 1.280e+02
Cosmology is on, scaling comoving quantities to the initial redshift, which is z0 = 12.000...
Cosmological parameters used:
h = 0.6766, Tcmb0 = 2.725e+00
Om0 = 0.3097, Ob0 = 0.0490
Using power-law opacity with 10,000 table points between tau=10^(-20) and tau=10^(4)
Using Black-Body sources with effective temperature T = 5.0e+04 K and Radius 1.437e-11 rsun
Spectrum Frequency Range: 3.288e+15 to 1.316e+17 Hz
This is Energy: 1.360e+01 to 5.442e+02 eV
Integrating photoionization rates tables...
INFO: No heating rates
Successfully copied radiation tables to GPU memory.
---- Calculated Clumping Factor (constant model):
min, mean and max clumping : 1.000e+00 1.000e+00 1.000e+00
---- Calculated Mean-Free Path (constant model):
Maximum comoving distance for photons from source mfp = 15.00 cMpc (constant model).
This corresponds to 30.000 grid cells.
Using ASORA Raytracing ( q_max = 52 )
Running in non-MPI (single-GPU/CPU) mode
Starting simulation...
Here is an example of how to access the cosmology class (from astropy) from the pyC\(^2\)Ray simulation class.
[4]:
sim.cosmology.H0
[4]:
$67.66 \; \mathrm{\frac{km}{Mpc\,s}}$
Here is another example of how to set the time step between two redshift steps.
[12]:
z_1, z_2 = 11.5, 11.0
dt = sim.set_timestep(z1=z_1, z2=z_2, num_timesteps=1) * u.s
print(dt.to('Myr'))
t_1 = sim.cosmology.age(z_1).to('Myr')
print(t_1)
t_2 = sim.cosmology.age(z_2).to('Myr')
print(t_1+dt, t_2)
24.678582778205058 Myr
389.6853979573013 Myr
414.36398073550635 Myr 414.36398073550635 Myr
[14]:
sim.write_output??
Signature: sim.write_output(z, ext='.dat')
Docstring:
Write ionization fraction & ionization rates as C2Ray binary files
Parameters
----------
z : float
Redshift (used to name the file)
Source:
def write_output(self, z, ext='.dat'):
"""Write ionization fraction & ionization rates as C2Ray binary files
Parameters
----------
z : float
Redshift (used to name the file)
ext : string
extension of the output file. If '.dat' save a binary file (with tools21cm), otherwise '.npy'.
"""
if(self.rank == 0):
suffix = f"_z{z:.3f}"+ext
if(suffix.endswith('.dat')):
t2c.save_cbin(filename=self.results_basename + "xfrac" + suffix, data=self.xh, bits=64, order='F')
t2c.save_cbin(filename=self.results_basename + "IonRates" + suffix, data=self.phi_ion, bits=32, order='F')
#t2c.save_cbin(filename=self.results_basename + "coldens" + suffix, data=self.coldens, bits=64, order='F')
elif(suffix.endswith('.npy')):
np.save(file=self.results_basename + "xfrac" + suffix, arr=self.xh)
np.save(file=self.results_basename + "IonRates" + suffix, arr=self.phi_ion)
#np.save(file=self.results_basename + "coldens" + suffix, arr=self.coldens)
# print min, max and average quantities
self.printlog('\n--- Reionization History ----')
self.printlog(' min, mean, max xHII : %.5e %.5e %.5e' %(self.xh.min(), self.xh.mean(), self.xh.max()))
self.printlog(' min, mean, max Irate : %.5e %.5e %.5e [1/s]' %(self.phi_ion.min(), self.phi_ion.mean(), self.phi_ion.max()))
self.printlog(' min, mean, max density : %.5e %.5e %.5e [1/cm3]' %(self.ndens.min(), self.ndens.mean(), self.ndens.max()))
# write summary output file
summary_exist = os.path.exists(self.results_basename+'PhotonCounts2.txt')
with open(self.results_basename+'PhotonCounts2.txt', 'a') as f:
if not (summary_exist):
header = '# z\ttot HI atoms\ttot phots\t mean ndens [1/cm3]\t mean Irate [1/s]\tR_mfp [cMpc]\tmean ionization fraction (by volume and mass)\n'
f.write(header)
# mass-average neutral faction
massavrg_ion_frac = np.sum(self.xh*self.ndens)/np.sum(self.ndens)
# calculate total number of neutral hydrogen atoms
tot_nHI = np.sum(self.ndens * (1-self.xh) * self.dr**3)
text = '%.3f\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n' %(z, tot_nHI, self.tot_phots, np.mean(self.ndens), np.mean(self.phi_ion), self.R_max_LLS/self.N*self.boxsize, np.mean(self.xh), massavrg_ion_frac)
f.write(text)
else:
# this is for the other ranks
pass
File: ~/codes/pyC2Ray/pyc2ray/c2ray_base.py
Type: method
Existing Sub-class¶
This tutorial is all about changing the methods of the basic class of the pyC\(^2\)Ray run.
We provide a series of standard class can be C2Ray_Test class. This subclass of the basic class C2Ray is a version used for test simulations and which don’t read N-body input and use simple generated source files.
All the sub-class require a parameter file parameters.yml as input.
c2ray_base.py: implemented the basic functionc2ray_cubep3m.py: specific for CUBEP3M N-bodyc2ray_ramses.py: specific for Ramses hyro N-body simulation… more to come
[2]:
sim = pc2r.C2Ray_Test(paramfile='parameters.yml')
GPU Device 0: "NVIDIA RTX A1000 6GB Laptop GPU" with compute capability 8.6
Succesfully allocated 67.1089 Mb of device memory for grid of size N = 128, with source batch size 1
_________ ____
____ __ __/ ____/__ \ / __ \____ ___ __
/ __ \/ / / / / __/ // /_/ / __ `/ / / /
/ /_/ / /_/ / /___ / __// _, _/ /_/ / /_/ /
/ .___/\__, /\____//____/_/ |_|\__,_/\__, /
/_/ /____/ /____/
Welcome! Mesh size is N = 128.
Simulation Box size (comoving Mpc): 1.400e-02
Cosmology is off.
Using power-law opacity with 10000 table points between tau=10^(-20) and tau=10^(4)
Using Black-Body sources with effective temperature T = 5.0e+04 K and Radius 1.437e-11 rsun
Spectrum Frequency Range: 3.289e+15 to 1.316e+17 Hz
This is Energy: 1.360e+01 to 5.442e+02 eV
Integrating photoionization rates tables...
INFO: No heating rates
Successfully copied radiation tables to GPU memory.
---- Calculated Clumping Factor (constant model):
min, mean and max clumping : 1.000e+00 1.000e+00 1.000e+00
---- Calculated Mean-Free Path (constant model):
Maximum comoving distance for photons from source mfp = 15.00 cMpc (constant model).
This corresponds to 137142.857 grid cells.
Using ASORA Raytracing ( q_max = 193 )
Running in non-MPI (single-GPU/CPU) mode
Starting simulation...
Running: "C2Ray Test"
Write a Sub-class for your Simulation¶
Ideally, you would like to be able to define your simulation class and personalize it based on simulation requirements.
For this reason, we suggest you create your own sub-class by inheriting the pyC\(^2\)Ray base class.
The example below shows how to define a new simulation class named C2Ray_tutorial, which inherits from the base class C2Ray. Here, we add two methods that load a source list, read_sources, and the density field, read_density.
[ ]:
class C2Ray_tutorial(pc2r.c2ray_base.C2Ray):
def __init__(self, paramfile):
"""Basis class for a C2Ray Simulation
Parameters
----------
paramfile : str
Name of a YAML file containing parameters for the C2Ray simulation
"""
super().__init__(paramfile)
self.printlog('Running: "C2Ray tutorial for %d Mpc/h volume"' %self.boxsize)
# ===========================================
# HEREAFTER: USER DEFINED METHODS
# ===========================================
def read_sources(self, z, nsrc, dt):
np.random.seed(918)
# Read random sources (e.g.: *.npy, *.h5, etc.)
pos_halo = np.random.uniform(low=0, high=sim.boxsize, size=(nsrc, 3))
mhalo = np.random.uniform(1e8, 1e14, nsrc)*u.Msun
# Define stellar-to-halo relation
fstar = 0.1
# Define escaping fraction
fesc = 0.1
# sum togheter the star mass for sources within the same voxel
pos_star, mstar = pc2r.other_utils.bin_sources(srcpos_mpc=pos_halo, mstar_msun=mhalo*fstar*fesc, boxsize=sim.boxsize, meshsize=sim.N)
"""
pos_star = np.array([sim.N//2, sim.N//2, sim.N//2])
pos_star = pos_star[None,...]
mstar = np.array([1e14])
"""
# this reference flux is necessary only for a numercial reason
S_star_ref = 1e48
# The normalize flux in CGS units
dotN = (mstar*u.Msun/(cst.m_p*dt)).cgs.value
# calculate some quantity thtat you want to print (e.g. total number of ionizing photons)
self.tot_phots = np.sum(dotN * dt)
return pos_star, dotN/S_star_ref
def read_density(self, z):
# Read the density field
self.ndens = 1e-6 * np.ones((sim.N, sim.N, sim.N))
return self.ndens
[9]:
paramfile = './parameters.yml'
# init the C2Ray class for the tutorial
sim = C2Ray_tutorial(paramfile)
Number of GPUS 1
_________ ____
____ __ __/ ____/__ \ / __ \____ ___ __
/ __ \/ / / / / __/ // /_/ / __ `/ / / /
/ /_/ / /_/ / /___ / __// _, _/ /_/ / /_/ /
/ .___/\__, /\____//____/_/ |_|\__,_/\__, /
/_/ /____/ /____/
GPU Device ID 0: "NVIDIA RTX A1000 6GB Laptop GPU" with compute capability 8.6
Successfully allocated 536.871 Mb of device memory for grid of size N = 256, with source batch size 1
Welcome! Mesh size is N = 256.
Simulation Box size (comoving Mpc): 1.280e+02
Cosmology is on, scaling comoving quantities to the initial redshift, which is z0 = 12.000...
Cosmological parameters used:
h = 0.6766, Tcmb0 = 2.725e+00
Om0 = 0.3097, Ob0 = 0.0490
Using power-law opacity with 10,000 table points between tau=10^(-20) and tau=10^(4)
Using Black-Body sources with effective temperature T = 5.0e+04 K and Radius 1.437e-11 rsun
Spectrum Frequency Range: 3.288e+15 to 1.316e+17 Hz
This is Energy: 1.360e+01 to 5.442e+02 eV
Integrating photoionization rates tables...
INFO: No heating rates
Successfully copied radiation tables to GPU memory.
---- Calculated Clumping Factor (constant model):
min, mean and max clumping : 1.000e+00 1.000e+00 1.000e+00
---- Calculated Mean-Free Path (constant model):
Maximum comoving distance for photons from source mfp = 15.00 cMpc (constant model).
This corresponds to 30.000 grid cells.
Using ASORA Raytracing ( q_max = 52 )
Running in non-MPI (single-GPU/CPU) mode
Starting simulation...
Running: "C2Ray tutorial for 128 Mpc/h volume"
[12]:
# Read homogeneous density field
ndens = sim.read_density(z=7.0)
[17]:
# Read source files
srcpos, normflux = sim.read_sources(nsrc=10, z=7.0, dt=10.*u.Myr)
print(srcpos)
[[ 15 187 45]
[ 15 218 1]
[ 33 93 31]
[ 50 28 170]
[ 96 156 213]
[185 199 122]
[192 232 217]
[224 121 197]
[243 150 229]
[255 90 252]]
[ ]: