[yt-svn] commit/yt: 20 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon May 1 12:38:23 PDT 2017
20 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/a5764a02ce44/
Changeset: a5764a02ce44
Branch: yt
User: atmyers
Date: 2017-04-03 20:46:38+00:00
Summary: initial work on nodal field handling.
Affected #: 4 files
diff -r 3fee41c2ea35fbe17c97ba1dc286d6e71f088716 -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -90,7 +90,9 @@
return tr
finfo = self.ds._get_field_info(*fields[0])
if not finfo.particle_type:
- return tr.reshape(self.ActiveDimensions)
+ num_nodes = 2**sum(finfo.nodal_flag)
+ new_shape = list(self.ActiveDimensions) + [num_nodes]
+ return tr.reshape(new_shape)
return tr
def convert(self, datatype):
@@ -374,11 +376,35 @@
self._last_count = mask.sum()
return mask
+ def _get_nodal_slices(self, shape, nodal_flag):
+ slices = []
+ dir_slices = [[] for _ in range(3)]
+
+ for i in range(3):
+ if nodal_flag[i]:
+ dir_slices[i] = [slice(0, shape[i]-1), slice(1, shape[i])]
+ else:
+ dir_slices[i] = [slice(0, shape[i])]
+
+ for i, sl_i in enumerate(dir_slices[0]):
+ for j, sl_j in enumerate(dir_slices[1]):
+ for k, sl_k in enumerate(dir_slices[2]):
+ slices.append([sl_i, sl_j, sl_k])
+
+ return slices
+
def select(self, selector, source, dest, offset):
+ nodal_flag = source.shape - self.ActiveDimensions
mask = self._get_selector_mask(selector)
count = self.count(selector)
if count == 0: return 0
- dest[offset:offset+count] = source[mask]
+ nodal_flag = source.shape - self.ActiveDimensions
+ if sum(nodal_flag) == 0:
+ dest[offset:offset+count] = source[mask]
+ else:
+ slices = self._get_nodal_slices(source.shape, nodal_flag)
+ for i , sl in enumerate(slices):
+ dest[offset:offset+count, i] = source[sl][mask]
return count
def count(self, selector):
@@ -386,6 +412,15 @@
if mask is None: return 0
return self._last_count
+ def count_nodal(self, selector, nodal_flag):
+ # We don't cache the selector results
+ mask = self._get_selector_mask(selector)
+ if mask is None: return 0
+ pad_width = [(0, nodal) for nodal in nodal_flag]
+ mask = np.pad(mask, pad_width, 'edge')
+ count = mask.sum()
+ return count
+
def count_particles(self, selector, x, y, z):
# We don't cache the selector results
count = selector.count_points(x,y,z, 0.0)
diff -r 3fee41c2ea35fbe17c97ba1dc286d6e71f088716 -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf yt/fields/derived_field.py
--- a/yt/fields/derived_field.py
+++ b/yt/fields/derived_field.py
@@ -110,6 +110,8 @@
self.vector_field = vector_field
self.ds = ds
+ self.nodal_flag = [0, 0, 0]
+
self._function = function
if validators:
diff -r 3fee41c2ea35fbe17c97ba1dc286d6e71f088716 -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf yt/frontends/stream/io.py
--- a/yt/frontends/stream/io.py
+++ b/yt/frontends/stream/io.py
@@ -48,7 +48,14 @@
raise NotImplementedError
rv = {}
for field in fields:
- rv[field] = self.ds.arr(np.empty(size, dtype="float64"))
+ nodal_flag = self.ds.nodal_flags[field[1]]
+ num_nodes = 2**sum(nodal_flag)
+ if num_nodes > 0:
+ rv[field] = self.ds.arr(np.empty((size, num_nodes),
+ dtype="float64"))
+ else:
+ rv[field] = self.ds.arr(np.empty(size, dtype="float64"))
+
ng = sum(len(c.objs) for c in chunks)
mylog.debug("Reading %s cells of %s fields in %s blocks",
size, [f2 for f1, f2 in fields], ng)
diff -r 3fee41c2ea35fbe17c97ba1dc286d6e71f088716 -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -128,10 +128,10 @@
period = period.in_units("code_length").d
buff = np.zeros((size[1], size[0]), dtype="f8")
pixelize_cartesian(buff, data_source['px'], data_source['py'],
- data_source['pdx'], data_source['pdy'],
- data_source[field],
- bounds, int(antialias),
- period, int(periodic))
+ data_source['pdx'], data_source['pdy'],
+ data_source[field],
+ bounds, int(antialias),
+ period, int(periodic))
return buff
def _oblique_pixelize(self, data_source, field, bounds, size, antialias):
https://bitbucket.org/yt_analysis/yt/commits/f8dc3569308d/
Changeset: f8dc3569308d
Branch: yt
User: atmyers
Date: 2017-04-04 23:24:13+00:00
Summary: a pixelizer for nodal data.
Affected #: 3 files
diff -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf -r f8dc3569308d3406468f7813b9c6634388db2161 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -268,7 +268,7 @@
_top_node = "/Slices"
_type_name = "slice"
_con_args = ('axis', 'coord')
- _container_fields = ("px", "py", "pdx", "pdy")
+ _container_fields = ("px", "py", "pz", "pdx", "pdy", "pdz")
def __init__(self, axis, coord, center=None, ds=None,
field_parameters=None, data_source=None):
YTSelectionContainer2D.__init__(self, axis, ds,
@@ -285,13 +285,59 @@
return self._current_chunk.fcoords[:,xax]
elif field == "py":
return self._current_chunk.fcoords[:,yax]
+ elif field == "pz":
+ return self._current_chunk.fcoords[:,self.axis]
elif field == "pdx":
return self._current_chunk.fwidth[:,xax] * 0.5
elif field == "pdy":
return self._current_chunk.fwidth[:,yax] * 0.5
+ elif field == "pdz":
+ return self._current_chunk.fwidth[:,self.axis] * 0.5
else:
raise KeyError(field)
+ _index_map = np.array([[0, 0, 0, 0, 0, 0, 0, 0],
+ [0, 1, 0, 1, 0, 1, 0, 1],
+ [0, 0, 1, 1, 0, 0, 1, 1],
+ [0, 1, 2, 3, 0, 1, 2, 3],
+ [0, 0, 0, 0, 1, 1, 1, 1],
+ [0, 1, 0, 1, 2, 3, 2, 3],
+ [0, 0, 1, 1, 2, 2, 3, 3],
+ [0, 1, 2, 3, 4, 5, 6, 7]])
+
+ def _get_linear_index(self, nodal_flag):
+ return 1*nodal_flag[2] + 2*nodal_flag[1] + 4*nodal_flag[0]
+
+ def _get_indices(self, nodal_flag):
+ li = self._get_linear_index(nodal_flag)
+ return self._index_map[li]
+
+ def _get_nodal_data(self, field):
+ finfo = self.ds._get_field_info(field)
+ nodal_flag = finfo.nodal_flag
+ field_data = self.__getitem__(field)
+ inds = self._get_indices(nodal_flag)
+ return field_data[:, inds]
+
+ def _get_nodal_px_py(self, nodal_flag):
+ xax = self.ds.coordinates.x_axis[self.axis]
+ yax = self.ds.coordinates.y_axis[self.axis]
+ axes = [xax, yax]
+
+ ps = []
+ for ax in axes:
+ coords = self.fcoords[:, ax]
+ if nodal_flag[ax]:
+ fwidth = self.fwidth[:, ax]
+ p = np.empty((coords.shape[0], 2), dtype='float64')
+ p[:, 0] = coords - 0.5*fwidth
+ p[:, 1] = coords + 0.5*fwidth
+ else:
+ p = np.atleast_2d(coords).transpose()
+ ps.append(p)
+
+ return ps
+
@property
def _mrep(self):
return MinimalSliceData(self)
diff -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf -r f8dc3569308d3406468f7813b9c6634388db2161 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -24,7 +24,7 @@
from yt.funcs import mylog
from yt.utilities.lib.pixelization_routines import \
pixelize_element_mesh, pixelize_off_axis_cartesian, \
- pixelize_cartesian
+ pixelize_cartesian, pixelize_cartesian_nodal
from yt.data_objects.unstructured_mesh import SemiStructuredMesh
@@ -120,6 +120,36 @@
def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
dim, periodic):
+ finfo = self.ds._get_field_info(field)
+ nodal_flag = finfo.nodal_flag
+ if np.any(nodal_flag):
+ return self._ortho_pixelize_nodal(data_source, field, bounds, size,
+ antialias, dim, periodic)
+ else:
+ return self._ortho_pixelize_cell(data_source, field, bounds, size,
+ antialias, dim, periodic)
+
+ def _ortho_pixelize_nodal(self, data_source, field, bounds, size, antialias,
+ dim, periodic):
+ # We should be using fcoords
+ period = self.period[:2].copy() # dummy here
+ period[0] = self.period[self.x_axis[dim]]
+ period[1] = self.period[self.y_axis[dim]]
+ if hasattr(period, 'in_units'):
+ period = period.in_units("code_length").d
+ buff = np.zeros((size[1], size[0]), dtype="f8")
+
+ nodal_data = data_source._get_nodal_data(field)
+ coord = data_source.coord.d
+
+ pixelize_cartesian_nodal(buff, data_source['px'], data_source['py'], data_source['pz'],
+ data_source['pdx'], data_source['pdy'], data_source['pdz'],
+ nodal_data, coord, bounds, int(antialias),
+ period, int(periodic))
+ return buff
+
+ def _ortho_pixelize_cell(self, data_source, field, bounds, size, antialias,
+ dim, periodic):
# We should be using fcoords
period = self.period[:2].copy() # dummy here
period[0] = self.period[self.x_axis[dim]]
diff -r a5764a02ce4480fa0bea0c0e3c2e7f06812958bf -r f8dc3569308d3406468f7813b9c6634388db2161 yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -233,6 +233,171 @@
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
+def pixelize_cartesian_nodal(np.float64_t[:,:] buff,
+ np.float64_t[:] px,
+ np.float64_t[:] py,
+ np.float64_t[:] pz,
+ np.float64_t[:] pdx,
+ np.float64_t[:] pdy,
+ np.float64_t[:] pdz,
+ np.float64_t[:, :] data,
+ np.float64_t coord,
+ bounds,
+ int antialias = 1,
+ period = None,
+ int check_period = 1):
+ cdef np.float64_t x_min, x_max, y_min, y_max
+ cdef np.float64_t period_x = 0.0, period_y = 0.0
+ cdef np.float64_t width, height, px_dx, px_dy, ipx_dx, ipx_dy
+ cdef np.float64_t ld_x, ld_y, cx, cy, cz
+ cdef int i, j, p, xi, yi
+ cdef int lc, lr, rc, rr
+ cdef np.float64_t lypx, rypx, lxpx, rxpx, overlap1, overlap2
+ # These are the temp vars we get from the arrays
+ cdef np.float64_t oxsp, oysp, ozsp
+ cdef np.float64_t xsp, ysp, zsp
+ cdef np.float64_t dxsp, dysp, dzsp
+ # Some periodicity helpers
+ cdef int xiter[2]
+ cdef int yiter[2]
+ cdef int ii, jj, kk, ind
+ cdef np.float64_t xiterv[2]
+ cdef np.float64_t yiterv[2]
+ if period is not None:
+ period_x = period[0]
+ period_y = period[1]
+ x_min = bounds[0]
+ x_max = bounds[1]
+ y_min = bounds[2]
+ y_max = bounds[3]
+ width = x_max - x_min
+ height = y_max - y_min
+ px_dx = width / (<np.float64_t> buff.shape[1])
+ px_dy = height / (<np.float64_t> buff.shape[0])
+ ipx_dx = 1.0 / px_dx
+ ipx_dy = 1.0 / px_dy
+ if px.shape[0] != py.shape[0] or \
+ px.shape[0] != pz.shape[0] or \
+ px.shape[0] != pdx.shape[0] or \
+ px.shape[0] != pdy.shape[0] or \
+ px.shape[0] != pdz.shape[0] or \
+ px.shape[0] != data.shape[0]:
+ raise YTPixelizeError("Arrays are not of correct shape.")
+ xiter[0] = yiter[0] = 0
+ xiterv[0] = yiterv[0] = 0.0
+ # Here's a basic outline of what we're going to do here. The xiter and
+ # yiter variables govern whether or not we should check periodicity -- are
+ # we both close enough to the edge that it would be important *and* are we
+ # periodic?
+ #
+ # The other variables are all either pixel positions or data positions.
+ # Pixel positions will vary regularly from the left edge of the window to
+ # the right edge of the window; px_dx and px_dy are the dx (cell width, not
+ # half-width). ipx_dx and ipx_dy are the inverse, for quick math.
+ #
+ # The values in xsp, dxsp, x_min and their y counterparts, are the
+ # data-space coordinates, and are related to the data fed in. We make some
+ # modifications for periodicity.
+ #
+ # Inside the finest loop, we compute the "left column" (lc) and "lower row"
+ # (lr) and then iterate up to "right column" (rc) and "uppeR row" (rr),
+ # depositing into them the data value. Overlap computes the relative
+ # overlap of a data value with a pixel.
+ #
+ # NOTE ON ROWS AND COLUMNS:
+ #
+ # The way that images are plotting in matplotlib is somewhat different
+ # from what most might expect. The first axis of the array plotted is
+ # what varies along the x axis. So for instance, if you supply
+ # origin='lower' and plot the results of an mgrid operation, at a fixed
+ # 'y' value you will see the results of that array held constant in the
+ # first dimension. Here is some example code:
+ #
+ # import matplotlib.pyplot as plt
+ # import numpy as np
+ # x, y = np.mgrid[0:1:100j,0:1:100j]
+ # plt.imshow(x, interpolation='nearest', origin='lower')
+ # plt.imshow(y, interpolation='nearest', origin='lower')
+ #
+ # The values in the image:
+ # lower left: arr[0,0]
+ # lower right: arr[0,-1]
+ # upper left: arr[-1,0]
+ # upper right: arr[-1,-1]
+ #
+ # So what we want here is to fill an array such that we fill:
+ # first axis : y_min .. y_max
+ # second axis: x_min .. x_max
+ with nogil:
+ for p in range(px.shape[0]):
+ xiter[1] = yiter[1] = 999
+ oxsp = px[p]
+ oysp = py[p]
+ ozsp = pz[p]
+ dxsp = pdx[p]
+ dysp = pdy[p]
+ dzsp = pdz[p]
+ if check_period == 1:
+ if (oxsp - dxsp < x_min):
+ xiter[1] = +1
+ xiterv[1] = period_x
+ elif (oxsp + dxsp > x_max):
+ xiter[1] = -1
+ xiterv[1] = -period_x
+ if (oysp - dysp < y_min):
+ yiter[1] = +1
+ yiterv[1] = period_y
+ elif (oysp + dysp > y_max):
+ yiter[1] = -1
+ yiterv[1] = -period_y
+ overlap1 = overlap2 = 1.0
+ zsp = ozsp
+ for xi in range(2):
+ if xiter[xi] == 999: continue
+ xsp = oxsp + xiterv[xi]
+ if (xsp + dxsp < x_min) or (xsp - dxsp > x_max): continue
+ for yi in range(2):
+ if yiter[yi] == 999: continue
+ ysp = oysp + yiterv[yi]
+ if (ysp + dysp < y_min) or (ysp - dysp > y_max): continue
+ lc = <int> fmax(((xsp-dxsp-x_min)*ipx_dx),0)
+ lr = <int> fmax(((ysp-dysp-y_min)*ipx_dy),0)
+ # NOTE: This is a different way of doing it than in the C
+ # routines. In C, we were implicitly casting the
+ # initialization to int, but *not* the conditional, which
+ # was allowed an extra value:
+ # for(j=lc;j<rc;j++)
+ # here, when assigning lc (double) to j (int) it got
+ # truncated, but no similar truncation was done in the
+ # comparison of j to rc (double). So give ourselves a
+ # bonus row and bonus column here.
+ rc = <int> fmin(((xsp+dxsp-x_min)*ipx_dx + 1), buff.shape[1])
+ rr = <int> fmin(((ysp+dysp-y_min)*ipx_dy + 1), buff.shape[0])
+ # Note that we're iterating here over *y* in the i
+ # direction. See the note above about this.
+ for i in range(lr, rr):
+ lypx = px_dy * i + y_min
+ rypx = px_dy * (i+1) + y_min
+ for j in range(lc, rc):
+ lxpx = px_dx * j + x_min
+ rxpx = px_dx * (j+1) + x_min
+
+ cx = (rxpx+lxpx)*0.5
+ cy = (rypx+lypx)*0.5
+ cz = coord
+
+ ii = <int> (cx - xsp + dxsp)
+ jj = <int> (cy - ysp + dysp)
+ kk = <int> (cz - zsp + dzsp)
+
+ ind = 4*ii + 2*jj + kk
+
+ buff[i,j] = data[p, ind]
+
+
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
def pixelize_off_axis_cartesian(
np.float64_t[:,:] buff,
np.float64_t[:] x,
https://bitbucket.org/yt_analysis/yt/commits/bd65a61be0bd/
Changeset: bd65a61be0bd
Branch: yt
User: atmyers
Date: 2017-04-05 00:04:01+00:00
Summary: Move all the nodal utilities into a file.
Affected #: 4 files
diff -r f8dc3569308d3406468f7813b9c6634388db2161 -r bd65a61be0bd73785de1533681d3d35de3141d27 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -29,6 +29,8 @@
YTParticleDepositionNotImplemented
from yt.utilities.lib.interpolators import \
ghost_zone_interpolate
+from yt.utilities.nodal_data_utils import \
+ get_nodal_slices
class AMRGridPatch(YTSelectionContainer):
_spatial = True
@@ -376,23 +378,6 @@
self._last_count = mask.sum()
return mask
- def _get_nodal_slices(self, shape, nodal_flag):
- slices = []
- dir_slices = [[] for _ in range(3)]
-
- for i in range(3):
- if nodal_flag[i]:
- dir_slices[i] = [slice(0, shape[i]-1), slice(1, shape[i])]
- else:
- dir_slices[i] = [slice(0, shape[i])]
-
- for i, sl_i in enumerate(dir_slices[0]):
- for j, sl_j in enumerate(dir_slices[1]):
- for k, sl_k in enumerate(dir_slices[2]):
- slices.append([sl_i, sl_j, sl_k])
-
- return slices
-
def select(self, selector, source, dest, offset):
nodal_flag = source.shape - self.ActiveDimensions
mask = self._get_selector_mask(selector)
@@ -402,7 +387,7 @@
if sum(nodal_flag) == 0:
dest[offset:offset+count] = source[mask]
else:
- slices = self._get_nodal_slices(source.shape, nodal_flag)
+ slices = get_nodal_slices(source.shape, nodal_flag)
for i , sl in enumerate(slices):
dest[offset:offset+count, i] = source[sl][mask]
return count
diff -r f8dc3569308d3406468f7813b9c6634388db2161 -r bd65a61be0bd73785de1533681d3d35de3141d27 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -296,48 +296,6 @@
else:
raise KeyError(field)
- _index_map = np.array([[0, 0, 0, 0, 0, 0, 0, 0],
- [0, 1, 0, 1, 0, 1, 0, 1],
- [0, 0, 1, 1, 0, 0, 1, 1],
- [0, 1, 2, 3, 0, 1, 2, 3],
- [0, 0, 0, 0, 1, 1, 1, 1],
- [0, 1, 0, 1, 2, 3, 2, 3],
- [0, 0, 1, 1, 2, 2, 3, 3],
- [0, 1, 2, 3, 4, 5, 6, 7]])
-
- def _get_linear_index(self, nodal_flag):
- return 1*nodal_flag[2] + 2*nodal_flag[1] + 4*nodal_flag[0]
-
- def _get_indices(self, nodal_flag):
- li = self._get_linear_index(nodal_flag)
- return self._index_map[li]
-
- def _get_nodal_data(self, field):
- finfo = self.ds._get_field_info(field)
- nodal_flag = finfo.nodal_flag
- field_data = self.__getitem__(field)
- inds = self._get_indices(nodal_flag)
- return field_data[:, inds]
-
- def _get_nodal_px_py(self, nodal_flag):
- xax = self.ds.coordinates.x_axis[self.axis]
- yax = self.ds.coordinates.y_axis[self.axis]
- axes = [xax, yax]
-
- ps = []
- for ax in axes:
- coords = self.fcoords[:, ax]
- if nodal_flag[ax]:
- fwidth = self.fwidth[:, ax]
- p = np.empty((coords.shape[0], 2), dtype='float64')
- p[:, 0] = coords - 0.5*fwidth
- p[:, 1] = coords + 0.5*fwidth
- else:
- p = np.atleast_2d(coords).transpose()
- ps.append(p)
-
- return ps
-
@property
def _mrep(self):
return MinimalSliceData(self)
diff -r f8dc3569308d3406468f7813b9c6634388db2161 -r bd65a61be0bd73785de1533681d3d35de3141d27 yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -26,6 +26,7 @@
pixelize_element_mesh, pixelize_off_axis_cartesian, \
pixelize_cartesian, pixelize_cartesian_nodal
from yt.data_objects.unstructured_mesh import SemiStructuredMesh
+from yt.utilities.nodal_data_utils import get_nodal_data
class CartesianCoordinateHandler(CoordinateHandler):
@@ -120,50 +121,33 @@
def _ortho_pixelize(self, data_source, field, bounds, size, antialias,
dim, periodic):
- finfo = self.ds._get_field_info(field)
- nodal_flag = finfo.nodal_flag
- if np.any(nodal_flag):
- return self._ortho_pixelize_nodal(data_source, field, bounds, size,
- antialias, dim, periodic)
- else:
- return self._ortho_pixelize_cell(data_source, field, bounds, size,
- antialias, dim, periodic)
-
- def _ortho_pixelize_nodal(self, data_source, field, bounds, size, antialias,
- dim, periodic):
# We should be using fcoords
period = self.period[:2].copy() # dummy here
period[0] = self.period[self.x_axis[dim]]
period[1] = self.period[self.y_axis[dim]]
if hasattr(period, 'in_units'):
period = period.in_units("code_length").d
+
buff = np.zeros((size[1], size[0]), dtype="f8")
-
- nodal_data = data_source._get_nodal_data(field)
- coord = data_source.coord.d
-
- pixelize_cartesian_nodal(buff, data_source['px'], data_source['py'], data_source['pz'],
- data_source['pdx'], data_source['pdy'], data_source['pdz'],
- nodal_data, coord, bounds, int(antialias),
- period, int(periodic))
+
+ finfo = self.ds._get_field_info(field)
+ nodal_flag = finfo.nodal_flag
+ if np.any(nodal_flag):
+ nodal_data = get_nodal_data(data_source, field)
+ coord = data_source.coord.d
+ pixelize_cartesian_nodal(buff,
+ data_source['px'], data_source['py'], data_source['pz'],
+ data_source['pdx'], data_source['pdy'], data_source['pdz'],
+ nodal_data, coord, bounds, int(antialias),
+ period, int(periodic))
+ else:
+ pixelize_cartesian(buff, data_source['px'], data_source['py'],
+ data_source['pdx'], data_source['pdy'],
+ data_source[field],
+ bounds, int(antialias),
+ period, int(periodic))
return buff
-
- def _ortho_pixelize_cell(self, data_source, field, bounds, size, antialias,
- dim, periodic):
- # We should be using fcoords
- period = self.period[:2].copy() # dummy here
- period[0] = self.period[self.x_axis[dim]]
- period[1] = self.period[self.y_axis[dim]]
- if hasattr(period, 'in_units'):
- period = period.in_units("code_length").d
- buff = np.zeros((size[1], size[0]), dtype="f8")
- pixelize_cartesian(buff, data_source['px'], data_source['py'],
- data_source['pdx'], data_source['pdy'],
- data_source[field],
- bounds, int(antialias),
- period, int(periodic))
- return buff
-
+
def _oblique_pixelize(self, data_source, field, bounds, size, antialias):
indices = np.argsort(data_source['pdx'])[::-1]
buff = np.zeros((size[1], size[0]), dtype="f8")
diff -r f8dc3569308d3406468f7813b9c6634388db2161 -r bd65a61be0bd73785de1533681d3d35de3141d27 yt/utilities/nodal_data_utils.py
--- /dev/null
+++ b/yt/utilities/nodal_data_utils.py
@@ -0,0 +1,41 @@
+import numpy as np
+
+_index_map = np.array([[0, 0, 0, 0, 0, 0, 0, 0],
+ [0, 1, 0, 1, 0, 1, 0, 1],
+ [0, 0, 1, 1, 0, 0, 1, 1],
+ [0, 1, 2, 3, 0, 1, 2, 3],
+ [0, 0, 0, 0, 1, 1, 1, 1],
+ [0, 1, 0, 1, 2, 3, 2, 3],
+ [0, 0, 1, 1, 2, 2, 3, 3],
+ [0, 1, 2, 3, 4, 5, 6, 7]])
+
+def _get_linear_index(nodal_flag):
+ return 1*nodal_flag[2] + 2*nodal_flag[1] + 4*nodal_flag[0]
+
+def _get_indices(nodal_flag):
+ li = _get_linear_index(nodal_flag)
+ return _index_map[li]
+
+def get_nodal_data(data_source, field):
+ finfo = data_source.ds._get_field_info(field)
+ nodal_flag = finfo.nodal_flag
+ field_data = data_source[field]
+ inds = _get_indices(nodal_flag)
+ return field_data[:, inds]
+
+def get_nodal_slices(shape, nodal_flag):
+ slices = []
+ dir_slices = [[] for _ in range(3)]
+
+ for i in range(3):
+ if nodal_flag[i]:
+ dir_slices[i] = [slice(0, shape[i]-1), slice(1, shape[i])]
+ else:
+ dir_slices[i] = [slice(0, shape[i])]
+
+ for i, sl_i in enumerate(dir_slices[0]):
+ for j, sl_j in enumerate(dir_slices[1]):
+ for k, sl_k in enumerate(dir_slices[2]):
+ slices.append([sl_i, sl_j, sl_k])
+
+ return slices
https://bitbucket.org/yt_analysis/yt/commits/3c2c470c05e9/
Changeset: 3c2c470c05e9
Branch: yt
User: atmyers
Date: 2017-04-05 19:25:34+00:00
Summary: teach the boxlib frontend how to read the non-cell-averaged fields.
Affected #: 3 files
diff -r bd65a61be0bd73785de1533681d3d35de3141d27 -r 3c2c470c05e9724cf5276c25c573ea4fe9840688 yt/frontends/boxlib/data_structures.py
--- a/yt/frontends/boxlib/data_structures.py
+++ b/yt/frontends/boxlib/data_structures.py
@@ -1332,6 +1332,55 @@
return vals
+def _read_raw_field_names(raw_file):
+ header_files = glob.glob(raw_file + "*_H")
+ return [hf.split("/")[-1][:-2] for hf in header_files]
+
+
+def _string_to_numpy_array(s):
+ return np.array([int(v) for v in s[1:-1].split(",")], dtype=np.int64)
+
+
+def _line_to_numpy_arrays(line):
+ lo_corner = _string_to_numpy_array(line[0][1:])
+ hi_corner = _string_to_numpy_array(line[1][:])
+ node_type = _string_to_numpy_array(line[2][:-1])
+ return lo_corner, hi_corner, node_type
+
+
+def _get_active_dimensions(box):
+ return box[1] - box[2] - box[0] + 1
+
+
+def _read_header(raw_file, field):
+ header_file = raw_file + field + "_H"
+ with open(header_file, "r") as f:
+
+ # skip the first five lines
+ for _ in range(5):
+ f.readline()
+
+ # read boxes
+ boxes = []
+ for line in f:
+ clean_line = line.strip().split()
+ if clean_line == [')']:
+ break
+ lo_corner, hi_corner, node_type = _line_to_numpy_arrays(clean_line)
+ boxes.append((lo_corner, hi_corner, node_type))
+
+ # read the file and offset position for the corresponding box
+ file_names = []
+ offsets = []
+ for line in f:
+ if line.startswith("FabOnDisk:"):
+ clean_line = line.strip().split()
+ file_names.append(clean_line[1])
+ offsets.append(int(clean_line[2]))
+
+ return boxes, file_names, offsets
+
+
class WarpXHierarchy(BoxlibHierarchy):
def __init__(self, ds, dataset_type="boxlib_native"):
@@ -1363,6 +1412,22 @@
line = f.readline()
species_id += 1
+ def _detect_output_fields(self):
+ super(WarpXHierarchy, self)._detect_output_fields()
+
+ # now detect the optional, non-cell-centered fields
+ self.raw_file = self.ds.output_dir + "/raw_fields/Level_0/"
+ self.ds.fluid_types += ("raw",)
+ self.raw_fields = _read_raw_field_names(self.raw_file)
+ self.field_list += [('raw', f) for f in self.raw_fields]
+ self.raw_field_map = {}
+ self.ds.nodal_flags = {}
+ for field_name in self.raw_fields:
+ boxes, file_names, offsets = _read_header(self.raw_file, field_name)
+ self.raw_field_map[field_name] = (boxes, file_names, offsets)
+ self.ds.nodal_flags[field_name] = np.array(boxes[0][2])
+
+
def _skip_line(line):
if len(line) == 0:
return True
diff -r bd65a61be0bd73785de1533681d3d35de3141d27 -r 3c2c470c05e9724cf5276c25c573ea4fe9840688 yt/frontends/boxlib/fields.py
--- a/yt/frontends/boxlib/fields.py
+++ b/yt/frontends/boxlib/fields.py
@@ -81,6 +81,14 @@
("", "particle_ones"),
)
+ def __init__(self, ds, field_list):
+ super(WarpXFieldInfo, self).__init__(ds, field_list)
+
+ # setup nodal flag information
+ for field in ds.index.raw_fields:
+ finfo = self.__getitem__(('raw', field))
+ finfo.nodal_flag = ds.nodal_flags[field]
+
def setup_particle_fields(self, ptype):
def get_mass(field, data):
diff -r bd65a61be0bd73785de1533681d3d35de3141d27 -r 3c2c470c05e9724cf5276c25c573ea4fe9840688 yt/frontends/boxlib/io.py
--- a/yt/frontends/boxlib/io.py
+++ b/yt/frontends/boxlib/io.py
@@ -22,6 +22,12 @@
from yt.funcs import mylog
from yt.frontends.chombo.io import parse_orion_sinks
+def _remove_raw(all_fields, raw_fields):
+ centered_fields = set(all_fields)
+ for raw in raw_fields:
+ centered_fields.discard(raw)
+ return list(centered_fields)
+
class IOHandlerBoxlib(BaseIOHandler):
_dataset_type = "boxlib_native"
@@ -31,25 +37,58 @@
def _read_fluid_selection(self, chunks, selector, fields, size):
chunks = list(chunks)
- if any((ftype != "boxlib" for ftype, fname in fields)):
+ if any(( not (ftype == "boxlib" or ftype == 'raw') for ftype, fname in fields)):
raise NotImplementedError
rv = {}
+ raw_fields = []
for field in fields:
- rv[field] = np.empty(size, dtype="float64")
+ if field[0] == "raw":
+ nodal_flag = self.ds.nodal_flags[field[1]]
+ num_nodes = 2**sum(nodal_flag)
+ rv[field] = np.empty((size, num_nodes), dtype="float64")
+ raw_fields.append(field)
+ else:
+ rv[field] = np.empty(size, dtype="float64")
+ centered_fields = _remove_raw(fields, raw_fields)
ng = sum(len(c.objs) for c in chunks)
mylog.debug("Reading %s cells of %s fields in %s grids",
size, [f2 for f1, f2 in fields], ng)
ind = 0
for chunk in chunks:
- data = self._read_chunk_data(chunk, fields)
+ data = self._read_chunk_data(chunk, centered_fields)
for g in chunk.objs:
for field in fields:
- ds = data[g.id].pop(field)
+ if field in centered_fields:
+ ds = data[g.id].pop(field)
+ else:
+ ds = self._read_raw_field(g, field)
nd = g.select(selector, ds, rv[field], ind) # caches
ind += nd
data.pop(g.id)
return rv
+ def _read_raw_field(self, grid, field):
+ field_name = field[1]
+ base_dir = self.ds.index.raw_file
+
+ box_list = self.ds.index.raw_field_map[field_name][0]
+ fn_list = self.ds.index.raw_field_map[field_name][1]
+ offset_list = self.ds.index.raw_field_map[field_name][2]
+
+ filename = base_dir + fn_list[grid.id]
+ offset = offset_list[grid.id]
+ box = box_list[grid.id]
+
+ lo = box[0]
+ hi = box[1]
+ shape = hi - lo + 1
+ with open(filename, "rb") as f:
+ f.seek(offset)
+ f.readline() # always skip the first line
+ arr = np.fromfile(f, 'float64', np.product(shape))
+ arr = arr.reshape(shape, order='F')
+ return arr
+
def _read_chunk_data(self, chunk, fields):
data = {}
grids_by_file = defaultdict(list)
https://bitbucket.org/yt_analysis/yt/commits/f6ef76633810/
Changeset: f6ef76633810
Branch: yt
User: atmyers
Date: 2017-04-05 19:36:30+00:00
Summary: allow derived fields to specify a nodal_flag.
Affected #: 1 file
diff -r 3c2c470c05e9724cf5276c25c573ea4fe9840688 -r f6ef76633810c48e0c27c84e46320be855fe1348 yt/fields/derived_field.py
--- a/yt/fields/derived_field.py
+++ b/yt/fields/derived_field.py
@@ -95,7 +95,7 @@
take_log=True, validators=None,
particle_type=None, vector_field=False, display_field=True,
not_in_all=False, display_name=None, output_units=None,
- dimensions=None, ds=None):
+ dimensions=None, ds=None, nodal_flag=None):
self.name = name
self.take_log = take_log
self.display_name = display_name
@@ -110,7 +110,10 @@
self.vector_field = vector_field
self.ds = ds
- self.nodal_flag = [0, 0, 0]
+ if nodal_flag is None:
+ self.nodal_flag = [0, 0, 0]
+ else:
+ self.nodal_flag = nodal_flag
self._function = function
https://bitbucket.org/yt_analysis/yt/commits/47590beae734/
Changeset: 47590beae734
Branch: yt
User: atmyers
Date: 2017-04-05 20:36:04+00:00
Summary: raise an exception indicating that projections don't work for nodal fields.
Affected #: 1 file
diff -r f6ef76633810c48e0c27c84e46320be855fe1348 -r 47590beae7341d7c1d59acf45debb96beeed8231 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -260,6 +260,11 @@
field = field or []
field = self._determine_fields(ensure_list(field))
+ for f in field:
+ nodal_flag = self.ds._get_field_info(f).nodal_flag
+ if any(nodal_flag):
+ raise RuntimeError("Nodal fields are currently not supported for projections.")
+
if not self.deserialize(field):
self.get_data(field)
self.serialize()
https://bitbucket.org/yt_analysis/yt/commits/15a434cbd3a4/
Changeset: 15a434cbd3a4
Branch: yt
User: atmyers
Date: 2017-04-10 16:48:21+00:00
Summary: remove duplicated line.
Affected #: 1 file
diff -r 47590beae7341d7c1d59acf45debb96beeed8231 -r 15a434cbd3a427d7366208826980f676cbad9e07 yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -379,7 +379,6 @@
return mask
def select(self, selector, source, dest, offset):
- nodal_flag = source.shape - self.ActiveDimensions
mask = self._get_selector_mask(selector)
count = self.count(selector)
if count == 0: return 0
https://bitbucket.org/yt_analysis/yt/commits/777cc2b74778/
Changeset: 777cc2b74778
Branch: yt
User: atmyers
Date: 2017-04-10 16:49:58+00:00
Summary: remove the now un-used count_nodal method of grid_patch
Affected #: 1 file
diff -r 15a434cbd3a427d7366208826980f676cbad9e07 -r 777cc2b7477869a2a3a10dcf6bb478e8b5782f6c yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -396,15 +396,6 @@
if mask is None: return 0
return self._last_count
- def count_nodal(self, selector, nodal_flag):
- # We don't cache the selector results
- mask = self._get_selector_mask(selector)
- if mask is None: return 0
- pad_width = [(0, nodal) for nodal in nodal_flag]
- mask = np.pad(mask, pad_width, 'edge')
- count = mask.sum()
- return count
-
def count_particles(self, selector, x, y, z):
# We don't cache the selector results
count = selector.count_points(x,y,z, 0.0)
https://bitbucket.org/yt_analysis/yt/commits/e222b91c3b86/
Changeset: e222b91c3b86
Branch: yt
User: atmyers
Date: 2017-04-10 16:54:18+00:00
Summary: remove nodal field logic from stream frontend.
Affected #: 1 file
diff -r 777cc2b7477869a2a3a10dcf6bb478e8b5782f6c -r e222b91c3b86dfe4484c7b44272b9434c4335d16 yt/frontends/stream/io.py
--- a/yt/frontends/stream/io.py
+++ b/yt/frontends/stream/io.py
@@ -48,13 +48,7 @@
raise NotImplementedError
rv = {}
for field in fields:
- nodal_flag = self.ds.nodal_flags[field[1]]
- num_nodes = 2**sum(nodal_flag)
- if num_nodes > 0:
- rv[field] = self.ds.arr(np.empty((size, num_nodes),
- dtype="float64"))
- else:
- rv[field] = self.ds.arr(np.empty(size, dtype="float64"))
+ rv[field] = self.ds.arr(np.empty(size, dtype="float64"))
ng = sum(len(c.objs) for c in chunks)
mylog.debug("Reading %s cells of %s fields in %s blocks",
https://bitbucket.org/yt_analysis/yt/commits/55ec50d706a3/
Changeset: 55ec50d706a3
Branch: yt
User: atmyers
Date: 2017-04-10 20:39:47+00:00
Summary: maintain the shape of cell-centered fields.
Affected #: 1 file
diff -r e222b91c3b86dfe4484c7b44272b9434c4335d16 -r 55ec50d706a3f1be5e4e06314db943f8fdd8668b yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -94,7 +94,7 @@
if not finfo.particle_type:
num_nodes = 2**sum(finfo.nodal_flag)
new_shape = list(self.ActiveDimensions) + [num_nodes]
- return tr.reshape(new_shape)
+ return np.squeeze(tr.reshape(new_shape))
return tr
def convert(self, datatype):
https://bitbucket.org/yt_analysis/yt/commits/1f2d7afb717b/
Changeset: 1f2d7afb717b
Branch: yt
User: atmyers
Date: 2017-04-16 20:44:48+00:00
Summary: adding a section to the docs about nodal fields.
Affected #: 1 file
diff -r 55ec50d706a3f1be5e4e06314db943f8fdd8668b -r 1f2d7afb717b0e0f00b73ef7e66ddef1f5213292 doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -325,6 +325,50 @@
``mach_number`` will always use the on-disk value, and not have yt
derive it, due to the complex interplay of the base state velocity.
+Viewing raw fields in WarpX
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Most BoxLib codes output cell-centered data. If the underlying discretization
+is not cell-centered, then fields are typically averaged to cell centers before
+they are written to plot files for visualization. WarpX, however, has the option
+to output the raw (i.e., not averaged to cell centers) data as well. If you
+run your WarpX simulation with ``warpx.plot_raw_fields = 1`` in your inputs
+file, then you should get an additional ``raw_fields`` subdirectory inside your
+plot file. When you load this dataset, yt will have additional on-disk fields
+defined, with the "raw" field type:
+
+.. code-block:: python
+ import yt
+ ds = yt.load("Laser/plt00015/")
+ print(ds.field_list)
+
+Instead of returning one field value per cell selected, raw fields return a
+number of values, depending on their centering. This centering is marked by
+a `nodal_flag` that describes whether the fields is nodal in each dimension.
+``nodal_flag = [0, 0, 0]`` means that the field is cell-centered, while
+``nodal_flag = [0, 0, 1]`` means that the field is nodal in the z direction
+and cell centered in the others, i.e. it is defined on the z faces of each cell.
+``nodal_flag = [1, 1, 0]`` would mean that the field is centered in the z direction,
+but nodal in the other two, i.e. it lives on the four cell edges that are normal
+to the z direction.
+
+..code-block:: python
+ ad = ds.all_data()
+ print(ds.field_info[('raw', 'Ex')].nodal_flag)
+ print(ad['raw', 'Ex'].shape)
+ print(ds.field_info[('raw', 'Bx')].nodal_flag)
+ print(ad['raw', 'Bx'].shape)
+ print(ds.field_info[('boxlib', 'Bx')].nodal_flag)
+ print(ad['boxlib', 'Bx'].shape)
+
+Here, the field ``('raw', 'Ex')`` is nodal in two directions, so four values per cell
+are returned, corresponding to the four edges in each cell on which the variable
+is defined. ``('raw', 'Bx')`` is nodal in one direction, so two values are returned
+per cell. The standard, averaged-to-cell-centers fields are still available.
+
+Currently, slices and data selection are implemented for nodal fields. Projections,
+volume rendering, and many of the analysis modules will not work.
+
.. _loading-pluto-data:
Pluto Data
https://bitbucket.org/yt_analysis/yt/commits/9cb5ede1a052/
Changeset: 9cb5ede1a052
Branch: yt
User: atmyers
Date: 2017-04-16 21:26:02+00:00
Summary: adding tests for raw fields.
Affected #: 2 files
diff -r 1f2d7afb717b0e0f00b73ef7e66ddef1f5213292 -r 9cb5ede1a0521348e130a2bbb7c5ace9faf6fa77 yt/frontends/boxlib/tests/test_outputs.py
--- a/yt/frontends/boxlib/tests/test_outputs.py
+++ b/yt/frontends/boxlib/tests/test_outputs.py
@@ -20,7 +20,9 @@
from yt.utilities.answer_testing.framework import \
requires_ds, \
small_patch_amr, \
- data_dir_load
+ data_dir_load, \
+ GridValuesTest, \
+ FieldValuesTest
from yt.frontends.boxlib.api import \
OrionDataset, \
NyxDataset, \
@@ -205,6 +207,17 @@
assert(np.all(np.logical_and(reg['particle_position_z'] <= right_edge[2],
reg['particle_position_z'] >= left_edge[2])))
+
+_raw_fields = [('raw', 'Bx'), ('raw', 'Ey'), ('raw', 'jz')]
+
+raw_fields = "Laser/plt00015"
+ at requires_ds(raw_fields)
+def test_raw_fields():
+ ds_fn = raw_fields
+ for field in _raw_fields:
+ yield GridValuesTest(ds_fn, field)
+
+
@requires_file(rt)
def test_OrionDataset():
assert isinstance(data_dir_load(rt), OrionDataset)
diff -r 1f2d7afb717b0e0f00b73ef7e66ddef1f5213292 -r 9cb5ede1a0521348e130a2bbb7c5ace9faf6fa77 yt/visualization/tests/test_raw_field_slices.py
--- /dev/null
+++ b/yt/visualization/tests/test_raw_field_slices.py
@@ -0,0 +1,60 @@
+"""
+Tests for making slices through raw fields
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2017, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+import os
+import os.path
+import tempfile
+import shutil
+import numpy as np
+import yt
+from yt.testing import fake_tetrahedral_ds
+from yt.testing import fake_hexahedral_ds
+from yt.utilities.answer_testing.framework import \
+ requires_ds, \
+ data_dir_load, \
+ GenericImageTest
+
+
+def setup():
+ """Test specific setup."""
+ from yt.config import ytcfg
+ ytcfg["yt", "__withintesting"] = "True"
+
+def compare(ds, field, test_prefix, decimals=12):
+ def slice_image(filename_prefix):
+ sl = yt.SlicePlot(ds, 'z', field)
+ sl.set_log('all', False)
+ image_file = sl.save(filename_prefix)
+ return image_file
+
+ slice_image.__name__ = "slice_{}".format(test_prefix)
+ test = GenericImageTest(ds, slice_image, decimals)
+ test.prefix = test_prefix
+ return test
+
+
+raw_fields = "Laser/plt00015"
+_raw_field_names = [('raw', 'Bx'),
+ ('raw', 'By'),
+ ('raw', 'Bz'),
+ ('raw', 'Ex'),
+ ('raw', 'Ey'),
+ ('raw', 'Ez'),
+ ('raw', 'jx'),
+ ('raw', 'jy'),
+ ('raw', 'jz')]
+ at requires_ds(raw_fields)
+def test_tri2():
+ ds = data_dir_load(raw_fields)
+ for field in _raw_field_names:
+ yield compare(ds, field, "answers_raw_%s" % field[1])
+
https://bitbucket.org/yt_analysis/yt/commits/9919c7d6d16d/
Changeset: 9919c7d6d16d
Branch: yt
User: atmyers
Date: 2017-04-16 22:38:30+00:00
Summary: remove unused imports
Affected #: 1 file
diff -r 9cb5ede1a0521348e130a2bbb7c5ace9faf6fa77 -r 9919c7d6d16d341f194e48b73072fb89a11cefe5 yt/visualization/tests/test_raw_field_slices.py
--- a/yt/visualization/tests/test_raw_field_slices.py
+++ b/yt/visualization/tests/test_raw_field_slices.py
@@ -10,14 +10,8 @@
#
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
-import os
-import os.path
-import tempfile
-import shutil
-import numpy as np
+
import yt
-from yt.testing import fake_tetrahedral_ds
-from yt.testing import fake_hexahedral_ds
from yt.utilities.answer_testing.framework import \
requires_ds, \
data_dir_load, \
https://bitbucket.org/yt_analysis/yt/commits/ad9ed208436d/
Changeset: ad9ed208436d
Branch: yt
User: atmyers
Date: 2017-04-16 23:31:50+00:00
Summary: remove unused import
Affected #: 1 file
diff -r 9919c7d6d16d341f194e48b73072fb89a11cefe5 -r ad9ed208436d3b74dd8c163cd2f898a495859b20 yt/frontends/boxlib/tests/test_outputs.py
--- a/yt/frontends/boxlib/tests/test_outputs.py
+++ b/yt/frontends/boxlib/tests/test_outputs.py
@@ -21,8 +21,7 @@
requires_ds, \
small_patch_amr, \
data_dir_load, \
- GridValuesTest, \
- FieldValuesTest
+ GridValuesTest
from yt.frontends.boxlib.api import \
OrionDataset, \
NyxDataset, \
https://bitbucket.org/yt_analysis/yt/commits/d522b72637c3/
Changeset: d522b72637c3
Branch: yt
User: atmyers
Date: 2017-04-19 18:27:39+00:00
Summary: Hopefully clarifying the nodal field docs
Affected #: 1 file
diff -r ad9ed208436d3b74dd8c163cd2f898a495859b20 -r d522b72637c3fdde92030c3a40294b8b37b20fac doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -342,7 +342,10 @@
ds = yt.load("Laser/plt00015/")
print(ds.field_list)
-Instead of returning one field value per cell selected, raw fields return a
+The raw fields in WarpX are nodal in at least one direction. We define a field
+to be "nodal" in a given direction if the field data is defined at the "low"
+and "high" sides of the cell in that direction, rather than at the cell center.
+Instead of returning one field value per cell selected, nodal fields return a
number of values, depending on their centering. This centering is marked by
a `nodal_flag` that describes whether the fields is nodal in each dimension.
``nodal_flag = [0, 0, 0]`` means that the field is cell-centered, while
https://bitbucket.org/yt_analysis/yt/commits/bbbbd8d43fa0/
Changeset: bbbbd8d43fa0
Branch: yt
User: atmyers
Date: 2017-04-19 18:32:50+00:00
Summary: Adding a docstring for the nodal_flag argument to derived_field.
Affected #: 1 file
diff -r d522b72637c3fdde92030c3a40294b8b37b20fac -r bbbbd8d43fa0dcd6da31876dec0c12e0f8a5c76a yt/fields/derived_field.py
--- a/yt/fields/derived_field.py
+++ b/yt/fields/derived_field.py
@@ -90,6 +90,15 @@
dimensions : str or object from yt.units.dimensions
The dimensions of the field, only needed if units="auto" and only used
for error checking.
+ nodal_flag : array-like with three components
+ This describes how the field is centered within a cell. If nodal_flag
+ is [0, 0, 0], then the field is cell-centered. If any of the components
+ of nodal_flag are 1, then the field is nodal in that direction, meaning
+ it is defined at the lo and hi sides of the cell rather than at the center.
+ For example, a field with nodal_flag = [1, 0, 0] would be defined at the
+ middle of the 2 x-faces of each cell. nodal_flag = [0, 1, 1] would mean the
+ that the field defined at the centers of the 4 edges that are normal to the
+ x axis, while nodal_flag = [1, 1, 1] would be defined at the 8 cell corners.
"""
def __init__(self, name, sampling_type, function, units=None,
take_log=True, validators=None,
https://bitbucket.org/yt_analysis/yt/commits/37ff7ec366ee/
Changeset: 37ff7ec366ee
Branch: yt
User: atmyers
Date: 2017-04-19 18:32:58+00:00
Summary: fix test name.
Affected #: 1 file
diff -r bbbbd8d43fa0dcd6da31876dec0c12e0f8a5c76a -r 37ff7ec366ee14773a581197430d0b100344fe1e yt/visualization/tests/test_raw_field_slices.py
--- a/yt/visualization/tests/test_raw_field_slices.py
+++ b/yt/visualization/tests/test_raw_field_slices.py
@@ -47,7 +47,7 @@
('raw', 'jy'),
('raw', 'jz')]
@requires_ds(raw_fields)
-def test_tri2():
+def test_raw_field_slices():
ds = data_dir_load(raw_fields)
for field in _raw_field_names:
yield compare(ds, field, "answers_raw_%s" % field[1])
https://bitbucket.org/yt_analysis/yt/commits/101488ef5657/
Changeset: 101488ef5657
Branch: yt
User: atmyers
Date: 2017-05-01 16:51:17+00:00
Summary: add new tests to tests.yaml
Affected #: 1 file
diff -r 37ff7ec366ee14773a581197430d0b100344fe1e -r 101488ef5657c9f943c4aa3107e3f71a4b2cfa3c tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -49,6 +49,7 @@
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_answers
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_filter
- yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers
+ - yt/visualization/tests/test_raw_field_slices.py:test_raw_field_slices
local_tipsy_002:
- yt/frontends/tipsy/tests/test_outputs.py
@@ -77,6 +78,7 @@
- yt/frontends/boxlib/tests/test_outputs.py:test_CastroDataset
- yt/frontends/boxlib/tests/test_outputs.py:test_RT_particles
- yt/frontends/boxlib/tests/test_outputs.py:test_units_override
+ - yt/frontends/boxlib/tests/test_outputs.py:test_raw_fields
local_boxlib_particles_003:
- yt/frontends/boxlib/tests/test_outputs.py:test_LyA
https://bitbucket.org/yt_analysis/yt/commits/94ba84bcd2fa/
Changeset: 94ba84bcd2fa
Branch: yt
User: atmyers
Date: 2017-05-01 17:51:53+00:00
Summary: Incrementing test numbers.
Affected #: 1 file
diff -r 101488ef5657c9f943c4aa3107e3f71a4b2cfa3c -r 94ba84bcd2fa9183d7f996986e4827f0fe96759a tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -42,7 +42,7 @@
local_owls_001:
- yt/frontends/owls/tests/test_outputs.py
- local_pw_013:
+ local_pw_014:
- yt/visualization/tests/test_plotwindow.py:test_attributes
- yt/visualization/tests/test_plotwindow.py:test_attributes_wt
- yt/visualization/tests/test_profile_plots.py:test_phase_plot_attributes
@@ -70,7 +70,7 @@
- yt/visualization/tests/test_mesh_slices.py:test_quad2
- yt/visualization/tests/test_mesh_slices.py:test_multi_region
- local_boxlib_003:
+ local_boxlib_004:
- yt/frontends/boxlib/tests/test_outputs.py:test_radadvect
- yt/frontends/boxlib/tests/test_outputs.py:test_radtube
- yt/frontends/boxlib/tests/test_outputs.py:test_star
https://bitbucket.org/yt_analysis/yt/commits/8f01bbe64d90/
Changeset: 8f01bbe64d90
Branch: yt
User: xarthisius
Date: 2017-05-01 19:38:19+00:00
Summary: Merged in atmyers/yt (pull request #2575)
Support face- and vertex- centered fields for WarpX datasets.
Approved-by: Nathan Goldbaum <ngoldbau at illinois.edu>
Approved-by: Michael Zingale <michael.zingale at stonybrook.edu>
Approved-by: yt-fido <yt.fido at gmail.com>
Approved-by: Kacper Kowalik <xarthisius.kk at gmail.com>
Affected #: 15 files
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -325,6 +325,53 @@
``mach_number`` will always use the on-disk value, and not have yt
derive it, due to the complex interplay of the base state velocity.
+Viewing raw fields in WarpX
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Most BoxLib codes output cell-centered data. If the underlying discretization
+is not cell-centered, then fields are typically averaged to cell centers before
+they are written to plot files for visualization. WarpX, however, has the option
+to output the raw (i.e., not averaged to cell centers) data as well. If you
+run your WarpX simulation with ``warpx.plot_raw_fields = 1`` in your inputs
+file, then you should get an additional ``raw_fields`` subdirectory inside your
+plot file. When you load this dataset, yt will have additional on-disk fields
+defined, with the "raw" field type:
+
+.. code-block:: python
+ import yt
+ ds = yt.load("Laser/plt00015/")
+ print(ds.field_list)
+
+The raw fields in WarpX are nodal in at least one direction. We define a field
+to be "nodal" in a given direction if the field data is defined at the "low"
+and "high" sides of the cell in that direction, rather than at the cell center.
+Instead of returning one field value per cell selected, nodal fields return a
+number of values, depending on their centering. This centering is marked by
+a `nodal_flag` that describes whether the fields is nodal in each dimension.
+``nodal_flag = [0, 0, 0]`` means that the field is cell-centered, while
+``nodal_flag = [0, 0, 1]`` means that the field is nodal in the z direction
+and cell centered in the others, i.e. it is defined on the z faces of each cell.
+``nodal_flag = [1, 1, 0]`` would mean that the field is centered in the z direction,
+but nodal in the other two, i.e. it lives on the four cell edges that are normal
+to the z direction.
+
+..code-block:: python
+ ad = ds.all_data()
+ print(ds.field_info[('raw', 'Ex')].nodal_flag)
+ print(ad['raw', 'Ex'].shape)
+ print(ds.field_info[('raw', 'Bx')].nodal_flag)
+ print(ad['raw', 'Bx'].shape)
+ print(ds.field_info[('boxlib', 'Bx')].nodal_flag)
+ print(ad['boxlib', 'Bx'].shape)
+
+Here, the field ``('raw', 'Ex')`` is nodal in two directions, so four values per cell
+are returned, corresponding to the four edges in each cell on which the variable
+is defined. ``('raw', 'Bx')`` is nodal in one direction, so two values are returned
+per cell. The standard, averaged-to-cell-centers fields are still available.
+
+Currently, slices and data selection are implemented for nodal fields. Projections,
+volume rendering, and many of the analysis modules will not work.
+
.. _loading-pluto-data:
Pluto Data
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -42,13 +42,14 @@
local_owls_001:
- yt/frontends/owls/tests/test_outputs.py
- local_pw_013:
+ local_pw_014:
- yt/visualization/tests/test_plotwindow.py:test_attributes
- yt/visualization/tests/test_plotwindow.py:test_attributes_wt
- yt/visualization/tests/test_profile_plots.py:test_phase_plot_attributes
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_answers
- yt/visualization/tests/test_particle_plot.py:test_particle_projection_filter
- yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers
+ - yt/visualization/tests/test_raw_field_slices.py:test_raw_field_slices
local_tipsy_002:
- yt/frontends/tipsy/tests/test_outputs.py
@@ -69,7 +70,7 @@
- yt/visualization/tests/test_mesh_slices.py:test_quad2
- yt/visualization/tests/test_mesh_slices.py:test_multi_region
- local_boxlib_003:
+ local_boxlib_004:
- yt/frontends/boxlib/tests/test_outputs.py:test_radadvect
- yt/frontends/boxlib/tests/test_outputs.py:test_radtube
- yt/frontends/boxlib/tests/test_outputs.py:test_star
@@ -77,6 +78,7 @@
- yt/frontends/boxlib/tests/test_outputs.py:test_CastroDataset
- yt/frontends/boxlib/tests/test_outputs.py:test_RT_particles
- yt/frontends/boxlib/tests/test_outputs.py:test_units_override
+ - yt/frontends/boxlib/tests/test_outputs.py:test_raw_fields
local_boxlib_particles_003:
- yt/frontends/boxlib/tests/test_outputs.py:test_LyA
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -260,6 +260,11 @@
field = field or []
field = self._determine_fields(ensure_list(field))
+ for f in field:
+ nodal_flag = self.ds._get_field_info(f).nodal_flag
+ if any(nodal_flag):
+ raise RuntimeError("Nodal fields are currently not supported for projections.")
+
if not self.deserialize(field):
self.get_data(field)
self.serialize()
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -29,6 +29,8 @@
YTParticleDepositionNotImplemented
from yt.utilities.lib.interpolators import \
ghost_zone_interpolate
+from yt.utilities.nodal_data_utils import \
+ get_nodal_slices
class AMRGridPatch(YTSelectionContainer):
_spatial = True
@@ -90,7 +92,9 @@
return tr
finfo = self.ds._get_field_info(*fields[0])
if not finfo.particle_type:
- return tr.reshape(self.ActiveDimensions)
+ num_nodes = 2**sum(finfo.nodal_flag)
+ new_shape = list(self.ActiveDimensions) + [num_nodes]
+ return np.squeeze(tr.reshape(new_shape))
return tr
def convert(self, datatype):
@@ -378,7 +382,13 @@
mask = self._get_selector_mask(selector)
count = self.count(selector)
if count == 0: return 0
- dest[offset:offset+count] = source[mask]
+ nodal_flag = source.shape - self.ActiveDimensions
+ if sum(nodal_flag) == 0:
+ dest[offset:offset+count] = source[mask]
+ else:
+ slices = get_nodal_slices(source.shape, nodal_flag)
+ for i , sl in enumerate(slices):
+ dest[offset:offset+count, i] = source[sl][mask]
return count
def count(self, selector):
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -268,7 +268,7 @@
_top_node = "/Slices"
_type_name = "slice"
_con_args = ('axis', 'coord')
- _container_fields = ("px", "py", "pdx", "pdy")
+ _container_fields = ("px", "py", "pz", "pdx", "pdy", "pdz")
def __init__(self, axis, coord, center=None, ds=None,
field_parameters=None, data_source=None):
YTSelectionContainer2D.__init__(self, axis, ds,
@@ -285,10 +285,14 @@
return self._current_chunk.fcoords[:,xax]
elif field == "py":
return self._current_chunk.fcoords[:,yax]
+ elif field == "pz":
+ return self._current_chunk.fcoords[:,self.axis]
elif field == "pdx":
return self._current_chunk.fwidth[:,xax] * 0.5
elif field == "pdy":
return self._current_chunk.fwidth[:,yax] * 0.5
+ elif field == "pdz":
+ return self._current_chunk.fwidth[:,self.axis] * 0.5
else:
raise KeyError(field)
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/fields/derived_field.py
--- a/yt/fields/derived_field.py
+++ b/yt/fields/derived_field.py
@@ -90,12 +90,21 @@
dimensions : str or object from yt.units.dimensions
The dimensions of the field, only needed if units="auto" and only used
for error checking.
+ nodal_flag : array-like with three components
+ This describes how the field is centered within a cell. If nodal_flag
+ is [0, 0, 0], then the field is cell-centered. If any of the components
+ of nodal_flag are 1, then the field is nodal in that direction, meaning
+ it is defined at the lo and hi sides of the cell rather than at the center.
+ For example, a field with nodal_flag = [1, 0, 0] would be defined at the
+ middle of the 2 x-faces of each cell. nodal_flag = [0, 1, 1] would mean the
+ that the field defined at the centers of the 4 edges that are normal to the
+ x axis, while nodal_flag = [1, 1, 1] would be defined at the 8 cell corners.
"""
def __init__(self, name, sampling_type, function, units=None,
take_log=True, validators=None,
particle_type=None, vector_field=False, display_field=True,
not_in_all=False, display_name=None, output_units=None,
- dimensions=None, ds=None):
+ dimensions=None, ds=None, nodal_flag=None):
self.name = name
self.take_log = take_log
self.display_name = display_name
@@ -110,6 +119,11 @@
self.vector_field = vector_field
self.ds = ds
+ if nodal_flag is None:
+ self.nodal_flag = [0, 0, 0]
+ else:
+ self.nodal_flag = nodal_flag
+
self._function = function
if validators:
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/frontends/boxlib/data_structures.py
--- a/yt/frontends/boxlib/data_structures.py
+++ b/yt/frontends/boxlib/data_structures.py
@@ -1333,6 +1333,55 @@
return vals
+def _read_raw_field_names(raw_file):
+ header_files = glob.glob(raw_file + "*_H")
+ return [hf.split("/")[-1][:-2] for hf in header_files]
+
+
+def _string_to_numpy_array(s):
+ return np.array([int(v) for v in s[1:-1].split(",")], dtype=np.int64)
+
+
+def _line_to_numpy_arrays(line):
+ lo_corner = _string_to_numpy_array(line[0][1:])
+ hi_corner = _string_to_numpy_array(line[1][:])
+ node_type = _string_to_numpy_array(line[2][:-1])
+ return lo_corner, hi_corner, node_type
+
+
+def _get_active_dimensions(box):
+ return box[1] - box[2] - box[0] + 1
+
+
+def _read_header(raw_file, field):
+ header_file = raw_file + field + "_H"
+ with open(header_file, "r") as f:
+
+ # skip the first five lines
+ for _ in range(5):
+ f.readline()
+
+ # read boxes
+ boxes = []
+ for line in f:
+ clean_line = line.strip().split()
+ if clean_line == [')']:
+ break
+ lo_corner, hi_corner, node_type = _line_to_numpy_arrays(clean_line)
+ boxes.append((lo_corner, hi_corner, node_type))
+
+ # read the file and offset position for the corresponding box
+ file_names = []
+ offsets = []
+ for line in f:
+ if line.startswith("FabOnDisk:"):
+ clean_line = line.strip().split()
+ file_names.append(clean_line[1])
+ offsets.append(int(clean_line[2]))
+
+ return boxes, file_names, offsets
+
+
class WarpXHierarchy(BoxlibHierarchy):
def __init__(self, ds, dataset_type="boxlib_native"):
@@ -1364,6 +1413,22 @@
line = f.readline()
species_id += 1
+ def _detect_output_fields(self):
+ super(WarpXHierarchy, self)._detect_output_fields()
+
+ # now detect the optional, non-cell-centered fields
+ self.raw_file = self.ds.output_dir + "/raw_fields/Level_0/"
+ self.ds.fluid_types += ("raw",)
+ self.raw_fields = _read_raw_field_names(self.raw_file)
+ self.field_list += [('raw', f) for f in self.raw_fields]
+ self.raw_field_map = {}
+ self.ds.nodal_flags = {}
+ for field_name in self.raw_fields:
+ boxes, file_names, offsets = _read_header(self.raw_file, field_name)
+ self.raw_field_map[field_name] = (boxes, file_names, offsets)
+ self.ds.nodal_flags[field_name] = np.array(boxes[0][2])
+
+
def _skip_line(line):
if len(line) == 0:
return True
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/frontends/boxlib/fields.py
--- a/yt/frontends/boxlib/fields.py
+++ b/yt/frontends/boxlib/fields.py
@@ -81,6 +81,14 @@
("", "particle_ones"),
)
+ def __init__(self, ds, field_list):
+ super(WarpXFieldInfo, self).__init__(ds, field_list)
+
+ # setup nodal flag information
+ for field in ds.index.raw_fields:
+ finfo = self.__getitem__(('raw', field))
+ finfo.nodal_flag = ds.nodal_flags[field]
+
def setup_particle_fields(self, ptype):
def get_mass(field, data):
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/frontends/boxlib/io.py
--- a/yt/frontends/boxlib/io.py
+++ b/yt/frontends/boxlib/io.py
@@ -22,6 +22,12 @@
from yt.funcs import mylog
from yt.frontends.chombo.io import parse_orion_sinks
+def _remove_raw(all_fields, raw_fields):
+ centered_fields = set(all_fields)
+ for raw in raw_fields:
+ centered_fields.discard(raw)
+ return list(centered_fields)
+
class IOHandlerBoxlib(BaseIOHandler):
_dataset_type = "boxlib_native"
@@ -31,25 +37,58 @@
def _read_fluid_selection(self, chunks, selector, fields, size):
chunks = list(chunks)
- if any((ftype != "boxlib" for ftype, fname in fields)):
+ if any(( not (ftype == "boxlib" or ftype == 'raw') for ftype, fname in fields)):
raise NotImplementedError
rv = {}
+ raw_fields = []
for field in fields:
- rv[field] = np.empty(size, dtype="float64")
+ if field[0] == "raw":
+ nodal_flag = self.ds.nodal_flags[field[1]]
+ num_nodes = 2**sum(nodal_flag)
+ rv[field] = np.empty((size, num_nodes), dtype="float64")
+ raw_fields.append(field)
+ else:
+ rv[field] = np.empty(size, dtype="float64")
+ centered_fields = _remove_raw(fields, raw_fields)
ng = sum(len(c.objs) for c in chunks)
mylog.debug("Reading %s cells of %s fields in %s grids",
size, [f2 for f1, f2 in fields], ng)
ind = 0
for chunk in chunks:
- data = self._read_chunk_data(chunk, fields)
+ data = self._read_chunk_data(chunk, centered_fields)
for g in chunk.objs:
for field in fields:
- ds = data[g.id].pop(field)
+ if field in centered_fields:
+ ds = data[g.id].pop(field)
+ else:
+ ds = self._read_raw_field(g, field)
nd = g.select(selector, ds, rv[field], ind) # caches
ind += nd
data.pop(g.id)
return rv
+ def _read_raw_field(self, grid, field):
+ field_name = field[1]
+ base_dir = self.ds.index.raw_file
+
+ box_list = self.ds.index.raw_field_map[field_name][0]
+ fn_list = self.ds.index.raw_field_map[field_name][1]
+ offset_list = self.ds.index.raw_field_map[field_name][2]
+
+ filename = base_dir + fn_list[grid.id]
+ offset = offset_list[grid.id]
+ box = box_list[grid.id]
+
+ lo = box[0]
+ hi = box[1]
+ shape = hi - lo + 1
+ with open(filename, "rb") as f:
+ f.seek(offset)
+ f.readline() # always skip the first line
+ arr = np.fromfile(f, 'float64', np.product(shape))
+ arr = arr.reshape(shape, order='F')
+ return arr
+
def _read_chunk_data(self, chunk, fields):
data = {}
grids_by_file = defaultdict(list)
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/frontends/boxlib/tests/test_outputs.py
--- a/yt/frontends/boxlib/tests/test_outputs.py
+++ b/yt/frontends/boxlib/tests/test_outputs.py
@@ -20,7 +20,8 @@
from yt.utilities.answer_testing.framework import \
requires_ds, \
small_patch_amr, \
- data_dir_load
+ data_dir_load, \
+ GridValuesTest
from yt.frontends.boxlib.api import \
OrionDataset, \
NyxDataset, \
@@ -205,6 +206,17 @@
assert(np.all(np.logical_and(reg['particle_position_z'] <= right_edge[2],
reg['particle_position_z'] >= left_edge[2])))
+
+_raw_fields = [('raw', 'Bx'), ('raw', 'Ey'), ('raw', 'jz')]
+
+raw_fields = "Laser/plt00015"
+ at requires_ds(raw_fields)
+def test_raw_fields():
+ ds_fn = raw_fields
+ for field in _raw_fields:
+ yield GridValuesTest(ds_fn, field)
+
+
@requires_file(rt)
def test_OrionDataset():
assert isinstance(data_dir_load(rt), OrionDataset)
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/frontends/stream/io.py
--- a/yt/frontends/stream/io.py
+++ b/yt/frontends/stream/io.py
@@ -49,6 +49,7 @@
rv = {}
for field in fields:
rv[field] = self.ds.arr(np.empty(size, dtype="float64"))
+
ng = sum(len(c.objs) for c in chunks)
mylog.debug("Reading %s cells of %s fields in %s blocks",
size, [f2 for f1, f2 in fields], ng)
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/geometry/coordinates/cartesian_coordinates.py
--- a/yt/geometry/coordinates/cartesian_coordinates.py
+++ b/yt/geometry/coordinates/cartesian_coordinates.py
@@ -24,8 +24,9 @@
from yt.funcs import mylog
from yt.utilities.lib.pixelization_routines import \
pixelize_element_mesh, pixelize_off_axis_cartesian, \
- pixelize_cartesian
+ pixelize_cartesian, pixelize_cartesian_nodal
from yt.data_objects.unstructured_mesh import SemiStructuredMesh
+from yt.utilities.nodal_data_utils import get_nodal_data
class CartesianCoordinateHandler(CoordinateHandler):
@@ -126,14 +127,27 @@
period[1] = self.period[self.y_axis[dim]]
if hasattr(period, 'in_units'):
period = period.in_units("code_length").d
+
buff = np.zeros((size[1], size[0]), dtype="f8")
- pixelize_cartesian(buff, data_source['px'], data_source['py'],
- data_source['pdx'], data_source['pdy'],
- data_source[field],
- bounds, int(antialias),
- period, int(periodic))
+
+ finfo = self.ds._get_field_info(field)
+ nodal_flag = finfo.nodal_flag
+ if np.any(nodal_flag):
+ nodal_data = get_nodal_data(data_source, field)
+ coord = data_source.coord.d
+ pixelize_cartesian_nodal(buff,
+ data_source['px'], data_source['py'], data_source['pz'],
+ data_source['pdx'], data_source['pdy'], data_source['pdz'],
+ nodal_data, coord, bounds, int(antialias),
+ period, int(periodic))
+ else:
+ pixelize_cartesian(buff, data_source['px'], data_source['py'],
+ data_source['pdx'], data_source['pdy'],
+ data_source[field],
+ bounds, int(antialias),
+ period, int(periodic))
return buff
-
+
def _oblique_pixelize(self, data_source, field, bounds, size, antialias):
indices = np.argsort(data_source['pdx'])[::-1]
buff = np.zeros((size[1], size[0]), dtype="f8")
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/utilities/lib/pixelization_routines.pyx
--- a/yt/utilities/lib/pixelization_routines.pyx
+++ b/yt/utilities/lib/pixelization_routines.pyx
@@ -233,6 +233,171 @@
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
+def pixelize_cartesian_nodal(np.float64_t[:,:] buff,
+ np.float64_t[:] px,
+ np.float64_t[:] py,
+ np.float64_t[:] pz,
+ np.float64_t[:] pdx,
+ np.float64_t[:] pdy,
+ np.float64_t[:] pdz,
+ np.float64_t[:, :] data,
+ np.float64_t coord,
+ bounds,
+ int antialias = 1,
+ period = None,
+ int check_period = 1):
+ cdef np.float64_t x_min, x_max, y_min, y_max
+ cdef np.float64_t period_x = 0.0, period_y = 0.0
+ cdef np.float64_t width, height, px_dx, px_dy, ipx_dx, ipx_dy
+ cdef np.float64_t ld_x, ld_y, cx, cy, cz
+ cdef int i, j, p, xi, yi
+ cdef int lc, lr, rc, rr
+ cdef np.float64_t lypx, rypx, lxpx, rxpx, overlap1, overlap2
+ # These are the temp vars we get from the arrays
+ cdef np.float64_t oxsp, oysp, ozsp
+ cdef np.float64_t xsp, ysp, zsp
+ cdef np.float64_t dxsp, dysp, dzsp
+ # Some periodicity helpers
+ cdef int xiter[2]
+ cdef int yiter[2]
+ cdef int ii, jj, kk, ind
+ cdef np.float64_t xiterv[2]
+ cdef np.float64_t yiterv[2]
+ if period is not None:
+ period_x = period[0]
+ period_y = period[1]
+ x_min = bounds[0]
+ x_max = bounds[1]
+ y_min = bounds[2]
+ y_max = bounds[3]
+ width = x_max - x_min
+ height = y_max - y_min
+ px_dx = width / (<np.float64_t> buff.shape[1])
+ px_dy = height / (<np.float64_t> buff.shape[0])
+ ipx_dx = 1.0 / px_dx
+ ipx_dy = 1.0 / px_dy
+ if px.shape[0] != py.shape[0] or \
+ px.shape[0] != pz.shape[0] or \
+ px.shape[0] != pdx.shape[0] or \
+ px.shape[0] != pdy.shape[0] or \
+ px.shape[0] != pdz.shape[0] or \
+ px.shape[0] != data.shape[0]:
+ raise YTPixelizeError("Arrays are not of correct shape.")
+ xiter[0] = yiter[0] = 0
+ xiterv[0] = yiterv[0] = 0.0
+ # Here's a basic outline of what we're going to do here. The xiter and
+ # yiter variables govern whether or not we should check periodicity -- are
+ # we both close enough to the edge that it would be important *and* are we
+ # periodic?
+ #
+ # The other variables are all either pixel positions or data positions.
+ # Pixel positions will vary regularly from the left edge of the window to
+ # the right edge of the window; px_dx and px_dy are the dx (cell width, not
+ # half-width). ipx_dx and ipx_dy are the inverse, for quick math.
+ #
+ # The values in xsp, dxsp, x_min and their y counterparts, are the
+ # data-space coordinates, and are related to the data fed in. We make some
+ # modifications for periodicity.
+ #
+ # Inside the finest loop, we compute the "left column" (lc) and "lower row"
+ # (lr) and then iterate up to "right column" (rc) and "uppeR row" (rr),
+ # depositing into them the data value. Overlap computes the relative
+ # overlap of a data value with a pixel.
+ #
+ # NOTE ON ROWS AND COLUMNS:
+ #
+ # The way that images are plotting in matplotlib is somewhat different
+ # from what most might expect. The first axis of the array plotted is
+ # what varies along the x axis. So for instance, if you supply
+ # origin='lower' and plot the results of an mgrid operation, at a fixed
+ # 'y' value you will see the results of that array held constant in the
+ # first dimension. Here is some example code:
+ #
+ # import matplotlib.pyplot as plt
+ # import numpy as np
+ # x, y = np.mgrid[0:1:100j,0:1:100j]
+ # plt.imshow(x, interpolation='nearest', origin='lower')
+ # plt.imshow(y, interpolation='nearest', origin='lower')
+ #
+ # The values in the image:
+ # lower left: arr[0,0]
+ # lower right: arr[0,-1]
+ # upper left: arr[-1,0]
+ # upper right: arr[-1,-1]
+ #
+ # So what we want here is to fill an array such that we fill:
+ # first axis : y_min .. y_max
+ # second axis: x_min .. x_max
+ with nogil:
+ for p in range(px.shape[0]):
+ xiter[1] = yiter[1] = 999
+ oxsp = px[p]
+ oysp = py[p]
+ ozsp = pz[p]
+ dxsp = pdx[p]
+ dysp = pdy[p]
+ dzsp = pdz[p]
+ if check_period == 1:
+ if (oxsp - dxsp < x_min):
+ xiter[1] = +1
+ xiterv[1] = period_x
+ elif (oxsp + dxsp > x_max):
+ xiter[1] = -1
+ xiterv[1] = -period_x
+ if (oysp - dysp < y_min):
+ yiter[1] = +1
+ yiterv[1] = period_y
+ elif (oysp + dysp > y_max):
+ yiter[1] = -1
+ yiterv[1] = -period_y
+ overlap1 = overlap2 = 1.0
+ zsp = ozsp
+ for xi in range(2):
+ if xiter[xi] == 999: continue
+ xsp = oxsp + xiterv[xi]
+ if (xsp + dxsp < x_min) or (xsp - dxsp > x_max): continue
+ for yi in range(2):
+ if yiter[yi] == 999: continue
+ ysp = oysp + yiterv[yi]
+ if (ysp + dysp < y_min) or (ysp - dysp > y_max): continue
+ lc = <int> fmax(((xsp-dxsp-x_min)*ipx_dx),0)
+ lr = <int> fmax(((ysp-dysp-y_min)*ipx_dy),0)
+ # NOTE: This is a different way of doing it than in the C
+ # routines. In C, we were implicitly casting the
+ # initialization to int, but *not* the conditional, which
+ # was allowed an extra value:
+ # for(j=lc;j<rc;j++)
+ # here, when assigning lc (double) to j (int) it got
+ # truncated, but no similar truncation was done in the
+ # comparison of j to rc (double). So give ourselves a
+ # bonus row and bonus column here.
+ rc = <int> fmin(((xsp+dxsp-x_min)*ipx_dx + 1), buff.shape[1])
+ rr = <int> fmin(((ysp+dysp-y_min)*ipx_dy + 1), buff.shape[0])
+ # Note that we're iterating here over *y* in the i
+ # direction. See the note above about this.
+ for i in range(lr, rr):
+ lypx = px_dy * i + y_min
+ rypx = px_dy * (i+1) + y_min
+ for j in range(lc, rc):
+ lxpx = px_dx * j + x_min
+ rxpx = px_dx * (j+1) + x_min
+
+ cx = (rxpx+lxpx)*0.5
+ cy = (rypx+lypx)*0.5
+ cz = coord
+
+ ii = <int> (cx - xsp + dxsp)
+ jj = <int> (cy - ysp + dysp)
+ kk = <int> (cz - zsp + dzsp)
+
+ ind = 4*ii + 2*jj + kk
+
+ buff[i,j] = data[p, ind]
+
+
+ at cython.cdivision(True)
+ at cython.boundscheck(False)
+ at cython.wraparound(False)
def pixelize_off_axis_cartesian(
np.float64_t[:,:] buff,
np.float64_t[:] x,
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/utilities/nodal_data_utils.py
--- /dev/null
+++ b/yt/utilities/nodal_data_utils.py
@@ -0,0 +1,41 @@
+import numpy as np
+
+_index_map = np.array([[0, 0, 0, 0, 0, 0, 0, 0],
+ [0, 1, 0, 1, 0, 1, 0, 1],
+ [0, 0, 1, 1, 0, 0, 1, 1],
+ [0, 1, 2, 3, 0, 1, 2, 3],
+ [0, 0, 0, 0, 1, 1, 1, 1],
+ [0, 1, 0, 1, 2, 3, 2, 3],
+ [0, 0, 1, 1, 2, 2, 3, 3],
+ [0, 1, 2, 3, 4, 5, 6, 7]])
+
+def _get_linear_index(nodal_flag):
+ return 1*nodal_flag[2] + 2*nodal_flag[1] + 4*nodal_flag[0]
+
+def _get_indices(nodal_flag):
+ li = _get_linear_index(nodal_flag)
+ return _index_map[li]
+
+def get_nodal_data(data_source, field):
+ finfo = data_source.ds._get_field_info(field)
+ nodal_flag = finfo.nodal_flag
+ field_data = data_source[field]
+ inds = _get_indices(nodal_flag)
+ return field_data[:, inds]
+
+def get_nodal_slices(shape, nodal_flag):
+ slices = []
+ dir_slices = [[] for _ in range(3)]
+
+ for i in range(3):
+ if nodal_flag[i]:
+ dir_slices[i] = [slice(0, shape[i]-1), slice(1, shape[i])]
+ else:
+ dir_slices[i] = [slice(0, shape[i])]
+
+ for i, sl_i in enumerate(dir_slices[0]):
+ for j, sl_j in enumerate(dir_slices[1]):
+ for k, sl_k in enumerate(dir_slices[2]):
+ slices.append([sl_i, sl_j, sl_k])
+
+ return slices
diff -r 09d6e114abd2fed4d527359f6c9ed87a59a1e1ac -r 8f01bbe64d90b35ea7d65e2918791cdbec362ccb yt/visualization/tests/test_raw_field_slices.py
--- /dev/null
+++ b/yt/visualization/tests/test_raw_field_slices.py
@@ -0,0 +1,54 @@
+"""
+Tests for making slices through raw fields
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2017, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import yt
+from yt.utilities.answer_testing.framework import \
+ requires_ds, \
+ data_dir_load, \
+ GenericImageTest
+
+
+def setup():
+ """Test specific setup."""
+ from yt.config import ytcfg
+ ytcfg["yt", "__withintesting"] = "True"
+
+def compare(ds, field, test_prefix, decimals=12):
+ def slice_image(filename_prefix):
+ sl = yt.SlicePlot(ds, 'z', field)
+ sl.set_log('all', False)
+ image_file = sl.save(filename_prefix)
+ return image_file
+
+ slice_image.__name__ = "slice_{}".format(test_prefix)
+ test = GenericImageTest(ds, slice_image, decimals)
+ test.prefix = test_prefix
+ return test
+
+
+raw_fields = "Laser/plt00015"
+_raw_field_names = [('raw', 'Bx'),
+ ('raw', 'By'),
+ ('raw', 'Bz'),
+ ('raw', 'Ex'),
+ ('raw', 'Ey'),
+ ('raw', 'Ez'),
+ ('raw', 'jx'),
+ ('raw', 'jy'),
+ ('raw', 'jz')]
+ at requires_ds(raw_fields)
+def test_raw_field_slices():
+ ds = data_dir_load(raw_fields)
+ for field in _raw_field_names:
+ yield compare(ds, field, "answers_raw_%s" % field[1])
+
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list