ring sturcture shows up in large-volume simulations

[qyx268]

Describe the bug:

When we have large-volume simulations, we get ring structures on scales of Gpc, in particularly in the photoionization rate box. I'm thinking if this is due to the

To Reproduce: You just need a large-volume simulation... You can reproduce with this script. run.txt

Expected behavior:

Details:

• OS: Linux nefertem.sns.it 2.6.32-504.hz100.el6.x86_64 #1 SMP Fri Oct 31 20:08:27 CDT 2014 x86_64 x86_64 x86_64 GNU/Linux
• Python Version: 3.6.5

Additional context

Has anyone had large-volume simulations performed with the old 21cmfast or 21CMMC driver? It will be helpful to see if this is a new thing in v3 or it has always been there but missed.

[qyx268]

@andreimesinger any further thoughts here? apart from R_BUBBLE_MAX, which is set to 50 atm when INHOMO_RECO is on, I don't even know where to begin... I'm thinking of setting it to 200 and see if the rings disappear or look differently...

[andreimesinger]

this one is very strange… i wouldn’t know where to begin.. your test sounds good, though i cannot see how that could be the issue… it looks like some sampling problem.

i would need to look at the code to be more useful here. maybe monday i can see, if it hasn’t been fixed by then

On 05.12.2020, at 23:54, Yuxiang Qin notifications@github.com wrote:

@andreimesinger https://github.com/andreimesinger any further thoughts here? apart from R_BUBBLE_MAX, which is set to 50 atm when INHOMO_RECO is on, I don't even know where to begin... I'm thinking of setting it to 200 and see if the rings disappear or look differently...

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-738413092, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAUAWKHCXM2IGAI2NS3STAJKBANCNFSM4UMM4GLA.

[BradGreig]

Hmm, this is going to be problematic to discover given the size of the simulation. I've had a brief mess around but haven't come up with anything yet. I've tried a L = 500 Mpc, HII_DIM = 400 box (evolving the density linearly to try and minimise compute time) but that didn't show up anything.

I have seen ring structure before, however, that was around each source not centred on one location like this is. Switching the HII_FILTER from k-space to real-space resolved that, however, I don't think it'll help in this instance.

If you have stored all the intermediate boxes maybe look at a single co-eval box and see if its visible in any slice. Maybe you'll be able to identify the first redshift where it appears and then be able to only run a single co-eval box to establish what is going on.

[steven-murray]

@qyx268 is this one run the only time you've seen it occur? Or does it occur in every run?

[qyx268]

@steven-murray no, this is the only large-volume run i have ever had...

@BradGreig yes, i checked the co-eval boxes. they are visible.

I think they are there even at higher redshifts. But they dont show up, because Gamma is low at high-z. But if you flip between these redshifts, you can sort of tell they are present at high z too....

[andreimesinger]

i don’t see it in high redshifts. only post eor where we don’t trust the calculation so much. we should look into this, but the velocity bug is more urgent. imo

On 04.12.2020., at 18:17, Yuxiang Qin notifications@github.com wrote:

 @steven-murray no, this is the only large-volume run i have ever had...

@BradGreig yes, i checked the co-eval boxes. they are visible.

I think they are there even at higher redshifts. But they dont show up, because Gamma is low at high-z. But if you flip between these redshifts, you can sort of tell they are present at high z too....

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[steven-murray]

@qyx268 Might it be a good idea to run a few with different options switched on and off to hone in on the likely cause?

[qyx268]

yeah, the thing is we don't really know what is likely causing this... The only thing I am betting on is actually in GlobalParams. If this doesn't give me any hint, I will do a grid search I guess...

[steven-murray]

yeah, that could be an onerous process. Is it just in the photoionization rate? Or are there other quantities that show it?

[qyx268]

looks like it's only in Gamma. This plot shows from left to right, Gamma, z_re, kinetic T, log dNrec, log Fcoll.

[andreimesinger]

what was the max mean free path used here? there is a clear imprint of the Rmax here, which comes from using the code beyond the regime in which it was intended…

does the ring like thing change if you pick different slices?

