BUG: artificacts in vorticity for athena

[chrishavlin]

Plots of vorticity (and presumably any variable that needs a derivative) have artifacts aligned with grid boundaries where the resolution changes:

import yt 
ds = yt.load_sample("AM06")
slc = yt.SlicePlot(ds, 'x', ('gas', 'vorticity_x'), window_size=(4,4))
slc.save('vort_x.png')
slc.annotate_grids()
slc.save('vort_x_with_grids.png')

Also note that you do get a RuntimeWarning:

/home/chavlin/src/yt_general/yt/yt/data_objects/construction_data_containers.py:1531: RuntimeWarning: 
Something went wrong during field computation. This is likely due to missing ghost-zones support in 
class <class 'weakref.ProxyType'>

So presumably this is a known issue, but I didn't see a clear corresponding open issue so I thought it was worth opening a new report to document it. There are, however, a number of open and closed issues that are related to ghost zones and vorticity (in athena and other frontends):

• Handling athena file with ghost zones #2033
• fix: log error instead of crashing in smoothed grid computation if ghost-zones can't be reconstructed #2715
• Issue plotting vorticity from amrvac data #2710
• BUG: [AMRVAC] datafiles with ghost cells are not read correctly #3382
• The AMReX frontend crashes when attempting to read plotfiles with ghost cells #1419

[chrishavlin]

also, feel free to close this as a duplicate if I missed an open issue or if one of the linked ones handles this sufficiently :)

[matthewturk]

That is a super weird error I do not recognize. I believe that what we should look at is whether or not the "smoothed" covering grids are being used.

[chrishavlin]

I believe that what we should look at is whether or not the "smoothed" covering grids are being used.

That would be a yes, this call to retrieve_ghost_zones hardcodes a smoothed=True parameter, which results in using a smoothed covering grid to fill the ghost zones

yt/yt/geometry/grid_geometry_handler.py

Lines 380 to 384 in 9e483ce

	for og in giter:
	if ngz > 0:
	g = og.retrieve_ghost_zones(ngz, [], smoothed=True)
	else:
	g = og

[chrishavlin]

setting it to False to see what happens... it does run, and the magnitude of the artifacts is lower, but they're still there:

[pgrete]

This is equally applicable to Parthenon data.
Is there a default way to prolongate data from coarser grids into ghost zone of finer grids?
Alternatively, one could use one-sided derivatives at level boundaries.

[chrishavlin]

Is there a default way to prolongate data from coarser grids into ghost zone of finer grids?

I thought that's what would happen by setting that smoothed parameter to false up in my previous comment -- when false, it will use a covering grid to fill ghost cells for each grid, which should simply sample the coarser grid(s) containing a finer grid. Relevant code is at

yt/yt/data_objects/index_subobjects/grid_patch.py

Lines 239 to 268 in fb8d4dc

	def retrieve_ghost_zones(self, n_zones, fields, all_levels=False, smoothed=False):
	# We will attempt this by creating a datacube that is exactly bigger
	# than the grid by nZones*dx in each direction
	nl = self.get_global_startindex() - n_zones
	new_left_edge = nl * self.dds + self.ds.domain_left_edge
	# Something different needs to be done for the root grid, though
	level = self.Level
	if all_levels:
	level = self.index.max_level + 1
	kwargs = {
	"dims": self.ActiveDimensions + 2 * n_zones,
	"num_ghost_zones": n_zones,
	"use_pbar": False,
	"fields": fields,
	}
	# This should update the arguments to set the field parameters to be
	# those of this grid.
	field_parameters = {}
	field_parameters.update(self.field_parameters)
	if smoothed:
	cube = self.ds.smoothed_covering_grid(
	level, new_left_edge, field_parameters=field_parameters, **kwargs
	)
	else:
	cube = self.ds.covering_grid(
	level, new_left_edge, field_parameters=field_parameters, **kwargs
	)
	cube._base_grid = self
	return cube

I stepped through it all and didn't notice anything obvious going wrong, but it might be easier to test with a simpler test dataset.