Input and Output Files

COLT uses the Hierarchical Data Format (HDF5) library for reading and writing files (with C++ wrappers). This page explains the general structure and philosophy of initial conditions and simulation output files.

Supported Geometries

  • Plane Parallel Slab (slab)

  • Spherical Shells (spherical)

  • 3D Cartesian grid (cartesian)

  • Adaptive Resolution Octree (octree)

  • Unstructured Voronoi Mesh (voronoi)

  • Secondary Outputs (module specific)

Initial Conditions Files

The following example code outlines the general procedure for writing a COLT initial conditions file. Specific use cases can vary significantly, and many datasets are optional or can be overridden by configuration options. The presentation shows functional pseudo code with if statements to clarify geometry or module specific options.

# Write colt file (generally uses cgs units)
filename = f'{colt_dir}/colt_{snap:03d}.hdf5'
with h5py.File(filename, 'w') as f:
    # Simulation properties
    if cosmological:
        f.attrs['redshift'] = z     # Current simulation redshift
        f.attrs['Omega0'] = Omega0  # Matter density [rho_crit_0]
        f.attrs['OmegaB'] = OmegaB  # Baryon density [rho_crit_0]
        f.attrs['h100'] = h         # Hubble constant [100 km/s/Mpc]
    else:
        f.attrs['time'] = time      # Simulation time [s]

    # Geometry properties
    # Note: Due to camera positioning, it is most convenient to center the box at [0,0,0]
    f.attrs['n_cells'] = np.int32(n_cells) # Total number of cells (optional)
    if geometry == 'cartesian':
        f.attrs['nx'] = np.int32(nx) # Number of x cells
        f.attrs['ny'] = np.int32(ny) # Number of y cells
        f.attrs['nz'] = np.int32(nz) # Number of z cells
        R = ...                     # Box radius [cm]
        f.attrs['r_box'] = R        # For irregular boxes use bbox instead
        # bbox = np.array([[-R,-R,-R],[R,R,R]], dtype=np.float64)
        # f.create_dataset('bbox', data=bbox) # Bounding box [cm]
        # f['bbox'].attrs['units'] = b'cm'
    elif geometry == 'octree':
        r = ...                     # Cell left corner positions [cm]
        w = ...                     # Cell widths [cm]
        f.create_dataset('r', data=r, dtype=np.float64)
        f['r'].attrs['units'] = b'cm'
        f.create_dataset('w', data=w, dtype=np.float64)
        f['w'].attrs['units'] = b'cm'
        parent = ...                # Cell parent indices
        child = ...                 # Cell child indices
        child_check = ...           # True if cell has children
        f.create_dataset('parent', data=parent, dtype=np.int32)
        f.create_dataset('child', data=child, dtype=np.int32)
        f.create_dataset('child_check', data=child_check, dtype=np.int32)
    elif geometry == 'slab':
        r_edges = ...               # Cell edge positions (1D) [cm]
        f.create_dataset('r_edges', data=r_edges, dtype=np.float64)
        f['r_edges'].attrs['units'] = b'cm'
    elif geometry == 'spherical':
        r_edges = ...               # Shell edge positions (1D) [cm]
        f.create_dataset('r_edges', data=r_edges, dtype=np.float64)
        f['r_edges'].attrs['units'] = b'cm'
    elif geometry == 'voronoi':
        R = ...                     # Box radius [cm]
        f.attrs['r_box'] = R        # For irregular boxes use bbox instead
        # bbox = np.array([[-R,-R,-R],[R,R,R]], dtype=np.float64)
        # f.create_dataset('bbox', data=bbox) # Bounding box [cm]
        # f['bbox'].attrs['units'] = b'cm'
        # If neither is specified then bbox is inferred from particle positions
        r = ...                     # Mesh generating points [cm]
        f.create_dataset('r', data=r, dtype=np.float64)
        f['r'].attrs['units'] = b'cm'

    # Gas properties
    v = ...                         # Velocities [cm/s] (optional)
    f.create_dataset('v', data=v, dtype=np.float64)
    f['v'].attrs['units'] = b'cm/s'
    if set_density_from_mass:
        m = ...                     # Masses [g]
        f.create_dataset('m', data=m, dtype=np.float64)
        f['m'].attrs['units'] = b'g'
    else:
        rho = ...                   # Densities [g/cm^3]
        f.create_dataset('rho', data=rho, dtype=np.float64)
        f['rho'].attrs['units'] = b'g/cm^3'
    if use_internal_energy:
        e_int = ...                 # Internal energies [cm^2/s^2]
        f.create_dataset('e_int', data=e_int, dtype=np.float64)
        f['e_int'].attrs['units'] = b'cm^2/s^2'
    else:
        T = ...                     # Temperatures [K]
        f.create_dataset('T', data=T, dtype=np.float64)
        f['T'].attrs['units'] = b'K'
    x_HI = ...                      # Neutral hydrogen fraction (optional)
    x_HII = ...                     # Ionized hydrogen fraction (optional)
    x_H2 = ...                      # Molecular hydrogen fraction (optional)
    x_HeI = ...                     # Neutral helium fraction (optional)
    x_HeII = ...                    # Ionized helium fraction (optional)
    x_e = ...                       # Electron fraction (optional)
    f.create_dataset('x_HI', data=x_HI, dtype=np.float64)
    f.create_dataset('x_HII', data=x_HII, dtype=np.float64)
    f.create_dataset('x_H2', data=x_H2, dtype=np.float64)
    f.create_dataset('x_HeI', data=x_HeI, dtype=np.float64)
    f.create_dataset('x_HeII', data=x_HeII, dtype=np.float64)
    f.create_dataset('x_e', data=x_e, dtype=np.float64)
    Z = ...                         # Metallicity [mass fraction] (optional)
    D = ...                         # Dust-to-gas ratio [mass fraction] (optional)
    f.create_dataset('Z', data=Z, dtype=np.float64)
    f.create_dataset('D', data=D, dtype=np.float64)
    # B = ...                         # Magnetic field [Gauss] (optional)
    # f.create_dataset('B', data=B, dtype=np.float64)
    # f['B'].attrs['units'] = b'G'
    # e_int = ...                     # Specific internal energy [cm^2/s^2] (optional)
    # f.create_dataset('e_int', data=e_int, dtype=np.float64)
    # f['e_int'].attrs['units'] = b'cm^2/s^2'
    # SFR = ...                       # Star formation rate [Msun/yr]
    # f.create_dataset('SFR', data=SFR, dtype=np.float64)
    # f['SFR'].attrs['units'] = b'Msun/yr'

    # Star properties
    f.attrs['n_stars'] = n_stars
    r_star = ...                    # Star positions [cm]
    f.create_dataset('r_star', data=r_star, dtype=np.float64)
    f['r_star'].attrs['units'] = b'cm'
    v_star = ...                    # Star velocities [cm]
    f.create_dataset('v_star', data=v_star, dtype=np.float64)
    f['v_star'].attrs['units'] = b'cm/s'
    m_init_star = ...               # Star initial masses [Msun]
    f.create_dataset('m_init_star', data=m_init_star, dtype=np.float64)
    f['m_init_star'].attrs['units'] = b'Msun'
    Z_star = ...                    # Star metallicity [mass fraction]
    f.create_dataset('Z_star', data=Z_star, dtype=np.float64)
    age_star = ...                  # Star ages [Gyr]
    f.create_dataset('age_star', data=age_star, dtype=np.float64)
    f['age_star'].attrs['units'] = b'Gyr'