On 7 Dec 2020, at 19:05, Yuxiang Qin notifications@github.com wrote:

looks like it's only in Gamma. This plot shows from left to right, Gamma, z_re, kinetic T, log dNrec, log Fcoll.

https://user-images.githubusercontent.com/15994713/101387698-1bb9c300-38bf-11eb-9533-6a2f023e51c6.png — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740084560, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UASFX2BDPNXIZNTIP4TSTUKONANCNFSM4UMM4GLA.

[qyx268]

it changes

Rmax = 50 when inhomo_reco is enabled

[andreimesinger]

hmm… the fact that the center is always at the same place is really strange..

On 7 Dec 2020, at 19:18, Yuxiang Qin notifications@github.com wrote:

it changes https://user-images.githubusercontent.com/15994713/101388858-d0081900-38c0-11eb-9bf7-bc8705325019.png Rmax = 50 when inhomo_reco is enabled

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740091432, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAWRGSLCHP67A7APYBDSTUL6BANCNFSM4UMM4GLA.

[qyx268]

i assume that's the density peak

[andreimesinger]

yeah, but if you go far away with your slices you will lose that peak.. did you plot distant slices or neighboring ones?

On 7 Dec 2020, at 19:20, Yuxiang Qin notifications@github.com wrote:

i assume that's the density peak

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740092609, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAU57TUOE7HJI7NX6LDSTUMGFANCNFSM4UMM4GLA.

[qyx268]

Distant..

On 7 Dec 2020, at 7:21 pm, andreimesinger notifications@github.com wrote:

yeah, but if you go far away with your slices you will lose that peak.. did you plot distant slices or neighboring ones?

On 7 Dec 2020, at 19:20, Yuxiang Qin notifications@github.com wrote:

i assume that's the density peak

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740092609, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAU57TUOE7HJI7NX6LDSTUMGFANCNFSM4UMM4GLA.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740093348, or unsubscribe https://github.com/notifications/unsubscribe-auth/AD2A6WP3MEKS3Z5CTPVAVLLSTUMLBANCNFSM4UMM4GLA.

[steven-murray]

This is all weird. I can definitely see the ring structure (less pronounced) in the other fields.

[BradGreig]

@qyx268 it might be worth sorting the density field to see the high density peaks. To see if one or a few high values are near the central region of this structure.

In terms of debugging, it might be worth extracting a smaller cube centred around this area and running the code on the sub-volume. This will likely require hacking the code, but might assist debugging it faster by significantly reducing the simulation volume.

I still recommend switching the HII_FILTER to see if this goes away.

[qyx268]

Those are good suggestions. I will take a look at the HII_FILTER. Here are results with Rmax=200...

[qyx268]

                            // check if fully ionized!
                            if ( (f_coll * ION_EFF_FACTOR + f_coll_MINI * ION_EFF_FACTOR_MINI> (1. - xHII_from_xrays)*(1.0+rec)) ){ //IONIZED!!
                                // if this is the first crossing of the ionization barrier for this cell (largest R), record the gamma
                                // this assumes photon-starved growth of HII regions...  breaks down post EoR
                                if (flag_options->INHOMO_RECO && (box->xH_box[HII_R_INDEX(x,y,z)] > FRACT_FLOAT_ERR) ){
                                    box->Gamma12_box[HII_R_INDEX(x,y,z)] = Gamma_R_prefactor * f_coll + Gamma_R_prefactor_MINI * f_coll_MINI;
                                }

                                // keep track of the first time this cell is ionized (earliest time)
                                if (previous_ionize_box->z_re_box[HII_R_INDEX(x,y,z)] < 0){
                                    box->z_re_box[HII_R_INDEX(x,y,z)] = redshift;
                                } else{
                                    box->z_re_box[HII_R_INDEX(x,y,z)] = previous_ionize_box->z_re_box[HII_R_INDEX(x,y,z)];
                                }

                                // FLAG CELL(S) AS IONIZED
                                if (global_params.FIND_BUBBLE_ALGORITHM == 2) // center method
                                    box->xH_box[HII_R_INDEX(x,y,z)] = 0;
                                if (global_params.FIND_BUBBLE_ALGORITHM == 1) // sphere method
                                    update_in_sphere(box->xH_box, user_params->HII_DIM, R/(user_params->BOX_LEN), \
                                                     x/(user_params->HII_DIM+0.0), y/(user_params->HII_DIM+0.0), z/(user_params->HII_DIM+0.0));
                            } // end ionized

