Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Non cubic #289

Merged
merged 27 commits into from
May 16, 2023
Merged

Non cubic #289

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
ccb330d
Add in variable NON_CUBIC_FACTOR plus some checks to allow non-cubic …
BradGreig Mar 30, 2022
8b68d07
Add in new variables and update macros
BradGreig Mar 30, 2022
f14ca75
Update initial conditions, plus some relevant functions and outputs f…
BradGreig Mar 30, 2022
9946de9
Update ionisation box to have non-cubic boxes
BradGreig Mar 30, 2022
64838ef
Update the perturbed field (and haloes) quantities for non-cubic boxes
BradGreig Mar 30, 2022
1386e97
Partial update to find haloes. Needs more work for identifying haloes
BradGreig Mar 30, 2022
5f1d733
Add non-cubic to spin temperature
BradGreig Mar 30, 2022
dc23520
Add non-cubic functionality to brightness temperature box
BradGreig Mar 30, 2022
1cdb935
Update header to add new variable from UserParams
BradGreig Mar 30, 2022
9f1e29c
fix: merge master
steven-murray Jun 2, 2022
e8ffd76
merge master
steven-murray May 5, 2023
1b630e9
fix: stray merge error
steven-murray May 5, 2023
64bcedd
test: coverage for non_cubic_factor
steven-murray May 5, 2023
b505186
test: fix test invocation
steven-murray May 8, 2023
14da932
test: allow higher abs tolerance in coeval_vs_low_level
steven-murray May 8, 2023
40fa3f3
test: reduce required precision on equivalence between coeval/low-level
steven-murray May 9, 2023
82dd26c
Remove first TODO with respect to identifying haloes
BradGreig May 12, 2023
f1702a3
Remove TODO with painting ionisation map with sphere method
BradGreig May 12, 2023
69f4acc
Missing factors probably from merge
BradGreig May 12, 2023
d3202c5
Update summary functions for non-cubic boxes
BradGreig May 12, 2023
caca97b
NON_CUBIC factor is a float not an int
BradGreig May 12, 2023
3bfb4ba
Missing line-of-sight factor of HII_D_PARA
BradGreig May 12, 2023
386605a
Halo field tests were failing, so updated the data following changes …
BradGreig May 12, 2023
62309ad
Mathematical operations may require these quantities to be signed
BradGreig May 16, 2023
14a419a
Changing the data back following fix in the previous commit
BradGreig May 16, 2023
9611d41
Weirdly the test data did not get correctly updated. Now is
BradGreig May 16, 2023
8f054ac
The variable ct should be private
BradGreig May 16, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 18 additions & 2 deletions src/py21cmfast/inputs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -409,6 +409,9 @@ class UserParams(StructWithDefaults):
Number of cells for the high-res box (sampling ICs) along a principal axis. To avoid
sampling issues, DIM should be at least 3 or 4 times HII_DIM, and an integer multiple.
By default, it is set to 3*HII_DIM.
NON_CUBIC_FACTOR : float, optional
Factor which allows the creation of non-cubic boxes. It will shorten/lengthen the line
of sight dimension of all boxes. NON_CUBIC_FACTOR * DIM/HII_DIM must result in an integer
BOX_LEN : float, optional
Length of the box, in Mpc. Default 300 Mpc.
HMF: int or str, optional
Expand Down Expand Up @@ -464,6 +467,7 @@ class UserParams(StructWithDefaults):
"BOX_LEN": 300.0,
"DIM": None,
"HII_DIM": 200,
"NON_CUBIC_FACTOR": 1.0,
"USE_FFTW_WISDOM": False,
"HMF": 1,
"USE_RELATIVE_VELOCITIES": False,
Expand Down Expand Up @@ -500,15 +504,27 @@ def DIM(self):
"""Number of cells for the high-res box (sampling ICs) along a principal axis."""
return self._DIM or 3 * self.HII_DIM

@property
def NON_CUBIC_FACTOR(self):
"""Factor to shorten/lengthen the line-of-sight dimension (non-cubic boxes)."""
dcf = self.DIM * self._NON_CUBIC_FACTOR
hdcf = self.HII_DIM * self._NON_CUBIC_FACTOR
if dcf % int(dcf) or hdcf % int(hdcf):
raise ValueError(
"NON_CUBIC_FACTOR * DIM and NON_CUBIC_FACTOR * HII_DIM must be integers"
)
else:
return self._NON_CUBIC_FACTOR

