[Yt-dev] RMS mass overdensity method

[Stephen Skory]

Hi all,

I am currently attempting to write a (parallel) RMS mass overdensity method. This is the mass variance over the average density smoothed over a length scale. sigma_8 is this quantity smoothed over 8 Mpc at z=0 in a linearly collapsed universe, for example.

I have a method that works (mostly), shown in the link below. It first finds the total mass in the box in order to calculate the average density field. It then picks a random center and gets the particle mass inside of a sphere at that center of the chosen length scale, normalizing by the average mass in that sphere. It does this many times, and then the RMS is calculated from the samples.

http://paste.enzotools.org/show/60/

Some of you may notice a flaw in the code: if the random center is near to the box edge, part of the sphere will be outside the box, giving too low masses. I hope to address this, see below. If there are other flaws, let a brother know!

Since I wish to use this on the L7 dataset, doing this in parallel is a must. And luckily, this proceedure requires very similar logic to parallel HOP. Spatial decomposition along with periodic oversampling (padding) will work, I think. The padding, which will be the radius of the smoothing spheres, will make sure that any center chosen in the 'real' box will have it's full compliment of particles inside of it, very similar to the parallel HOP concept.

Below is a rough attempt at how I think a parallel implementation can work, but there are some issues I would like some advice from the yt cognesceti:

http://paste.enzotools.org/show/61/

- I calculate the total mass on lines 61-65 by partitioning the box just like parallel HOP does, and there's no over counting. But when I get into the individual threads, on line 18 I'm trying to oversample for each sub-box. I've modified _partition_hierarchy_3d (see below) so that in a special case (like this) we can oversample the entire box in serial mode, although this doesn't work yet. What I'd like for serial mode, for example, are particle lists with oversampled reflected particles such that the range of points is [-padding, 1+padding]. But I'm not sure this is the right way to go, modifying _parititon_hierarchy_3d.

 
http://paste.enzotools.org/show/62/

- in _get_over_density(), I am using the pf.h.sphere() method to cut out the particles inside the smoothing sphere, but I'm not sure I want to do this. I'm aware that as it is it won't work using the _partition_hierarchy_3d method. Can you give me a suggestion of the 'right' way to do this? I'd wonder if investing time in the KDtree scipy stuff Matt emailed about earlier is worth it, but we already have this domain decomposition plus oversampling machinery, and we definitely need the oversampling for this method. And domain decomposition makes data reading so much faster.

Thanks for your comments!


_______________________________________________________
sskory at physics.ucsd.edu           o__  Stephen Skory
http://physics.ucsd.edu/~sskory/ _.>/ _Graduate Student
________________________________(_)_\(_)_______________