This is how we assign values to Gamma12_box. I think at lower redshifts, the first if is True even at the largest filtering radius. Therefore, the ring structure depends on Rmax.

[andreimesinger]

i think you need to extend your color map for gamma… but you still get the “feature” that there is a “center” to the maps… i don’t know where this would come from… i think the code needs to be checked; doing these diagnostics is pretty random/inefficient without understanding what the code is doing… i can take a look in a bit

On 05.12.2020, at 15:19, Yuxiang Qin notifications@github.com wrote:

Those are good suggestions. I will take a look at the HII_FILTER. Here are results with Rmax=200... https://user-images.githubusercontent.com/15994713/101494197-4d866480-3967-11eb-88f6-e42631eb50f2.png https://user-images.githubusercontent.com/15994713/101494393-91796980-3967-11eb-99d7-85bbdbf6c3b0.png — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740646627, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UARLIIFYEM6VO5R46PDSTYYVPANCNFSM4UMM4GLA.

[qyx268]

the color bar is extended. otherwise, this is what u will see

[qyx268]

I'm not sure I follow the point about there is a center. We do filtering in 3D. The ring in the 2D maps are shells in 3D. It just says there is a density peak in 3D. What's wrong with this?

e.g. this is what happens if I slide on different plane

[andreimesinger]

btw, the Rmax is necessary in this calculation post EoR… because we originally were not intending to use the code post EoR, i took the value of 50 Mpc as a reasonable max at z~5.5-6. but post overlap, the code doesn’t make sense without having an Rmax as either a free parameter, or taking an observational estimate of mfp redshift evolution

On 05.12.2020, at 15:26, Yuxiang Qin notifications@github.com wrote:

the color bar is extended. otherwise, this is what u will see https://user-images.githubusercontent.com/15994713/101496155-b40c8200-3969-11eb-9072-1fbb72704f0f.png — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740651042, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAX2VROCZK6V2KN3U2DSTYZRRANCNFSM4UMM4GLA.

[andreimesinger]

why is there a center, even if does correspond to a density peak? if you filter the density field with some R_filter = 50Mpc, you just get a smoothed density field; you don’t get one sphere..

On 05.12.2020, at 15:31, Yuxiang Qin notifications@github.com wrote:

I'm not sure I follow the point about there is a center. We do filtering in 3D. The ring in the 2D maps are shells in 3D. It just says there is a density peak in 3D. What's wrong with this?

e.g. this is what happens if I slide on different plane

https://user-images.githubusercontent.com/15994713/101496644-4b71d500-396a-11eb-83c7-e03e6d23674e.png https://user-images.githubusercontent.com/15994713/101496659-4f9df280-396a-11eb-9949-622dbd8b4a24.png — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740653803, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UATMYNAIMN6ANKVVT3DSTY2DFANCNFSM4UMM4GLA.

[qyx268]

i'm not saying this is correct. I think this might be due to the filtering as brad said.

[qyx268]

looks up, some how, we update Gamma_12 using the sphere method... but i cant find it easily in the code

[andreimesinger]

yes, it is looks like it is something in the sampling, like that only some cells end up getting sampled or updated

On 05.12.2020, at 15:42, Yuxiang Qin notifications@github.com wrote:

looks up, some how, we update Gamma_12 using the sphere method... but i cant find it easily in the code

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740660301, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAS5TS34ZC56QL3UKC3STY3MBANCNFSM4UMM4GLA.

[qyx268]

speaking of using Rmax in post EoR, if this is a MFP argument, shouldnt Gamma be more or less a uniform background when the bubble overlapped?

[andreimesinger]

exactly… it should be something very similar to filtering the density field on the scale of the mfp (with some normalization of course)

