[Yt-dev] RMS mass overdensity method

Matthew Turk matthewturk at gmail.com
Sun Mar 1 19:57:45 PST 2009


Hi Stephen,

What you've written looks good.  But, I've been mulling over your
procedure a little bit, and I wonder if we can do a couple things
differently.

For starters, I'd prefer not to have any MPI calls in the object
itself.  This breaks the YT philosophy of transaction-based
parallelism; I believe we can get around this quite easily, however,
so that should not be a problem.  You definitely want to use spheres,
but I think you can achieve better parallelism if you use Derived
Quantities as well.

I don't think we need to partition the hierarchy.  You're right,
domain decomposition does make things go faster if you load the data
all at once.  From what I can tell, you don't do that here.  What I
would propose would be one of two solutions, depending on the
structure of the grids you think you will be dealing with (and the
structure of the grids we can reasonably expect people using this
method to have.)

(point of order: We need to implement a periodic sphere.  This is not
too hard; if you go over any one of the Domains, you just iterate over
that inside each of the cuttings in the AMRSphereBase object.  I can
take care of this if nobody else wants to.)

Possibility 1:

If you expect that the spheres you take will mostly be inside a
*single grid* at a time, then we can decompose spheres across
processors.  We use the ParallelGridIterator to do this, but instead
of iterating over grids, we iterate over spheres.  I'll actually
probably subclass this to make it work with arbitrary objects.  On
each sphere, you would want to call
quantities["TotalQuantity"]("ParticleMassMsun", lazy_reader=False).
Once done, you'd use one of the ParallelAnalysisInterface methods to
calculate the final sigma 8.

Possibility 2:

If you believe that each sphere will intersect a number of grids on
the order of the number of processors, a *much* simpler method can be
written.  You simply iterate over every single sphere you want to
generate, use the TotalQuantity method above, but this time set
lazy_reader=True.  This will take care of decomposition and joining
over processors and broadcast the result to all procs.  Then you can
just divide as appropriate.

To calculate the total mass in the box you definitely want to use the
quantity TotalMass, lazy_reader=True.  Additionally, because these
methods calculate the mass assuming the particle fully occupies a
cell, you should probably avoid using any radius^3 arguments to
calculate the volume, or you will get intersection errors.

I can try to write some of this code up, if it's not clear.  What do you think?

-Matt

> 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
> ________________________________(_)_\(_)_______________
>
> _______________________________________________
> Yt-dev mailing list
> Yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org
>



More information about the yt-dev mailing list