Skip to content

Commit

Permalink
fix: NCF updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@steven-murray
steven-murray committed Dec 19, 2023
1 parent 4580080 commit a9b3b72
Show file tree
Hide file tree
Showing 3 changed files with 21 additions and 25 deletions.
10 changes: 4 additions & 6 deletions src/py21cmfast/src/BrightnessTemperatureBox.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,6 @@ int ComputeBrightnessTemp(float redshift, struct UserParams *user_params, struct
for (i=0; i<user_params->HII_DIM; i++){
for (j=0; j<user_params->HII_DIM; j++){
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)];
for (k=0; k<user_params->HII_DIM; k++){
*((float *)v + HII_R_FFT_INDEX(i,j,k)) = perturb_field->velocity_z[HII_R_INDEX(i,j,k)];
}
}
Expand Down Expand Up @@ -127,7 +125,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++){
dvdx = clip(vel_gradient[HII_R_FFT_INDEX(i,j,k)], -max_v_deriv, max_v_deriv);
box->brightness_temp[HII_R_INDEX(i,j,k)] /= (dvdx/H + 1.0);
}
Expand Down Expand Up @@ -208,7 +206,7 @@ void get_velocity_gradient(struct UserParams *user_params, float *v, float *vel_
{
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);

float k_x, k_y, k_z;
int n_x, n_y, n_z;
Expand All @@ -227,7 +225,7 @@ void get_velocity_gradient(struct UserParams *user_params, float *v, float *vel_
else
k_y = n_y * DELTA_K;

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

// take partial deriavative along the line of sight
Expand All @@ -237,5 +235,5 @@ void get_velocity_gradient(struct UserParams *user_params, float *v, float *vel_
}
}

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);
}
4 changes: 1 addition & 3 deletions src/py21cmfast/src/PerturbField.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -984,8 +984,6 @@ int ComputePerturbField(
for (i=0; i<user_params->HII_DIM; i++){
for (j=0; j<user_params->HII_DIM; j++){
for (k=0; k<HII_D_PARA; k++){
*((float *)perturbed_field->velocity + HII_R_INDEX(i,j,k)) = *((float *)HIRES_density_perturb + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5), (unsigned long long)(j*f_pixel_factor+0.5), (unsigned long long)(k*f_pixel_factor+0.5)));
for (k=0; k<user_params->HII_DIM; k++){
*((float *)perturbed_field->velocity_z + HII_R_INDEX(i,j,k)) = *((float *)HIRES_density_perturb + R_FFT_INDEX((unsigned long long)(i*f_pixel_factor+0.5), (unsigned long long)(j*f_pixel_factor+0.5), (unsigned long long)(k*f_pixel_factor+0.5)));
}
}
Expand All @@ -1000,7 +998,7 @@ int ComputePerturbField(
#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<user_params->HII_D_PARA; k++){
*((float *)perturbed_field->velocity_z + HII_R_INDEX(i,j,k)) = *((float *)LOWRES_density_perturb + HII_R_FFT_INDEX(i,j,k));
}
}
Expand Down
32 changes: 16 additions & 16 deletions src/py21cmfast/src/subcell_rsds.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -70,17 +70,17 @@ void apply_subcell_rsds(
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 @@ -89,7 +89,7 @@ void apply_subcell_rsds(
}

// 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 @@ -138,7 +138,7 @@ void apply_subcell_rsds(

// 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 @@ -156,7 +156,7 @@ void apply_subcell_rsds(

// 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 @@ -165,7 +165,7 @@ void apply_subcell_rsds(
}
// 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 @@ -183,7 +183,7 @@ void apply_subcell_rsds(

// 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 @@ -203,7 +203,7 @@ void apply_subcell_rsds(
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 @@ -226,8 +226,8 @@ void apply_subcell_rsds(
// 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 @@ -243,17 +243,17 @@ void apply_subcell_rsds(
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 @@ -266,7 +266,7 @@ void apply_subcell_rsds(
}
}

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 Down

0 comments on commit a9b3b72

Please sign in to comment.