For completeness, looking more at the docs, I was able to get some reasonable numbers using arbitrary_grid.  I had to step through small sections of my box to keep the arrays a reasonable size, and carefully set my region size and cell number so that the arbitrary_grid cells were the same size as the cells in my best refinement level.  I also tried to make sure that my arbitrary_grid cells fell on my highly-refined cells so that I wouldn't end up with zeros, but I am not sure how necessary that was.  Then I just read the value I needed into numpy arrays and brute-forced it. 

Best,
Stephanie
--
Dr. Stephanie Tonnesen
Associate Research Scientist
CCA, Flatiron Institute
New York, NY

stonnes@gmail.com


On Thu, Jun 20, 2019 at 4:40 PM Stephanie Tonnesen <stonnes@gmail.com> wrote:
Hi all,

I am having some trouble figuring out how to calculate the mass flux as a function of z-position in a wind-tunnel simulation.  I am running Enzo with 5 AMR levels.  

If I trying and make a 1d-profile with create_profile, I would have to make my bin sizes as small as the highest level of refinement, which means I have plenty of z values with zeros.  

I tried to make a covering_grid, thinking I could just deal with everything as arrays, but as I expected I run out of memory if I make a covering grid of the whole domain at max resolution.  When I create a YTdisk region, I get an error when trying to make a covering_grid:  
AttributeError: 'YTDisk' object has no attribute 'covering_grid'

I then tried to use the surface and calculate_flux commands:

surf3 = ds.disk([0.5,0.5,0.5],[0., 0., 1.],(30,'kpc'),(100,'kpc'))
surfu = ds.surface(surf3,"z",(54,"kpc"))
fluxu = surfu.calculate_flux("x-velocity","y-velocity","z-velocity","density")

at which point I get this error:
____________________________________________________
Traceback (most recent call last):
  File "yt_surface_flux.py", line 85, in <module>
    fluxu = surfu.calculate_flux("x-velocity","y-velocity","z-velocity","density")
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/construction_data_containers.py", line 1269, in calculate_flux
    field_x, field_y, field_z, fluxing_field)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/construction_data_containers.py", line 1280, in _calculate_flux_in_grid
    vc_data = grid.get_vertex_centered_data(vc_fields)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/grid_patch.py", line 310, in get_vertex_centered_data
    cg = self.retrieve_ghost_zones(1, fields, smoothed=smoothed)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/grid_patch.py", line 272, in retrieve_ghost_zones
    **kwargs)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/construction_data_containers.py", line 927, in __init__
    YTCoveringGrid.__init__(self, *args, **kwargs)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/construction_data_containers.py", line 555, in __init__
    self.get_data(fields)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/construction_data_containers.py", line 642, in get_data
    if len(fill) > 0: self._fill_fields(fill)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/construction_data_containers.py", line 1003, in _fill_fields
    for chunk in ls.data_source.chunks(fields, "io"):
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/data_containers.py", line 1263, in chunks
    self.get_data() # Ensure we have built ourselves
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/data_containers.py", line 1306, in get_data
    self.index._identify_base_chunk(self)
  File "/mnt/home/stonnesen/yt-conda/yt/yt/geometry/grid_geometry_handler.py", line 289, in _identify_base_chunk
    gi = dobj.selector.select_grids(self.grid_left_edge,
  File "/mnt/home/stonnesen/yt-conda/yt/yt/data_objects/data_containers.py", line 1258, in selector
    self._selector = sclass(self)
  File "yt/geometry/selection_routines.pyx", line 860, in yt.geometry.selection_routines.RegionSelector.__init__
RuntimeError: Error: yt attempted to read outside the boundaries of a non-periodic domain along dimension 0.
Region left edge = -0.0078125 code_length, Region right edge = 1.0078125 code_length
Dataset left edge = 0.0 code_length, Dataset right edge = 1.0 code_length

This commonly happens when trying to compute ghost cells up to the domain boundary. Two possible solutions are to load a smaller regi
on that does not border the edge or override the periodicity for this dataset.
___________________________________________

However, my simulated region is 300 kpc across, so I am not at the edge in any axis.  

I suppose I have 2 questions:  can I make a covering_grid of a region of a simulation, or can I have some guidance on making the calculate_flux command work?

Thanks in advance!!
Best,
Stephanie

--
Dr. Stephanie Tonnesen
Associate Research Scientist
CCA, Flatiron Institute
New York, NY

stonnes@gmail.com