[yt-users] Projection issue

[Elizabeth Tasker]

Hi Britton,

Thanks! That also explains why I previously had everything working with the projected mass (unigrid data).

Elizabeth

On Apr 4, 2013, at 11:32 PM, Britton Smith <brittonsmith at gmail.com> wrote:

> Hi Elizabeth,
> 
> The problem with calculating the mass from a projection of the mass field is that you need to divide each projection pixel by the total path length that cell has been projected.  With AMR, this is trickier because cells with more refined child cells won't be used in the projection, so for each value in proj you don't know the total path.  This does work with unigrid data, so for example, this will give the right total mass:
> (proj['CellMassMsun'] / proj['pdx']).sum() / pf.units['cm'] / 2
> 
> However, the following will work:
> ((2 * pf.units['cm'] * proj['pdx'])**2 * proj['Density']).sum() / 1.989e33
> since the Density projection gives the column density.
> 
> I'm not sure what's going on with the fixed_res_proj.
> 
> Britton
> 
> On Thu, Apr 4, 2013 at 7:58 AM, Elizabeth Tasker <tasker at astro1.sci.hokudai.ac.jp> wrote:
> Hello,
> 
> I have 2 problems, but possibly I only need the solution to one :)
> 
> (1) If I calculate the total mass in a simulation via the projection:
> 
> projection = pf.h.proj(2, 'CellMass')
> sum(projection["CellMass"]/projection["pdx"]/2.0/3.08568025e18)
> 
> (i.e. sum(sum(m_i*dl)/dl) )
> 
> I get:  1.4855032022312827e+43
> 
> However, if I just sum it directly:
> 
> dd = pf.h.all_data()
> sum(dd["CellMass"])
> 
> I get: 6.3279045575257765e+42
> 
> I was expecting those two numbers to be approximately the same. Why are they very different? Have I misunderstood the value of pdx? I think it should be equal to half the integration step.
> 
> (2) The second problem occurred when I tried to fix this by using a fixed resolution projection (in case a variable dl was the issue).
> 
> I tried:
> 
> proj = pf.h.fixed_res_proj(2, 0, [0, 0, 0], [128,128,128])
> 
> (My root grid is 128^3)
> 
> but I get:
> 
> ValueError                                Traceback (most recent call last)
> /home/tasker/yt-new/src/yt-hg/scripts/iyt in <module>()
> ----> 1 proj = pf.h.fixed_res_proj(2, 0, [0, 0, 0], [128,128,128])
> 
> /home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in __init__(self, axis, level, left_edge, dims, fields, pf, **kwargs)
>    2460         self.domain_width = np.rint((self.pf.domain_right_edge -
>    2461                     self.pf.domain_left_edge)/self.dds).astype('int64')
> -> 2462         self._refresh_data()
>    2463 
>    2464     def _get_list_of_grids(self):
> 
> /home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in _refresh_data(self)
>     318         """
>     319         self.clear_data()
> --> 320         self.get_data()
>     321 
>     322     def keys(self):
> 
> /home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in get_data(self, fields)
>    2494         Iterates over the list of fields and generates/reads them all.
>    2495         """
> -> 2496         self._get_list_of_grids()
>    2497         if not self.has_key('pdx'):
>    2498             self._generate_coords()
> 
> /home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in _get_list_of_grids(self)
>    2470         else:
>    2471             grids,ind = self.pf.hierarchy.get_box_grids(
> -> 2472                             self.left_edge, self.right_edge)
>    2473         level_ind = (self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)
>    2474         sort_ind = np.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])
> 
> /home/tasker/yt-new/src/yt-hg/yt/data_objects/object_finding_mixin.pyc in get_box_grids(self, left_edge, right_edge)
>     201         eps = np.finfo(np.float64).eps
>     202         grid_i = np.where((np.all((self.grid_right_edge - left_edge) > eps, axis=1)
> --> 203                          & np.all((right_edge - self.grid_left_edge) > eps, axis=1)) == True)
>     204 
>     205         return self.grids[grid_i], grid_i
> 
> ValueError: operands could not be broadcast together with shapes (3,3) (7100,3) 
> 
> 
> What did I do wrong?
> 
> Thank you!
> 
> Elizabeth
> 
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org
> 
> 
> _______________________________________________
> yt-users mailing list
> yt-users at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20130405/5c0ac72a/attachment.html>