Example Configuration File

While the format for configuration files was already described, we now provide a working example configuration file for H-alpha Monte Carlo radiative transfer.

# COLT config file
--- !mcrt                     # Monte Carlo radiative transfer module

init_dir: ics                 # Initial conditions directory (default: "ics")
init_base: colt               # Initial conditions base name (optional)
output_dir: Ha                # Output directory name (default: "output")
output_base: Ha               # Output file base name (default: "colt")

recombinations: true          # Include recombination emission
collisions: true              # Include collisional excitation

cosmological: true            # Indicates whether the simulation is cosmological
line: Balmer-alpha            # Name of the line (default "Lyman-alpha")
v_turb_kms: 10                # Microturbulent velocity [km/s]
dust_model: MW                # Dust model: SMC, MW, etc.

# Information about sources
n_photons: 1000000            # Number of photon packets (increase for production runs)
j_exp: 0.75                   # Luminosity boosting exponent

# Information about escape
output_photons: false         # Output escaped photon packets
spherical_escape: true        # Photons escape from a sphere
escape_radius_bbox: 0.9       # Escape radius relative to the bbox
emission_radius_bbox: 0.75    # Emission radius relative to the bbox

# Information about cameras
freq_min: -500                # Frequency minimum [km/s]
freq_max: 500                 # Frequency maximum [km/s]
n_bins: 500                   # Number of frequency bins
image_radius_bbox: 0.5        # Image radius relative to the bbox (1 Rvir)
n_pixels: 512                 # Number of image pixels
cameras:                      # Add camera directions manually
  - [0,0,1]
  - [1,0,0]
output_freq_stds: true        # Output frequency standard deviations (default: false)
output_freq2_images: true     # Frequency moment images (default: false)

Output Files

The HDF5 file format allows self-describing output files. More specific details are coming soon.