21cmfast / 21cmFAST

Official repository for 21cmFAST: a code for generating fast simulations of the cosmological 21cm signal
MIT License
58 stars 37 forks source link

halos popping in/out from one redshift to next with N_POISSON < 0 [BUG] #294

Closed kenj19 closed 2 years ago

kenj19 commented 2 years ago

Describe the bug:

Halos stochastically popping in and out of existence from one redshift to the next, despite setting global_params.N_POISSON to a negative integer.

To Reproduce: Run determine_halo_list at z=7.0, 8.0 and populate box (of DIM^3) with halo counts. Found that some halos present at z=8.0 are no longer present at z=7.0.

import py21cmfast as p21c
import numpy as np

print(f"\n ============= Using 21cmFAST version {p21c.__version__} ============== ")

BOX_LEN = 50
HII_DIM = 50
DIM = 150
rseed = 233031

high_z, low_z = 8.0, 7.0

with p21c.global_params.use(N_POISSON=-1):

    print(f"\n ============= Using N_POISSON = {p21c.global_params.N_POISSON} ============== ")

    user_params = p21c.UserParams(
        BOX_LEN=BOX_LEN,
        DIM=DIM,
        HII_DIM=HII_DIM,
        USE_INTERPOLATION_TABLES=True
    )

    init_cond = p21c.initial_conditions(
        user_params=user_params,
        random_seed=rseed,
    )

    high_z_halo_field = p21c.determine_halo_list(
        redshift=high_z,
        init_boxes=init_cond,
        user_params=user_params,
    )

    low_z_halo_field = p21c.determine_halo_list(
        redshift=low_z, 
        init_boxes=init_cond, 
        user_params=user_params,
    )

high_z_halo_coords, low_z_halo_coords = high_z_halo_field.halo_coords, low_z_halo_field.halo_coords 

Now, populate boxes with number of halos (Nh).

high_z_Nh, low_z_Nh = np.zeros(shape=(DIM,DIM,DIM)), np.zeros(shape=(DIM,DIM,DIM))

for i in range(len(high_z_halo_coords)):
    high_z_Nh[tuple(high_z_halo_coords[i])] += 1

for j in range(len(low_z_halo_coords)):
    low_z_Nh[tuple(low_z_halo_coords[j])] += 1

Difference low_z and high_z boxes. If any voxels < 0, halos at higher redshift pop out of existence at lower redshift.

diff_Nh = low_z_Nh-high_z_Nh

print(f"\n ============= Number of halos present at z=8.0, but not z=7.0: {len(diff_Nh[diff_Nh<0])} ============== \n ")

The following outputs are produced:

============= Using 21cmFAST version 3.1.5 ============== 

 ============= Using N_POISSON = -1 ============== 

 ============= Number of halos present at z=8.0, but not z=7.0: 3616 ==============

Expected behavior:

I expected that there would be no voxels with negative values for N_POISSON < 0, as halos present at higher redshift should also be present at lower redshift.

Details: OS: macOS Monterey (v12.4) Python: v3.9.13 21cmFAST: v3.1.5

BradGreig commented 2 years ago

Hi @kenj19, I do not think this is the correct way to identify haloes across redshift snapshots. What I think you are finding are haloes that have moved into neighbouring cells between the two redshift snapshots.

The halo finder is applied on the density field at the specific redshift. This is the evolved density field (using 2LPT) thus across the two snapshots the density field will have evolved. Thus, the locations of the halos will also evolve.

This movement of the halos could then cause two things to happen. Either: (i) it moves to a voxel that previously did not have a halo or (ii) it moves into a voxel which already contains a halo. In the former case, you'll just have the halo with a position offset, in the latter it will be the mass of the most massive halo in the cell. Thus, if it is smaller that the existing halo it will not be identifiable.

I suspect this is what you are finding.

kenj19 commented 2 years ago

Hi @BradGreig, thank you for the response - this resolves my misunderstanding. Thanks!