@property
def tot_fft_num_pixels(self):
"""Total number of pixels in the high-res box."""
return self.DIM**3
return self.NON_CUBIC_FACTOR * self.DIM**3

@property
def HII_tot_num_pixels(self):
"""Total number of pixels in the low-res box."""
return self.HII_DIM**3
return self.NON_CUBIC_FACTOR * self.HII_DIM**3

@property
def POWER_SPECTRUM(self):
Expand Down
27 changes: 20 additions & 7 deletions src/py21cmfast/outputs.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,12 @@ def prepare_for_spin_temp(self, flag_options: FlagOptions, force: bool = False):
self.prepare(keep=keep, force=force)

def _get_box_structures(self) -> dict[str, dict | tuple[int]]:
shape = (self.user_params.HII_DIM,) * 3
hires_shape = (self.user_params.DIM,) * 3
shape = (self.user_params.HII_DIM,) * 2 + (
int(self.user_params.NON_CUBIC_FACTOR * self.user_params.HII_DIM),
)
hires_shape = (self.user_params.DIM,) * 2 + (
int(self.user_params.NON_CUBIC_FACTOR * self.user_params.DIM),
)

out = {
"lowres_density": shape,
Expand Down Expand Up @@ -204,8 +208,10 @@ class PerturbedField(_OutputStructZ):

def _get_box_structures(self) -> dict[str, dict | tuple[int]]:
return {
"density": (self.user_params.HII_DIM,) * 3,
"velocity": (self.user_params.HII_DIM,) * 3,
"density": (self.user_params.HII_DIM,) * 2
+ (int(self.user_params.NON_CUBIC_FACTOR * self.user_params.HII_DIM),),
"velocity": (self.user_params.HII_DIM,) * 2
+ (int(self.user_params.NON_CUBIC_FACTOR * self.user_params.HII_DIM),),
}

def get_required_input_arrays(self, input_box: _BaseOutputStruct) -> list[str]:
Expand Down Expand Up @@ -397,7 +403,9 @@ def __init__(
super().__init__(**kwargs)

def _get_box_structures(self) -> dict[str, dict | tuple[int]]:
shape = (self.user_params.HII_DIM,) * 3
shape = (self.user_params.HII_DIM,) * 2 + (
int(self.user_params.NON_CUBIC_FACTOR * self.user_params.HII_DIM),
)
return {
"Ts_box": shape,
"x_e_box": shape,
Expand Down Expand Up @@ -520,7 +528,9 @@ def _get_box_structures(self) -> dict[str, dict | tuple[int]]:
else:
n_filtering = 1

shape = (self.user_params.HII_DIM,) * 3
shape = (self.user_params.HII_DIM,) * 2 + (
int(self.user_params.NON_CUBIC_FACTOR * self.user_params.HII_DIM),
)
filter_shape = (n_filtering,) + shape

out = {
Expand Down Expand Up @@ -619,7 +629,10 @@ class BrightnessTemp(_AllParamsBox):
_filter_params = _OutputStructZ._filter_params

def _get_box_structures(self) -> dict[str, dict | tuple[int]]:
return {"brightness_temp": (self.user_params.HII_DIM,) * 3}
return {
"brightness_temp": (self.user_params.HII_DIM,) * 2
+ (int(self.user_params.NON_CUBIC_FACTOR * self.user_params.HII_DIM),)
}

@cached_property
def global_Tb(self):
Expand Down
1 change: 1 addition & 0 deletions src/py21cmfast/src/21cmFAST.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ struct UserParams{
int HII_DIM;
int DIM;
float BOX_LEN;
float NON_CUBIC_FACTOR;
bool USE_FFTW_WISDOM;
int HMF;
int USE_RELATIVE_VELOCITIES;
Expand Down
50 changes: 25 additions & 25 deletions src/py21cmfast/src/BrightnessTemperatureBox.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
float *x_pos_offset = calloc(astro_params->N_RSD_STEPS,sizeof(float));
float **delta_T_RSD_LOS = (float **)calloc(user_params->N_THREADS,sizeof(float *));
for(i=0;i<user_params->N_THREADS;i++) {
delta_T_RSD_LOS[i] = (float *)calloc(user_params->HII_DIM,sizeof(float));
delta_T_RSD_LOS[i] = (float *)calloc(HII_D_PARA,sizeof(float));
}


Expand All @@ -39,7 +39,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
#pragma omp for
for (i=0; i<user_params->HII_DIM; i++){
for (j=0; j<user_params->HII_DIM; j++){
for (k=0; k<user_params->HII_DIM; k++){
for (k=0; k<HII_D_PARA; k++){
*((float *)v + HII_R_FFT_INDEX(i,j,k)) = perturb_field->velocity[HII_R_INDEX(i,j,k)];
}
}
Expand Down Expand Up @@ -69,7 +69,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
#pragma omp for reduction(+:ave)
for (i=0; i<user_params->HII_DIM; i++){
for (j=0; j<user_params->HII_DIM; j++){
for (k=0; k<user_params->HII_DIM; k++){
for (k=0; k<HII_D_PARA; k++){

pixel_deltax = perturb_field->density[HII_R_INDEX(i,j,k)];
pixel_x_HI = ionized_box->xH_box[HII_R_INDEX(i,j,k)];
Expand Down Expand Up @@ -115,7 +115,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct

memcpy(vel_gradient, v, sizeof(fftwf_complex)*HII_KSPACE_NUM_PIXELS);

dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, user_params->N_THREADS, vel_gradient);
dft_r2c_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, vel_gradient);

#pragma omp parallel shared(vel_gradient) private(n_x,n_y,n_z,k_x,k_y,k_z) num_threads(user_params->N_THREADS)
{
Expand All @@ -132,8 +132,8 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
else
k_y = n_y * DELTA_K;

for (n_z=0; n_z<=HII_MIDDLE; n_z++){
k_z = n_z * DELTA_K;
for (n_z=0; n_z<=HII_MIDDLE_PARA; n_z++){
k_z = n_z * DELTA_K_PARA;

// take partial deriavative along the line of sight
*((fftwf_complex *) vel_gradient + HII_C_INDEX(n_x,n_y,n_z)) *= k_z*I/(float)HII_TOT_NUM_PIXELS;
Expand All @@ -142,7 +142,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
}
}

dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, user_params->N_THREADS, vel_gradient);
dft_c2r_cube(user_params->USE_FFTW_WISDOM, user_params->HII_DIM, HII_D_PARA, user_params->N_THREADS, vel_gradient);

// now add the velocity correction to the delta_T maps (only used for T_S >> T_CMB case).
max_v_deriv = fabs(global_params.MAX_DVDR*H);
Expand All @@ -158,7 +158,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
#pragma omp for
for (i=0; i<user_params->HII_DIM; i++){
for (j=0; j<user_params->HII_DIM; j++){
for (k=0; k<user_params->HII_DIM; k++){
for (k=0; k<HII_D_PARA; k++){

gradient_component = fabs(vel_gradient[HII_R_FFT_INDEX(i,j,k)]/H + 1.0);

Expand Down Expand Up @@ -217,17 +217,17 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
for (j=0; j<user_params->HII_DIM; j++){

// Generate the optical-depth for the specific line-of-sight with R.S.D
for(k=0;k<user_params->HII_DIM;k++) {
for(k=0;k<HII_D_PARA;k++) {
delta_T_RSD_LOS[omp_get_thread_num()][k] = 0.0;
}

for (k=0; k<user_params->HII_DIM; k++){
for (k=0; k<HII_D_PARA; k++){

if((fabs(box->brightness_temp[HII_R_INDEX(i,j,k)]) >= FRACT_FLOAT_ERR) && \
(ionized_box->xH_box[HII_R_INDEX(i,j,k)] >= FRACT_FLOAT_ERR)) {

if(k==0) {
d1_low = v[HII_R_FFT_INDEX(i,j,user_params->HII_DIM-1)]/H;
d1_low = v[HII_R_FFT_INDEX(i,j,HII_D_PARA-1)]/H;
d2_low = v[HII_R_FFT_INDEX(i,j,k)]/H;
}
else {
Expand All @@ -236,7 +236,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
}

// Displacements (converted from velocity) for the original cell centres straddling half of the sub-cells (cell after)
if(k==(user_params->HII_DIM-1)) {
if(k==(HII_D_PARA-1)) {
d1_high = v[HII_R_FFT_INDEX(i,j,k)]/H;
d2_high = v[HII_R_FFT_INDEX(i,j,0)]/H;
}
Expand Down Expand Up @@ -285,7 +285,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct

// check if the new cell position is at the edge of the box. If so, periodic boundary conditions
if(k<((int)cell_distance+1)) {
delta_T_RSD_LOS[omp_get_thread_num()][k-((int)cell_distance+1) + user_params->HII_DIM] += \
delta_T_RSD_LOS[omp_get_thread_num()][k-((int)cell_distance+1) + HII_D_PARA] += \
box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
else {
Expand All @@ -303,7 +303,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct

// Check if the first part of the sub-cell is at the box edge
if(k<(((int)cell_distance))) {
delta_T_RSD_LOS[omp_get_thread_num()][k-((int)cell_distance) + user_params->HII_DIM] += \
delta_T_RSD_LOS[omp_get_thread_num()][k-((int)cell_distance) + HII_D_PARA] += \
fraction_within*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
else {
Expand All @@ -312,7 +312,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
}
// Check if the second part of the sub-cell is at the box edge
if(k<(((int)cell_distance + 1))) {
delta_T_RSD_LOS[omp_get_thread_num()][k-((int)cell_distance+1) + user_params->HII_DIM] += \
delta_T_RSD_LOS[omp_get_thread_num()][k-((int)cell_distance+1) + HII_D_PARA] += \
fraction_outside*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
else {
Expand All @@ -330,7 +330,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct

// Check the periodic boundaries conditions and move the fraction of each sub-cell to the appropriate new cell
if(k==0) {
delta_T_RSD_LOS[omp_get_thread_num()][user_params->HII_DIM-1] += \
delta_T_RSD_LOS[omp_get_thread_num()][HII_D_PARA-1] += \
fraction_outside*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
delta_T_RSD_LOS[omp_get_thread_num()][k] += \
fraction_within*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
Expand All @@ -350,7 +350,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
fraction_within = 1. - fraction_outside;

// Check the periodic boundaries conditions and move the fraction of each sub-cell to the appropriate new cell
if(k==(user_params->HII_DIM-1)) {
if(k==(HII_D_PARA-1)) {
delta_T_RSD_LOS[omp_get_thread_num()][k] += \
fraction_within*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
delta_T_RSD_LOS[omp_get_thread_num()][0] += \
Expand All @@ -373,8 +373,8 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
// sub-cell is entirely contained within the new cell (just add it to the new cell)

// check if the new cell position is at the edge of the box. If so, periodic boundary conditions
if(k>(user_params->HII_DIM - 1 - (int)cell_distance)) {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance - user_params->HII_DIM] += \
if(k>(HII_D_PARA - 1 - (int)cell_distance)) {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance - HII_D_PARA] += \
box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
else {
Expand All @@ -390,17 +390,17 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
fraction_within = 1. - fraction_outside;

// Check if the first part of the sub-cell is at the box edge
if(k>(user_params->HII_DIM - 1 - ((int)cell_distance-1))) {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance-1 - user_params->HII_DIM] += \
if(k>(HII_D_PARA - 1 - ((int)cell_distance-1))) {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance-1 - HII_D_PARA] += \
fraction_within*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
else {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance-1] += \
fraction_within*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
// Check if the second part of the sub-cell is at the box edge
if(k>(user_params->HII_DIM - 1 - ((int)cell_distance))) {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance - user_params->HII_DIM] += \
if(k>(HII_D_PARA - 1 - ((int)cell_distance))) {
delta_T_RSD_LOS[omp_get_thread_num()][k+(int)cell_distance - HII_D_PARA] += \
fraction_outside*box->brightness_temp[HII_R_INDEX(i,j,k)]/((float)astro_params->N_RSD_STEPS);
}
else {
Expand All @@ -413,7 +413,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
}
}

for(k=0;k<user_params->HII_DIM;k++) {
for(k=0;k<HII_D_PARA;k++) {
box->brightness_temp[HII_R_INDEX(i,j,k)] = delta_T_RSD_LOS[omp_get_thread_num()][k];

ave += delta_T_RSD_LOS[omp_get_thread_num()][k];
Expand All @@ -430,7 +430,7 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
#pragma omp for reduction(+:ave)
for (i=0; i<user_params->HII_DIM; i++){
for (j=0; j<user_params->HII_DIM; j++){
for (k=0; k<user_params->HII_DIM; k++){
for (k=0; k<HII_D_PARA; k++){

dvdx = vel_gradient[HII_R_FFT_INDEX(i,j,k)];

Expand Down
Loading