On 05.12.2020, at 15:47, Yuxiang Qin notifications@github.com wrote:

shouldnt Gamma be more or less a uniform background when the bubble overlapped?

[andreimesinger]

looking at the code, i cannot see anything obviously wrong (at least not in my OC version). the only thing that comes to mind is that the fcoll interpolation tables might be inaccurate in some extreme regime, maybe producing non continuous outputs… because from the image, it looks like the ring is centered around a void and the rings come from concentric filterings around this void.

i think one needs to print out some same fcoll and gamma outputs from this if statement:

      // if this is the first crossing of the ionization barrier for this cell (largest R), record the gamma
      // this assumes photon-starved growth of HII regions...  breaks down post EoR
      if (INHOMO_RECO && (xH[HII_R_INDEX(x, y, z)] > FRACT_FLOAT_ERR) ){
    Gamma12[HII_R_INDEX(x, y, z)] = Gamma_R;

note the warning message i had put here when i wrote this… in the medium-long term, we should improve the physics here, so as to enable it to work better post EoR…

-a

On 8 Dec 2020, at 15:51, Andrei Albert Mesinger andrei.mesinger@sns.it wrote:

exactly… it should be something very similar to filtering the density field on the scale of the mfp (with some normalization of course)

On 05.12.2020, at 15:47, Yuxiang Qin <notifications@github.com mailto:notifications@github.com> wrote:

shouldnt Gamma be more or less a uniform background when the bubble overlapped?

[qyx268]

still haven't had a change to look into this. But the run @BradGreig suggested with HII_FILTER=0 has finished. Just wanted to post the result here if anyone has some insights...

This can be compared to https://github.com/21cmfast/21cmFAST/issues/194#issue-756476531

compared to https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740084560

compared to https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-740091432

[steven-murray]

What is that weird cross structure?!

[andreimesinger]

looks like it is aliasing around the same point (underdensity probably) that was causing the trouble with the other filter… i still suspect some of the fcoll interpolation tables are giving outlying results at these redshifts, driving these weird things

On 16 Dec 2020, at 17:56, Steven Murray notifications@github.com wrote:

What is that weird cross structure?!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-746629080, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAUGRJWQICB2ZKK5UDLSVDRDJANCNFSM4UMM4GLA.

[andreimesinger]

the cross is the center of the previous “wheel”… this is why i think it is a “hot pixel” or something

On 16 Dec 2020, at 17:56, Steven Murray notifications@github.com wrote:

What is that weird cross structure?!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-746629080, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UAUGRJWQICB2ZKK5UDLSVDRDJANCNFSM4UMM4GLA.

[steven-murray]

I agree, there's definitely something numerically strange going on with that pixel/small region.

On Wed, Dec 16, 2020 at 10:05 AM andreimesinger notifications@github.com wrote:

the cross is the center of the previous “wheel”… this is why i think it is a “hot pixel” or something

On 16 Dec 2020, at 17:56, Steven Murray notifications@github.com wrote:

What is that weird cross structure?!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub < https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-746629080>, or unsubscribe < https://github.com/notifications/unsubscribe-auth/ADH2UAUGRJWQICB2ZKK5UDLSVDRDJANCNFSM4UMM4GLA .

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-746642307, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJWRXW5MLSRJRIHNJPRFRLSVDSAPANCNFSM4UMM4GLA .

[qyx268]

some follow up on this... in case people have any further suggestions...

So this appears to exist only in high-res, >=1Gpc boxes.

e.g. this is Gamma slices at low-z in a 800Mpc simulation:

this is Gamma slices at low-z in a 1Gpc but low-res simulation:

The most important thing is that this also appears to be affecting the global eor history

I've also check simulation outputs from v2. Unless the seed is playing a role here, this bug only exists in v3... (p.s. the 1125Mpc run used a large Mturn)

[BradGreig]

Hey @qyx268, can you determine the min/max density for each filtering radius for any co-eval box that contains the ring-like structure. You should also consider outputting the cell location, to see if it is located near the centre of the ring (you roughly know the centre in the x-y plane). Or, extract a cuboid around roughly the centre of the ring (knowing the x-y coordinates, and look through the z direction) to search for extreme densities etc.

If you can determine this, then you might be able to track through this (or a few cells) and output various quantities. This way you might be able to establish what is going on, and where.

It might also be worth trying the exact same setup you have used except setting USE_MASS_DEPENDENT_ZETA = False. Setting this to False will remove the need for any interpolation tables, and use the error function for f_coll. This would highlight if any of the interpolation tables (for mass dependent zeta) are having issues.

[qyx268]

Thanks for these suggestions. I did some of then already. e.g. the ring center doesn't always match the location of the lowest density... also checked the min/max density before filtering (i.e. Perturb), 800Mpc which has no ring structures can span a larger density range than the 1Gpc run... I'll try to check the smoothed density, though I doubt exceeding the density range in the interpolation table is the case though...

I'm trying some switches at the moment, I will try USE_MASS_DEPENDENT_ZETA later. The RAM on nefertem is pretty full atm..

[BradGreig]

Hey @qyx268, just thinking about this some more. It might be worth outputting delta_R (smoothed) and f_coll for all smoothing radii.

Looking at your HII_FILTER = 0 case, you can see the aliasing in the xHI field. Thus, the source of this can only be coming from f_coll or delta_R. I suspect this aliasing will then be prevalent in either f_coll or delta_R (most likely f_coll). Where it appears (i.e. which R) would then allow you to narrow down the source of the problem.

[qyx268]

I did show f_coll at the last filtering radius. I'll look into this more.

looks like it's only in Gamma. This plot shows from left to right, Gamma, z_re, kinetic T, log dNrec, log Fcoll.

So far, I can confirm that with USE_MASS_DEPENDENT_ZETA off, we still see this.

I'm still not convinced this is due to existing codes that we just copy from v2. I wanted to see the result from PERTURB_ON_HIGH_RES=False

[BradGreig]

Interesting, the fact it appears in the USE_MASS_DEPENDENT_ZETA = False case means it is highly likely linked to the density field. Maybe try outputting the smoothed density field (see if it suddenly appears for the certain R).

Also, xHI looks very weird above. Maybe an issue with the lightcone generation?

[qyx268]

yeah, that's a good idea. I will check them.

the lightcone looks weird because I only have 7 coevolve boxes.

p21c.global_params.Z_HEAT_MAX = 15.0
p21c.global_params.ZPRIME_STEP_FACTOR = 1.2
redshift = 4.9

[qyx268]

some updates in case people have suggestions.

1. with USE_MASS_DEPENDENT_ZETA = False, USE_INTERPOLATION_TABLES=False, PERTURB_ON_HIGH_RES=False, USE_FFTW_WISDOM=False, the rings are still there... This does help shorten IonisationBox.c a lot. However, I still cant see anything wrong...
2. There is a small difference between v2 and v3 -- the size/structure of Fcoll, it is HII_TOT_FFT_NUM_PIXELS and using HII_R_FFT_INDEX in v2 while in v3, it's HII_TOT_NUM_PIXELS and HII_R_INDEX. However, this doesn't matter and having Fcoll in HII_R_INDEX system is actually better since then Fcoll can be reshaped into (HII_DIM,HII_DIM,HII_DIM).
3. Fcoll at different filtering radius looks okay...
4. Previously, I added a calculation of MFP because of the forest project. It's in the same if statement where we assign Gamma12

   if (flag_options->INHOMO_RECO && (box->xH_box[HII_R_INDEX(x,y,z)] > FRACT_FLOAT_ERR) ){
                                   box->Gamma12_box[HII_R_INDEX(x,y,z)] = Gamma_R_prefactor * f_coll;
                                   box->MFP_box[HII_R_INDEX(x,y,z)] = R;
                               }

   Then I found the ring structures also exist in MFP. (left to right: Gamma12_box, z_re_box, temp_kinetic_all_gas , dNrec_box, Fcoll[lagest R], and MFP_box).

I'm doing a run where rec is fixed at zero to see if something wrong in rec (i.e. N_rec_filtered)... which i really doubt though considering dNrec_box looks fine. let's see...

[qyx268]

fixing rec=0 does remove the rings...

[qyx268]

[BradGreig]

Hmmm, but I don't think it is the recombinations that is the actual problem. The structure might be appearing there, but the source must be elsewhere. For example, if it was in the recombinations, then it should have also appeared in the USE_MASS_DEPENDENT_ZETA = False case.

Can you try outputting the same result you have above (for point 4., with Gamma and MFP), except changing FRACT_FLOAT_ERR just in that specific line of code where Gamma and MFP are generated. Say trying 1e-10 and 1e-4 (default value is 1e-7). I'm wondering if the ionisation critera just happens to be right on the numerical accuracy threshold, and thus the ring structure might be caused by this criteria being met for different radii.

I think the ring structure has to be sourced directly by the ionisation criteria. It might also be worth plotting the ionisation critera as a function of smoothing radius. I would plot both f_coll * ION_EFF_FACTOR and (1. - xHII_from_xrays)*(1.0+rec) as a function of smoothing radii and also a binary map of when f_coll * ION_EFF_FACTOR > (1. - xHII_from_xrays)*(1.0+rec) for each smoothing radii.

I think you are making great progress in isolating this issue.

[qyx268]

It did appear in USE_MASS_DEPENDENT_ZETA = False. In fact, the recent tests were all done with USE_MASS_DEPENDENT_ZETA = False.

The problem seems indeed in N_rec_filtered, I found dNrec_box can be negative (in both the 1Gpc boxes and also the smaller boxes where rings dont appear). In case where rings appear, it has one cell with dNrec_box=106870880000 and xH=1.

I'm gonna go get lunch, a quick thought is that this might be caused by not initializing previous_ionize_box->dNrec_box when USE_MINI_HALOS=False... However, presumably, this is not the reason here since for the first redshift (box->mean_f_coll * ION_EFF_FACTOR + box->mean_f_coll_MINI * ION_EFF_FACTOR_MINI< global_params.HII_ROUND_ERR) is likely true, and then in the next snapshot, previous_ionize_box->dNrec_box will be initialized to 0 in the python part... But let's see...

I was also thinking about FRACT_FLOAT_ERR. I will test it after finding out why dNrec_box can be negative and >1e10...

[BradGreig]

Ah oops, I misread your earlier message about USE_MASS_DEPENDENT_ZETA = False.

Interesting, this single large value cell for dNrec_box is likely the problem. Unfortunately you'll need to track that specific cell all the way through to find out what causes it (if that is even trivial, I don't envy that task!). It is probably occurring when dNrec is calculated (outside of the R loop, near the end of IonisationBox). I wonder if it is stepping out of the interpolation table range for some specific value.

