[yt-dev] data on different AMR levels

[Matthew Turk]

On Wed, Jan 16, 2013 at 9:52 PM, Elizabeth Tasker
<tasker at astro1.sci.hokudai.ac.jp> wrote:
> Hi Mike and Matt,
>
> Thank you! Both those solutions look great.
>
> I like Mike's a lot for its transparency, but possibly Matt's might end up making slightly neater code.
>
> One question, Matt:
>
>> Okay, sorry for being dense.  But basically what I'm getting is that
>> you want the functionality of covering grids, but you want to be able
>> to restrict the grids from going above some level L, right?
>
> Or below. I just want cells on one level only.

Okay, the grid collection might work, but I still think the covering
grid is the best, with the masking I describe below.

gc = pf.h.grid_collection([0.5, 0.5, 0.5], pf.h.select_grids(L))

gc["Density"]
gc["dx"]
gc["x"]

etc.  You'll have to manually manipulate the dx stuff, cell mass, etc.

>
>>
>> And for the Level 1 case you mention above, do you want contributions
>> from Level 0?  (That's not so hard to do.)
>
> No. So I want an array that contains only cells from level 1, with no contributions from 0 or 2.

Okay, then mask or use the grid collection.  :)

>
>> To do this with a covering grid you could pre-allocate:
>>
>> Li = 0
>> Lo = 3
>>
>> cg = pf.h.covering_grid(Li, [0, 0, 0], pf.domain_dimensions * pf.refine_by**Li)
>> output = np.empty(pf.domain_dimensions * pf.refine_by**Lo)
>> nd = pf.refine_by**(Lo-li)
>> for i in range(nd):
>>    for j in range(nd):
>>        for k in range(nd):
>>            output[i::nd,j::nd,k::nd] = cg[field]
>
>
> Ah ok! So this is the same as Mike's method in that a covering_grid is created at the resolution of the cell size on the specified level, and then the array values are multiplied to give the correct number for a higher resolution? (Which works fine for the basic quantities; density, energy and velocity but won't for cell mass -- which is totally fine).

Cell mass is converted correctly, I believe.  This is accomplished
through a validator, NeedsGridType, which is specified if a field
*needs* to be calculated in the original grid (i.e., Level) but
otherwise fields are generated in the CG.  I've verified that
specifically for CellMassMsun this is correctly done:

In [12]: ms = []

In [13]: for L in range(4):
    cg = pf.h.covering_grid(L, [0,0,0], pf.domain_dimensions * pf.refine_by**L)
    ms.append(cg["CellMassMsun"].sum())
   ....:
yt : [DEBUG    ] 2013-01-16 22:00:26,392 Appending object to
galaxy0030 (type: <class 'yt.data_objects.hierarchy.AMRCoveringGrid'>)
yt : [DEBUG    ] 2013-01-16 22:00:26,436 Getting fields ['Density']
from 1 possible grids
yt : [DEBUG    ] 2013-01-16 22:00:26,439 Appending object to
galaxy0030 (type: <class 'yt.data_objects.hierarchy.AMRCoveringGrid'>)
yt : [DEBUG    ] 2013-01-16 22:00:26,442 Getting fields ['Density']
from 9 possible grids
yt : [DEBUG    ] 2013-01-16 22:00:26,465 Appending object to
galaxy0030 (type: <class 'yt.data_objects.hierarchy.AMRCoveringGrid'>)
yt : [DEBUG    ] 2013-01-16 22:00:26,518 Getting fields ['Density']
from 17 possible grids
yt : [DEBUG    ] 2013-01-16 22:00:26,703 Appending object to
galaxy0030 (type: <class 'yt.data_objects.hierarchy.AMRCoveringGrid'>)
yt : [DEBUG    ] 2013-01-16 22:00:27,320 Getting fields ['Density']
from 25 possible grids

In [14]: ms
Out[14]:
[18524964522.745506,
 18524964453.698158,
 18524964476.255985,
 18524964484.909931]

>
> My only question is:
>
> If I do this for level (e.g.) 2 and level 2 doesn't cover the whole domain -- what does the covering grid do when there are no level 2 cells over a region of space?

Fills 'em with level 1 or level 0, depending.  You can just mask these
out, though, which might be the best solution if the covering grid
(3D, handles the mass issue, etc) is the better solution:

bag_ind = cg["GridLevels"] != 2
cg[field][bad_ind] = -1.0 # or something

That should do it, and that would go somewhere inside the loop above.

-Matt

>
> Thanks again!
>
> Elizabeth
> _______________________________________________
> yt-dev mailing list
> yt-dev at lists.spacepope.org
> http://lists.spacepope.org/listinfo.cgi/yt-dev-spacepope.org