<div dir="ltr"><div><div><div><div>Hi Elizabeth,<br><br></div>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:<br>
(proj['CellMassMsun'] / proj['pdx']).sum() / pf.units['cm'] / 2<br></div></div><div><br></div>However, the following will work:<br>((2 * pf.units['cm'] * proj['pdx'])**2 * proj['Density']).sum() / 1.989e33<br>
</div>since the Density projection gives the column density.<br><div class="gmail_extra"><br></div><div class="gmail_extra">I'm not sure what's going on with the fixed_res_proj.<br><br></div><div class="gmail_extra">
Britton<br></div><div class="gmail_extra"><br><div class="gmail_quote">On Thu, Apr 4, 2013 at 7:58 AM, Elizabeth Tasker <span dir="ltr"><<a href="mailto:tasker@astro1.sci.hokudai.ac.jp" target="_blank">tasker@astro1.sci.hokudai.ac.jp</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hello,<div><br></div><div>I have 2 problems, but possibly I only need the solution to one :)</div><div><br>
</div><div>(1) If I calculate the total mass in a simulation via the projection:</div>
<div><br></div><div>projection = pf.h.proj(2, 'CellMass')<br></div><div>sum(projection["CellMass"]/projection["pdx"]/2.0/3.08568025e18)<br></div><div><br></div><div>
(i.e. sum(sum(m_i*dl)/dl) )</div><div><br></div><div>I get: 1.4855032022312827e+43</div><div><br></div><div>However, if I just sum it directly:</div><div><br></div><div>dd = pf.h.all_data()</div>
<div>sum(dd["CellMass"])<br></div><div><br></div><div>I get: 6.3279045575257765e+42</div><div><br></div><div>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.</div>
<div><br></div><div>(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).</div><div><br></div><div>I tried:</div><div>
<br></div><div>proj = pf.h.fixed_res_proj(2, 0, [0, 0, 0], [128,128,128])<br></div><div><br></div><div>(My root grid is 128^3)</div><div><br></div><div>but I get:</div><div><br></div><div>
<div>ValueError Traceback (most recent call last)</div><div>/home/tasker/yt-new/src/yt-hg/scripts/iyt in <module>()</div><div>----> 1 proj = pf.h.fixed_res_proj(2, 0, [0, 0, 0], [128,128,128])</div>
<div><br></div><div>/home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in __init__(self, axis, level, left_edge, dims, fields, pf, **kwargs)</div><div> 2460 self.domain_width = np.rint((self.pf.domain_right_edge -</div>
<div> 2461 self.pf.domain_left_edge)/self.dds).astype('int64')</div><div>-> 2462 self._refresh_data()</div><div> 2463 </div><div> 2464 def _get_list_of_grids(self):</div><div>
<br></div><div>/home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in _refresh_data(self)</div><div> 318 """</div><div> 319 self.clear_data()</div><div>--> 320 self.get_data()</div>
<div> 321 </div><div> 322 def keys(self):</div><div><br></div><div>/home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in get_data(self, fields)</div><div> 2494 Iterates over the list of fields and generates/reads them all.</div>
<div> 2495 """</div><div>-> 2496 self._get_list_of_grids()</div><div> 2497 if not self.has_key('pdx'):</div><div> 2498 self._generate_coords()</div><div>
<br></div><div>/home/tasker/yt-new/src/yt-hg/yt/data_objects/data_containers.pyc in _get_list_of_grids(self)</div><div> 2470 else:</div><div> 2471 grids,ind = self.pf.hierarchy.get_box_grids(</div>
<div>-> 2472 self.left_edge, self.right_edge)</div><div> 2473 level_ind = (self.pf.hierarchy.grid_levels.ravel()[ind] <= self.level)</div><div> 2474 sort_ind = np.argsort(self.pf.h.grid_levels.ravel()[ind][level_ind])</div>
<div><br></div><div>/home/tasker/yt-new/src/yt-hg/yt/data_objects/object_finding_mixin.pyc in get_box_grids(self, left_edge, right_edge)</div><div> 201 eps = np.finfo(np.float64).eps</div><div> 202 grid_i = np.where((np.all((self.grid_right_edge - left_edge) > eps, axis=1)</div>
<div>--> 203 & np.all((right_edge - self.grid_left_edge) > eps, axis=1)) == True)</div><div> 204 </div><div> 205 return self.grids[grid_i], grid_i</div><div><br></div><div>
ValueError: operands could not be broadcast together with shapes (3,3) (7100,3) </div><div><br></div><div><br></div><div>What did I do wrong?</div><div><br></div><div>Thank you!</div><span class="HOEnZb"><font color="#888888"><div>
<br></div>
<div>Elizabeth</div></font></span></div></div>
<br>_______________________________________________<br>
yt-users mailing list<br>
<a href="mailto:yt-users@lists.spacepope.org">yt-users@lists.spacepope.org</a><br>
<a href="http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org" target="_blank">http://lists.spacepope.org/listinfo.cgi/yt-users-spacepope.org</a><br>
<br></blockquote></div><br></div></div>