search

sbird / fake_spectra

A code for generating fake spectra from a cosmological simulation
MIT License
12 stars 13 forks source link

Error with Voronoi interpolation #20

Closed sbird closed 5 years ago

sbird commented 5 years ago

The code when run on one of the Illustris-TNG simulations gives:

i = 327682, something wrong left = 32759, right = 32768

Which seems to be related to the interpolation in index_table.cpp, within assign_cells.

There is a condition

        if(xp - arr2[2*min_ind+1] < 1.01*reso)
            arr2[2*min_ind+1] = xp;

which looks inverted: we want to increase the right upwards boundary, not decrease it. So I think maybe it should be 

        if(xp - arr2[2*min_ind+1] > 1.01*reso)
            arr2[2*min_ind+1] = xp;
sbird commented 5 years ago

@xiaohanzai what do you think?

sbird commented 5 years ago

@mahdiqezlou Does this work?

xiaohanzai commented 5 years ago

I think my original intention is that the increase in arr2[2 x min_ind+1] should not be larger than 1 reso. Suppose a segment of the sightline is already assigned to cell i. Then when you move rightward, you find the next segment belonging to cell j. Now, if in the next step we find the point xp belonging to cell i again, this shouldn't happen (unless we are close to the end of the sightline and the periodic boundary condition should be taken into account). So the condition if(xp - arr2[2 x min_ind+1] < 1.01 x reso) is to prevent this situation from happening.

There doesn't seem to be enough output to see where the problem comes from. Maybe try to print out min_ind and see how arr2[2 x min_ind] and arr2[2 x min_ind+1] evolve?

I did have problems with this if condition myself, when I ran the code with a 768^3 simulation. The precision of arr2 seems strange, but I'm not sure what's causing this. Changing from float to double didn't help. So I increased 1.01 x reso to 1.1 x reso and it solved my problem. I think this still gives the correct assignment to cells.

qezlou commented 5 years ago

@xiaohanzai Thank you for your response! I am going to run it again with those printings to see what will be the outcome.

sbird commented 5 years ago

I see. So this condition is really a consistency check, an assertion. I can see that this could be a floating point precision issue, but it doesn't make much sense to me either if the change from float to double does not help. Thanks for explaining your purpose here!

qezlou commented 5 years ago

I ran the fake_spectra on Illustris-TNG once more after modifying the 1.01 reso to 1.1 reso as @xiaohanzai mentioned above. The code seems working now. As she said, I assume this modification still gives the correct result.

xiaohanzai commented 5 years ago

@mahdiqezlou That's great to know! Thanks! @sbird Unfortunately I was never able to figure out why the if condition fails. It only became a problem for me when I ran the code with the larger sims. 512^3 and 256^3 both went fine. Maybe it's ok to put 1.1 x reso as a temporary fix for now?

sbird commented 5 years ago

@mahdiqezlou @xiaohanzai I pushed a slightly more comprehensive fix here: https://github.com/sbird/fake_spectra/pull/21 the main change is that if there is floating point roundoff, the condition will not be false when it should be true, instead of true when it should be false.

xiaohanzai commented 5 years ago

@sbird Thanks a lot! I'll test this with my simulations.