[andreimesinger]

great job isolating this!

On 05.02.2021, at 03:17, Yuxiang Qin notifications@github.com wrote:

It did appear in USE_MASS_DEPENDENT_ZETA = False. In fact, the recent tests were all done with USE_MASS_DEPENDENT_ZETA = False.

The problem seems indeed in N_rec_filtered, I found dNrec_box can be negative (in both the 1Gpc boxes and also the smaller boxes where rings dont appear). In case where rings appear, it has one cell with dNrec_box=106870880000 and xH=1.

I'm gonna go get lunch, a quick thought is that this might be caused by not initializing previous_ionize_box->dNrec_box when USE_MINI_HALOS=False... However, presumably, this is not the reason here since for the first redshift (box->mean_f_coll ION_EFF_FACTOR + box->mean_f_coll_MINI ION_EFF_FACTOR_MINI< global_params.HII_ROUND_ERR) is likely true, and then in the next snapshot, previous_ionize_box->dNrec_box will be initialized to 0 in the python part... But let's see...

I was also thinking about FRACT_FLOAT_ERR. I will test it after finding out why dNrec_box can be negative and >1e10...

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/21cmfast/21cmFAST/issues/194#issuecomment-775604981, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADH2UARLNJMJ767C3XJSROTS6CLNLANCNFSM4UMM4GLA.