Always filter emissivity fields - #755
Open
jordanflitter wants to merge 16 commits into
Open
Conversation
Prepared HaloBox to work also with Eulerian source models. If we work with Eulerian source models, then we use the Eulerian perturbed density field rather than the Lagrangian density field from the initial conditions, and we do not advect it into its final Eulerian positions as we do for the Lagrangian source models. For E-INTEGRAL, the SFRD and X-ray emissivity are computed very similarly to what is already found in RadiationFields.c. Note that the X-ray emissivity is computed differently between Lagrangian and Eulerian source models, so special care was given to differentiate between these two types of source models. Also, we don’t compute the SFRD directly for CONST-ION-EFF, but rather evaluate it via rho_b * F_STAR10 * df_coll/dt, so another attention is given there. For the Eulerian source models, we multiply the emissivity fields by (1+delta), except for N_ion, and we also do mean fixing, just like it was used to be done until now. In addition, I optimized the allocation of the fields of HaloBox: halo_sfr and halo_sfr_mini are required only if USE_TS_FLUCT=True, while whalo_sfr is used only for the Lagrangian source models (for now). n_ion is allocated for all source models, even though it is used only for the Lagrangian source models, as this might change in the future.
Two bugs were fixed in RadiationFields.c for CONST-ION-EFF. The first one, and the more major one, is that mean_sfr_zpp was always computed via EvaluateSFRD. This is wrong because EvaluateSFRD contains a mass-dependent term for star formation efficiency and also exponential damping, unlike the collapsed fraction which has a sharp cutoff. On top of that, the local collapsed fraction is computed under the assumption of the PS mass function, but EvaluateSFRD assumes the input HMF, which is usually Sheth-Tormen. The second bug was that the mean fixing did something like [df_coll/dz](x) *= f_coll_global / <f_coll(x)>, instead of [df_coll/dz](x) *= df_coll_dz_global / <[df_coll/dz](x)>. The latter of course is more correct as it fixes the mean of the fluctuating quantity, which is the derivative of the collapsed fraction.
Changed our integration tests that work with CONST-ION-EFF to work also with HMF=PS. The reason for that is that in CONST-ION-EFF the fluctuating collpased fraction is taken from extended PS theory, so it is the more appropriate halo mass function to consider when working with CONST-ION-EFF. It also silences a warning that we have about working with CONST-ION-EFF and HMF which is different than PS.
Introduced a temporary USE_NEW_CODE boolean flag. This flag was helpful in testing and comparing the old E-INTEGRAL source model vs the new E-INTEGRAL source model. After fixing two bugs in HaloBox.c, the global spin temperature from run_global_evolution agree very well, with and without mini-halos. As for the global spin temperature from run_lightcone, divergence at low redshifts can be seen. I confirmed that this divergence is real and not a consequence of an uknown bug by comparing between L-INTEGRAL and the new E-INTEGRAL, but forced to still use the lagrangian density field: in that configuration, both source models had the same input density field, as well as the same algorithm for computing the radiation fields, and indeed the same output (global spin temperature from run_lightcone) was obtained. This means that the divergence at low redshifts in the global spin temperature between the old and new E-INTEGRAL source models (from run_lightcone) is a genuine real effect, and not a bug. The reason for that divergence is because filtering and computing the emissivity fields given the density field are two non-linear operations, and as such, they do not commute. In addition to confirming the correctness of the new E-INTEGRAL source model, all the tests in the test suite pass.
A few commits ago I thought I fixed a bug in CONST-ION-EFF in RadiationFields.c since I thought that mean_sfr_zpp is determined by the mass-dependent version of EvaluateSFRD, namely an integral over the user's mass function (which is not necessarily PS), with mass-dependent star formation efficiency term, and an exponential cutoff, all seem not to match the call per cell in calculate_sfrd_from_grid. However, I didn't notice that EvaluateSFRD has a branch for the mass-independent source model, where Fcoll_General is calculated (or an interpolation of a table that was initialized with this function is performed). This is actually okay, since Fcoll_General executes the correct integral, but with the user's halo mass function, so the mean fixing corrects the mean by using the correct user's halo mass function (rather than PS). I therefore restored the code, the test, and the data files for CONST-ION-EFF from a few commits ago. Having said that, there is still a subtle bug of computing avg_fix_term for the plain collpased fraction, while that avg_fix_term is used in correcting the derivative of the collapsed fraction. This still needs to be fixed and will be done in the next commit.
As I said in my last commit, there was still a bug of computing the mean ratio for the collapsed fraction for SOURCE_MODE=CONST-ION-EFF, while this ratio was applied to fix the derivative of the collapsed fraction. I fixed that in this commit. Since I use in get_uhmf_averages the numerical derivative of Fcoll_General, and the latter function doesn't need sigma_min, I could remove the integral_cond input and its associated lines that I added a few commits ago. As this fix changes the spin temperature output, I had to also update the datafile for the ts_nomdz configuration. I can tell that before updating the datafile, I got an AssertionError due to a mismatch in the global signal which was as high as 28%, though I think it could have been higher because we stop the test only at z=18. This shows anyway the impact that this bug had on the output of the simulation for CONST-ION-EFF.
Removed the temporary USE_NEW_CODE flag, by assuming it is True all the time. This allowed to clear many lines from the old code in RadiationFields.c, as well as in SpinTemperatureBox.c. As a consequence of this change, we now always compute HaloBox, and we also always compute RadiationFields.c (if USE_TS_FLUCT=True). In addition, we can now use LYA_MULTIPLE_SCATTERING also for the Eulerian source models. Due to adopting USE_NEW_CODE=True, the spin temperature for the Eulerian source models have changed, it is now higher (mostly at low redshifts when the density field becomes highly non-linear) due to the Jensen inequality, as we now compute first the emissivity fields (convex functions of the density field) and then filter, unlike the opposite ordering in the old code of the Eulerian source models. I therefore reproduced the relevant data files (that work with Eulerian source models and USE_TS_FLUCT=True), in order to pass the integration tests. Also minor adjustments had to be done in some test modules that test the python code. After all these changes, all the tests now pass.
Included a new rst file in the docs named Radiation_Fields.rst. That file describes the history of the radiation fields in 21cmFAST, including the changes in this PR
Until this PR, MINIMIZE_MEMORY also affected how the computation of the spin temperature and the radiation fields is performed. After this PR, it's no longer the case. I therefore removed test_minimize_memory.py, as this test module was originally included to test the code of MINIMIZE_MEMORY in SpinTemperatureBox.c. I also removed USE_TS_FLUCT from the minimize_mem configuration and updated the corresponding datafile from the same reason. Finally, I noticed that mean_free_path is evaluated only if we do recombination calculations, so I optimized its allocation in outputs.py
Reviewer's GuideUnifies emissivity and radiation field computation for Eulerian and Lagrangian source models, reworking HaloBox/RadiationFields/TsBox pipelines, fixing CONST-ION-EFF SFRD consistency, and simplifying MINIMIZE_MEMORY behavior and Python wrappers/docs/tests accordingly. Sequence diagram for unified emissivity/radiation pipeline in spin temperature computationsequenceDiagram
actor User
participant drivers_single_field as compute_spin_temperature
participant HaloBox as ComputeHaloBox
participant RadiationFields as UpdateRadiationFields
participant TsBox as ComputeTsBox
participant PerturbedField
participant InitialConditions
User->>drivers_single_field: compute_spin_temperature(inputs, redshift, ...)
drivers_single_field->>PerturbedField: load_all()
drivers_single_field->>InitialConditions: load_all()
drivers_single_field->>HaloBox: ComputeHaloBox(redshift, InitialConditions, PerturbedField, HaloCatalog, previous_spin_temp, previous_ionize_box, HaloBox)
note over HaloBox,PerturbedField: Emissivity fields (halo_sfr, halo_xray, n_ion, ...) computed
drivers_single_field->>RadiationFields: compute_radiation_fields(redshift, hboxes, PerturbedField, InitialConditions, ...)
RadiationFields->>RadiationFields: UpdateRadiationFields(redshift, HaloBox, R_inner, R_outer, R_ct, R_star, mode, cleanup, RadiationFields)
note over RadiationFields: Radiation fields (xray_ionization_rate, lya_flux_*, lyw_flux, ...) updated from emissivities
drivers_single_field->>TsBox: ComputeTsBox(redshift, prev_redshift, perturbed_field_redshift, TsBox, PerturbedField, IonizedBox, RadiationFields, cleanup)
TsBox->>TsBox: use RadiationFields.Q_HI and radiation grids for spin temperature
TsBox-->>User: TsBox result
Flow diagram for unified emissivity/radiation computation across Eulerian and Lagrangian source modelsflowchart LR
subgraph SourceModelSelection
A[source_model]
A -->|source_model_uses_eulerian_grids| B[Eulerian branch]
A -->|source_model_uses_lagrangian_grids| C[Lagrangian branch]
end
subgraph EmissivityComputation
B --> B1[ComputeHaloBox redshift InitialConditions PerturbedField HaloCatalog etc]
C --> C1[ComputeHaloBox redshift InitialConditions PerturbedField HaloCatalog etc]
B1 --> D[HaloBox emissivity fields halo_sfr halo_xray n_ion etc]
C1 --> D
end
subgraph RadiationFieldsComputation
D --> E[UpdateRadiationFields redshift HaloBox R_inner R_outer R_ct R_star etc]
E --> F[RadiationFields grids xray_ionization_rate xray_lya_flux lyw_flux lya_flux_*]
end
subgraph SpinTemperatureEvaluation
F --> G[ComputeTsBox redshift prev_redshift perturbed_field_redshift TsBox PerturbedField IonizedBox RadiationFields etc]
G --> H[TsBox spin temperature Q_HI etc]
end
File-Level Changes
Possibly linked issues
Tips and commandsInteracting with Sourcery
Customizing Your ExperienceAccess your dashboard to:
Getting Help
|
jordanflitter
added
context: python wrapper
![sourcery-ai[bot]](https://avatars.githubusercontent.com/in/48477?s=60&v=4){kind=small}
There was a problem hiding this comment.
Hey - I've found 13 issues, and left some high level feedback:
- The prefactor setup for SFRD/Nion/X-ray (e.g., in
get_uhmf_averages,get_cell_integrals, andmove_grid_galprops) is now duplicated and quite intricate; consider centralising this logic in a shared helper to reduce the risk of the Eulerian/Lagrangian and mass‑dependent/mass‑independent branches drifting out of sync. - In
move_grid_galprops, the strictvol_ratio_out != 1.0check for Eulerian source models may trip on floating-point round-off; using a small tolerance or checking integer volumes directly would make this safety check more robust.
Prompt for AI Agents
Please address the comments from this code review:
## Overall Comments
- The prefactor setup for SFRD/Nion/X-ray (e.g., in `get_uhmf_averages`, `get_cell_integrals`, and `move_grid_galprops`) is now duplicated and quite intricate; consider centralising this logic in a shared helper to reduce the risk of the Eulerian/Lagrangian and mass‑dependent/mass‑independent branches drifting out of sync.
- In `move_grid_galprops`, the strict `vol_ratio_out != 1.0` check for Eulerian source models may trip on floating-point round-off; using a small tolerance or checking integer volumes directly would make this safety check more robust.
## Individual Comments
### Comment 1
<location path="docs/Radiation_Fields.rst" line_range="29" />
<code_context>
+ionized electron of energy :math:`E_e = h_{\rm P}\nu - E^i_{\rm th}` that is deposited as heat in the IGM (a function of :math:`E_e` and
+the local ionization fraction), :math:`f_i` is the number fraction of species :math:`i` in the IGM, and :math:`\sigma_i(\nu)` is the
+cross-section for photo-ionization of a particle from species :math:`i` due to the interaction with an X-ray photon of frequency :math:`\nu`.
+To understand the above experssion, let us consider the following components:
+
+* :math:`f_i \sigma_i(\nu) J_X(\nu)` is the production rate of electrons (with energy :math:`E_e`) per gas particle per
</code_context>
<issue_to_address>
**issue (typo):** Correct typo in 'experssion'.
Spelling: use 'expression' instead of 'experssion'.
```suggestion
To understand the above expression, let us consider the following components:
```
</issue_to_address>
### Comment 2
<location path="docs/Radiation_Fields.rst" line_range="47-48" />
<code_context>
+
+where :math:`k_B` is Boltzmann's constant and :math:`x_e` is the local ionization fraction. The X-ray heating rate is thus a physical
+quantity that can be represented as a 3D box at a given redshift, and is thus considered a "radiation field" in ``21cmFAST``.
+Similarly, the photionization rate due to X-ray radiation, as well as the contribution to the Lyman-alpha flux from X-ray excitation of HI
+atoms, are also considered "radiation fields".
+
</code_context>
<issue_to_address>
**issue (typo):** Fix spelling of 'photionization'.
Use 'photoionization' instead of 'photionization' in this sentence.
```suggestion
Similarly, the photoionization rate due to X-ray radiation, as well as the contribution to the Lyman-alpha flux from X-ray excitation of HI
atoms, are also considered "radiation fields".
```
</issue_to_address>
### Comment 3
<location path="docs/Radiation_Fields.rst" line_range="82" />
<code_context>
+traditional radiative transfer codes.
+
+The Lyman alpha flux from continuum photons (photons that were emitted with frequency below Lyman beta and redshifted into Lyman alpha) and
+injected photons (photons that were emitted with frequency above Lyman beta and went atomic cascades once they had redshifted into a Lyman
+resonance) is computed in a similar manner,
+
</code_context>
<issue_to_address>
**suggestion (typo):** Improve grammar in 'went atomic cascades'.
Please rephrase 'went atomic cascades' to 'went through atomic cascades' or 'underwent atomic cascades' to correct the grammar.
Suggested implementation:
```
injected photons (photons that were emitted with frequency above Lyman beta and underwent atomic cascades once they had redshifted into a Lyman
resonance) is computed in a similar manner,
```
If this sentence is wrapped differently in your actual file (e.g., different line breaks), adjust the SEARCH block to match the exact existing text so the replacement applies correctly.
</issue_to_address>
### Comment 4
<location path="docs/Radiation_Fields.rst" line_range="128" />
<code_context>
+
+ \epsilon_\alpha(z, \delta) = \bar\n_b(z=0) (1+\delta) N_{\gamma /{\rm b}} f_* \frac{d}{dt}f_{\rm coll}(z, M_{\rm turn} ; \delta, M_{\rm cond}),
+
+where :math:`\bar\n_b` is the mean baryon number density and :math:`N_{\gamma /{\rm b}}` is the number of stellar photons per baryon
+(the latter is a free parameter).
+
</code_context>
<issue_to_address>
**issue (typo):** Fix LaTeX command `\bar\n_b`.
`\n` is not a valid LaTeX command, so this won’t render correctly. Use `\bar n_b` within the math environment instead.
```suggestion
\epsilon_\alpha(z, \delta) = \bar n_b(z=0) (1+\delta) N_{\gamma /{\rm b}} f_* \frac{d}{dt}f_{\rm coll}(z, M_{\rm turn} ; \delta, M_{\rm cond}),
```
</issue_to_address>
### Comment 5
<location path="docs/Radiation_Fields.rst" line_range="144" />
<code_context>
+
+where :math:`\bar\rho_m` is the mean matter density and :math:`dn/dM_h(z, M_h ; \delta_{\rm cond}, M_{\rm cond})` is the conditional halo
+mass function in a region of size :math:`R(M_{\rm cond}) = [3M_{\rm cond}/(4\pi \bar\rho_m)]^{1/3}` and overdensity :math:`\delta_{\rm cond}`.
+For the calculation of the conditionl collapsed fraction, ``21cmFAST`` v1.0.0 assumed that the conditional halo mass function was given
+by the Sheth-Tormen conditional mass function. Under this assumption, the conditional collapsed fraction has an analytical form,
+
</code_context>
<issue_to_address>
**issue (typo):** Correct spelling of 'conditionl'.
Replace 'conditionl' with 'conditional' in this sentence.
```suggestion
For the calculation of the conditional collapsed fraction, ``21cmFAST`` v1.0.0 assumed that the conditional halo mass function was given
```
</issue_to_address>
### Comment 6
<location path="docs/Radiation_Fields.rst" line_range="157" />
<code_context>
+radius :math:`R(M) = [3M/(4\pi \bar\rho_m)]^{1/3}`.
+
+While ``21cmFAST`` v1.0.0 used the analytical result for the conditional collapsed fraction from the extended Press-Schechter formalism
+in order to assess the fluctuations in the emissivity fields, the code also normalized the mean collapsed fraction in the box to match
+the global collapsed fraction, as given by solving the above intergral numerically with the user's selected halo mass function
+(which by default was the Sheth-Tormen mass function).
+
</code_context>
<issue_to_address>
**issue (typo):** Correct spelling of 'intergral'.
'intergral' should be corrected to 'integral' in that sentence.
```suggestion
the global collapsed fraction, as given by solving the above integral numerically with the user's selected halo mass function
```
</issue_to_address>
### Comment 7
<location path="docs/Radiation_Fields.rst" line_range="165" />
<code_context>
+
+Given the above definitions for the emissivity fields, the effective emissivity fields are achieved by interpolating in time/redshift and
+filtering in space. Since ``21cmFAST`` v1.0.0 did not have the emissivity fields on a grid, the interpolation and filtering was done in
+the following way. Firstly, given they Eulerian density field :math:`\delta_{\rm E}(z, \mathbf{x})` at redshift :math:`z` (the current
+snapshot's redshift), the code evaluated the density field at the integrated redshift :math:`z'` by scaling with the linear growth factor,
+namely :math:`\delta_{\rm E}(z', \mathbf{x}) \approx \delta_{\rm E}(z, \mathbf{x}) D(z')/D(z)`. Then, instead of filtering the emissivity
</code_context>
<issue_to_address>
**issue (typo):** Fix grammar in 'given they Eulerian density field'.
Replace 'they' with 'the' so the sentence reads 'given the Eulerian density field'.
```suggestion
the following way. Firstly, given the Eulerian density field :math:`\delta_{\rm E}(z, \mathbf{x})` at redshift :math:`z` (the current
```
</issue_to_address>
### Comment 8
<location path="docs/Radiation_Fields.rst" line_range="254" />
<code_context>
+where :math:`\nu_\alpha` is the Lyman-alpha frequency. Several notes:
+
+* Note that unlike previous radiation fields, the LW flux contains units of energy, owned by the factor of :math:`h_{\rm P}\nu_\alpha`.
+* The contribution to the LW flux came from both atomic cooling galaxies (ACGs) and MCGs, as the SFRD in both populations is modedled a bit
+differently, mostly in the modeling of the star formation efficiency, the turnover mass and the duty fraction (see more details on the
+latter two in :doc:`M_TURN`). Likewise, the contribution to the X-ray and Lyman-alpha fluxes also came from both ACGs and MCGs.
</code_context>
<issue_to_address>
**issue (typo):** Correct spelling of 'modedled'.
Replace "modedled" with "modeled" in this sentence.
```suggestion
* The contribution to the LW flux came from both atomic cooling galaxies (ACGs) and MCGs, as the SFRD in both populations is modeled a bit
```
</issue_to_address>
### Comment 9
<location path="docs/Radiation_Fields.rst" line_range="258" />
<code_context>
+differently, mostly in the modeling of the star formation efficiency, the turnover mass and the duty fraction (see more details on the
+latter two in :doc:`M_TURN`). Likewise, the contribution to the X-ray and Lyman-alpha fluxes also came from both ACGs and MCGs.
+* It was assumed that ACGs contained only popII stars, while MCGs contained only popIII stars. Therefore,
+the shape of the SEDs in ACGs and MCGs were also different in the evaluation of Lyman-alpha and LW fluxes, as well their amplitudes
+(controlled by the free parameter :math:`N_{\gamma /{\rm b}}`, which was now split into two free parameters). However, for the evaluation
+of the X-ray flux, the SEDs in both ACGs and MCGs were assumed to be the same, and were modeled as in previous versions, namely
</code_context>
<issue_to_address>
**suggestion (typo):** Improve phrasing 'as well their amplitudes'.
Change the phrase to "as well as their amplitudes" for correct grammar.
```suggestion
the shape of the SEDs in ACGs and MCGs were also different in the evaluation of Lyman-alpha and LW fluxes, as well as their amplitudes
```
</issue_to_address>
### Comment 10
<location path="docs/Radiation_Fields.rst" line_range="283-286" />
<code_context>
+* ``"E-INTEGRAL"``:
+ This option was equivalent to setting ``USE_MASS_DEPNDENT_ZETA`` to True in previous versions.
+
+* ``"L-INTEGRAL"``:
+ This is a new source model that was introduced in v4.0.0. It was similar to ``"E-INTEGRAL"``, but the emissivity
+ fields were computed on the Lagrangian density grid, and then mapped to the Eulerian grid (see details below).
+
+* ``"DEXM-ESF"``:
+ This is a new source model that was introduced in v4.0.0, similar to ``"L-INTEGRAL"``, but the emissivity fields recieved
+ contributions from resolved discrete halos that were found via the excursion-set formalism (ESF), as was implemented by
+ the ``DEXM`` algorithm (see Mesinger & Furlanetto 2007, https://arxiv.org/pdf/0704.0946).
</code_context>
<issue_to_address>
**issue (typo):** Correct spelling of 'recieved'.
In this description, change 'recieved' to 'received'.
```suggestion
* ``"DEXM-ESF"``:
This is a new source model that was introduced in v4.0.0, similar to ``"L-INTEGRAL"``, but the emissivity fields received
contributions from resolved discrete halos that were found via the excursion-set formalism (ESF), as was implemented by
the ``DEXM`` algorithm (see Mesinger & Furlanetto 2007, https://arxiv.org/pdf/0704.0946).
```
</issue_to_address>
### Comment 11
<location path="docs/Radiation_Fields.rst" line_range="308-309" />
<code_context>
+
+ \epsilon(z, \mathbf{x}) \to \epsilon^{\rm eff}(z', \mathbf{x})=\epsilon[z', \delta_{\rm E}^R(z,\mathbf{x}) D(z')/D(z)].
+
+In contrast, for the three "Lagrangian" source models, the emissivity fields were computed on the Lagrangian density grid, and were
+then advected to the Eulerian grid by an algorithm similar to 2LPT. Moroever, the emissivity fields in the "Lagrangian" source models
+were computed as an interemediate output of the code (stored in the ``HaloBox`` class). This enabled to interpolate the emissivity fields
+at redshift :math:`z'`. Then, for each spherical integrated shell with innter radius :math:`R_{\rm i}` and outer radius :math:`R_{\rm o}`,
</code_context>
<issue_to_address>
**issue (typo):** Correct spelling of 'Moroever'.
Please correct the spelling to 'Moreover'.
Suggested implementation:
```
then advected to the Eulerian grid by an algorithm similar to 2LPT. Moreover, the emissivity fields in the "Lagrangian" source models
```
While not part of your specific request, you may also want to correct other typos in the same paragraph (e.g., "interemediate" → "intermediate", "innter" → "inner") for clarity and professionalism.
</issue_to_address>
### Comment 12
<location path="docs/Radiation_Fields.rst" line_range="273-277" />
<code_context>
+emissivity fields and the radiation fields was changed for several new configurations that were introduced in that version. This was
+controlled by the new enum parameter ``SOURCE_MODEL``, which had five options:
+
+* ``"CONST-ION-EFF"``:
+ This option was equivalent to setting ``USE_MASS_DEPNDENT_ZETA`` to False in previous versions.
+
+* ``"E-INTEGRAL"``:
</code_context>
<issue_to_address>
**issue (typo):** Fix typo in flag name `USE_MASS_DEPNDENT_ZETA`.
Earlier the flag is spelled `USE_MASS_DEPENDENT_ZETA`, but here it appears as `USE_MASS_DEPNDENT_ZETA` (missing `E`). Please correct this for consistency and to avoid confusion when users reference the flag name.
```suggestion
* ``"CONST-ION-EFF"``:
This option was equivalent to setting ``USE_MASS_DEPENDENT_ZETA`` to False in previous versions.
* ``"E-INTEGRAL"``:
This option was equivalent to setting ``USE_MASS_DEPENDENT_ZETA`` to True in previous versions.
```
</issue_to_address>
### Comment 13
<location path="docs/Radiation_Fields.rst" line_range="341" />
<code_context>
+that in the former, the emissivity fields were computed on the Eulerian density grid, while in the latter, the emissivity fields were
+computed on the Lagrangian density grid and then advected to the Eulerian grid.
+
+Furthermore, in v4.3.0, the effect of Lyman-alpha multiple scattering could have been applied via a new flag called ``LYA_MULTIPLE_SCATTERING``,
+regardless if the simulated source model was "Eulerian" or "Lagrangian". When this flag was set to True, the effective emissivity for
+Lyman-alpha photons was achieved by filtering the interpolated emissivity field with a generalization of the sphericall shell filter
+(see more details in Flitter, Munoz and Mesinger 2026, https://arxiv.org/pdf/2601.14360).
+
</code_context>
<issue_to_address>
**issue (typo):** Correct spelling of 'sphericall'.
```suggestion
Lyman-alpha photons was achieved by filtering the interpolated emissivity field with a generalization of the spherical shell filter
```
</issue_to_address>Help me be more useful! Please click 👍 or 👎 on each comment and I'll use the feedback to improve your reviews.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## release-v4.3 #755 +/- ##
================================================
- Coverage 91.45% 91.42% -0.04%
================================================
Files 33 33
Lines 5245 5224 -21
Branches 921 914 -7
================================================
- Hits 4797 4776 -21
Misses 286 286
Partials 162 162 ☔ View full report in Codecov by Harness. |
Corrected some typos/rendering issues on the docs. In addition, I fixed a test that somehow managed to pass in certain OS (weird...)
Contributor
Author
Validation testsIn order to validate the new code, I introduced a temporary I started first by running two quick 
As can be seen, the residuals of the global spin temperature are almost sub-percent at all redshifts and Eulerian source models, and because all the integration tests for the Lagrangian source models have passed without modifying the datafiles, this means that there is no bug in the prefactors that multiply the integrals, or in the mean fixing. Interestingly, if I repeat the same comparison, but for the global spin temperature from 
|
Contributor
Author
Origin of the discrepancy at low redshifts (Jensen's inequality!)The main difference between the old and new codes for the Eulerian source models is flipping the order of operations:
Because the emissivity function is a convex function of the density field, and filtering can be viewed as a weighted average, we can apply the Jensen's inequality: for a convex function At sufficiently low redshifts, the WF coupling saturates, and the spin temperature follows the gas kinetic temperature. Because of the increased X-ray flux in the new code, X-ray heating is more efficient, which means that the gas heats up faster and the spin temperature is higher. We see exactly that when we compare the coeval boxes of the spin temperature between the old and new codes! Interestingly, it also seems that the spin temperature field is smoother in the new code. 
Because the spin temperature is higher (but not much higher compared to CMB temperature), the brightness temperature is also higher. Because the spin temperature in the new code is smoother, also the brightness temperature becomes smoother, so power in small scales (high 
|
steven-murray
left a comment
steven-murray
left a comment
Member
There was a problem hiding this comment.
Thanks @jordanflitter for the code and for the tests and plots. This is looking good to me, with the caveat of some comments.
Comment on lines
+2546
to
+2551
| # TODO: the last warning below (that we also test), is only reachable if max_integral_mass is larger than | ||
| # the upper mass_limits, but if we have halos then the latter is max_dexm which is exactly the same as max_integral_mass. | ||
| # Thus the conclusion is that this warning is reachable only when we use L-INTEGRAL. Therefore, I am keeping the quick return | ||
| # if we don't use lagrangian_source_grid (instead of a quick return if we don't use has_discrete_halos), but I suggest to | ||
| # modify the code below to include also Eulerian source models, (and the comment above that says that "there are no problems | ||
| # if we are not using halos") |
Member
There was a problem hiding this comment.
Trying to follow this... should we not implement that in this PR then?
Contributor
Author
There was a problem hiding this comment.
I am not sure. I am bit confused about this function. We test that the last warning in this function is triggered in test_input_structs.py::test_halomass_ranges with SOURCE_MODEL="L-INTEGRAL", but I think it doesn't make sense because we are supposed to integrate over the entire mass-range in this source model, with no gaps. So I think this function/test should be revised, and I am not sure if it should be in this PR (because this check doesn't do anything with Eulerian source models, I just wonder if it's really relevant to L-INTEGRAL).
Contributor
Author
There was a problem hiding this comment.
I fixed check_halomass_range and test_halomass_ranges in the last commit.
|
|
||
| .. math:: | ||
|
|
||
| \epsilon_X(z, \delta) = \bar\rho_b(z=0) (1+\delta) \zeta_X f_* \frac{d}{dt}f_{\rm coll}(z, M_{\rm turn} ; \delta, M_{\rm cond}), |
Member
There was a problem hiding this comment.
For my own enrichment: this is saying that I generate X-ray photons at a rate proportional to the rate at which I form stars. I suppose this is empirically reasonable but it seems a bit strange: in my imagination if I had a galaxy that had stopped forming stars, it could still produce X-ray photons. Is there some empirical work that this assumption is based on?
Contributor
Author
There was a problem hiding this comment.
I completely agree with everything you said ("you hit the nail on the head!", as LLMs like to say).
Up to this day, we still assume that the emissivity fields are proportional to the SFRD. This is very clear from the modeling of the X-ray emissivity, which is proportional to L_X/SFR *SFRD, where L_X/SFR is our astro parameter (bad naming!), L_X/SFR ~ 10^40 erg/sec *yr/M_sun (it actually accounts for photons up to energy 2 keV). Anyway, the v2 paper where this modeling is first introduced is https://arxiv.org/pdf/1809.08995, where it seems like they took the empiric relation from https://arxiv.org/pdf/1702.00409. But I totally agree with you, I always thought that it's weird to assume that the star-formation timescale is the relevant timescale for the calculation of the emissivity (but maybe I'm missing something...). Anyway, this is the justification for the X-ray emissivity. As for the emissivity that is used for the Lya flux, I don't know what is the justification to assume that it is proportional to the SFRD.
And note that while you reached the correct conclusion, that part you were commenting about was regarding the modeling in v1. In that paper/version, there was no reference to the SFRD in the model, but if we compare the equations with v2 (which was the first 21cmFAST paper to use the term "SFRD" in the context of modeling the emissivity), then we can see that the SFRD in v1 was proportional the derivative of the collapsed fraction.
| where :math:`\bar\rho_m` is the mean matter density and :math:`dn/dM_h(z, M_h ; \delta_{\rm cond}, M_{\rm cond})` is the conditional halo | ||
| mass function in a region of size :math:`R(M_{\rm cond}) = [3M_{\rm cond}/(4\pi \bar\rho_m)]^{1/3}` and overdensity :math:`\delta_{\rm cond}`. | ||
| For the calculation of the conditional collapsed fraction, ``21cmFAST`` v1.0.0 assumed that the conditional halo mass function was given | ||
| by the Sheth-Tormen conditional mass function. Under this assumption, the conditional collapsed fraction has an analytical form, |
Member
There was a problem hiding this comment.
Is this the ST CHMF or the PS CHMF?
Contributor
Author
There was a problem hiding this comment.
So in theory, the conditional mass function that should have been used there is ST, but since the collapsed fraction only has analytic from with PS, then the PS mass function was used in computing the local conditional collapsed fraction. Then, as the text below explains, the ST mass function was used to normalize the collapsed fraction, such that the mean collapsed fraction in the box would match the theoretical global collapsed fraction (while assuming the ST mass function).
Comment on lines
+155
to
+158
| While ``21cmFAST`` v1.0.0 used the analytical result for the conditional collapsed fraction from the extended Press-Schechter formalism | ||
| in order to assess the fluctuations in the emissivity fields, the code also normalized the mean collapsed fraction in the box to match | ||
| the global collapsed fraction, as given by solving the above integral numerically with the user's selected halo mass function | ||
| (which by default was the Sheth-Tormen mass function). |
Member
There was a problem hiding this comment.
So, I think what you are saying here is that the CHMF equation here, if integrated over the PDF of delta, by construction gives the global Press-Schechter HMF. However, since the user is able to select a non-press-shechter global HMF, we still use this PS formula for the CHMF, but we apply a renormalization such that \int d delta p(delta) CHMF == Custom_HMF(sigma). Is that right?
Contributor
Author
There was a problem hiding this comment.
Yes, that's exactly correct ("you hit the nail on the head again!")
Today btw, we don't do mean fixing if we use ST mass function and Lagrangian source models. This is clear from the following line (in scaling_relations.c):
consts->fix_mean =
source_model_uses_eulerian_grids(matter_options_global->SOURCE_MODEL) ||
matter_options_global->HMF == HMF_WATSON || matter_options_global->HMF == HMF_WATSON_Z ||
matter_options_global->HMF == HMF_REED07 || matter_options_global->HMF == HMF_YUNG24;I don't like this logic and indeed we have an open issue about it, #669.
This merger includes cautious resolution of conflicts. All tests passed.
Overhaul of check_halomass_range. It now matches better with the C code. It works only if we have discrete halos in the simulation, otherwise everything is fine. The warning about minimum mass was removed as it wasn't reachable. We raise a warning if the largest halo in the simulation (derived from the box size) is not very big. Also, test_halomass_ranges was modified.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fix #668.
Before this PR, the Eulerian source models did not compute the emissivity fields (in
HaloBox) and filter them to get the radiation fields (inRadiationFields). Instead, the density field was filtered in order to compute the "effective emissivity field" for each shell, while the radiation fields were computed on the fly in the evaluation of the spin temperature. This PR changes that and unifies the prescription for computing the emissivity and radiation fields, regardless if the considered source models are Eulerian or Lagrangian. This change allowed to clean many redundant lines inRadiationFields.cthat were used for the "old" Eulerian source models, making the code much cleaner and more modular (but there's still some work to do...).In evaluating the spin temperature, the only difference between the two types of source models is that in the Eulerian source models we use the Eulerian density field to compute the emissivity fields on the Eulerian grid, while for the Lagrangian source models we use the Lagrangian field to compute the emissivity fields on the Lagrangian grid, and then we advect them to the Eulerian grid (there's also a multiplication by
1+deltaof the emissivity fields in the Eulerian source models, I am still not sure if it's really required...). There are however still differences in how the code computes the reionization field for the two types of source models, but I left this work for a future PR.Because this PR changes the output (i.e. realizations of the spin/brightness temperature field) for the Eulerian source models, I have made several tests to confirm that I did not introduce a new bug, see my next comment below.
Other more minor changes that were made:
CONST-ION-EFF. In this source model, it is assumed that the SFRD/emissivity fields follow the time derivative of the collapsed fraction. However, the mean fixing ratioavg_fix_termwas computed for the plain collapsed fraction, but was applied on the derivative of the collapsed fraction. Compared to what we do in the other source models, it is more consistent to computeavg_fix_termfor the derivative of the collapsed fraction (which is proportional to the SFRD). That was fixed in this PR.Radiation_Fields.rstwas added to the docs. This file describes the history of how emissivity fields and radiation fields were computed in 21cmFAST, including the changes that were made in v4.3 and were part of this PR.MINIMIZE_MEMORYaffected how the radiation fields were computed in the Eulerian source models. After this PR, its effect on the code has been reduced significantly, which allowed to delete/optimize some tests.Summary by Sourcery
Unify emissivity generation and radiation-field filtering across Eulerian and Lagrangian source models.
New Features:
Bug Fixes:
Enhancements:
Documentation:
Tests: