[Yt-svn] yt-commit r977 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Wed Dec 3 19:03:06 PST 2008
Author: mturk
Date: Wed Dec 3 19:03:06 2008
New Revision: 977
URL: http://yt.spacepope.org/changeset/977
Log:
Speedups to the Covering Grids -- but not the smoothed covering grids, which I
think largely need to be either replaced or reworked.
(And Jeff said PointCombine.c was 'read-only'!)
Now the exceptions for getting ghost zones ACTUALLY pass upwards the fields
they need; furthermore, the refinement ACTUALLY sets appropriate min/max bounds
rather than just iterating over the entire grid every time.
The grids iterated over for covering grids start at the highest level and work
backwards now.
The fields are now grabbed several-at-a-time instead of one-at-a-time.
I think that's all...
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/FieldInfoContainer.py
trunk/yt/lagos/PointCombine.c
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Wed Dec 3 19:03:06 2008
@@ -1647,7 +1647,7 @@
self.left_edge, self.right_edge)
level_ind = na.where(self.pf.hierarchy.gridLevels.ravel()[ind] <= self.level)
sort_ind = na.argsort(self.pf.h.gridLevels.ravel()[ind][level_ind])
- self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)]
+ self._grids = self.pf.hierarchy.grids[ind][level_ind][(sort_ind,)][::-1]
def extract_region(self, indices):
mylog.error("Sorry, dude, do it yourself, it's already in 3-D.")
@@ -1662,27 +1662,29 @@
self._get_list_of_grids()
# We don't generate coordinates here.
if field == None:
- fields_to_get = self.fields
+ _fields_to_get = self.fields
else:
- fields_to_get = ensure_list(field)
+ _fields_to_get = ensure_list(field)
+ fields_to_get = [f for f in _fields_to_get if f not in self.data]
+ if len(fields_to_get) == 0: return
for field in fields_to_get:
- if self.data.has_key(field):
- continue
- mylog.debug("Getting field %s from %s possible grids",
- field, len(self._grids))
self[field] = na.zeros(self.ActiveDimensions, dtype='float64') -999
- if self._use_pbar: pbar = \
- get_pbar('Searching grids for values ', len(self._grids))
- for i,grid in enumerate(self._grids):
- if self._use_pbar: pbar.update(i)
- self._get_data_from_grid(grid, field)
- if not na.any(self[field] == -999): break
- if self._use_pbar: pbar.finish()
- if na.any(self[field] == -999):# and self.dx < self.hierarchy.grids[0].dx:
- print "COVERING PROBLEM", na.where(self[field]==-999)[0].size
- print na.where(self[field]==-999)
- return
- raise KeyError
+ mylog.debug("Getting fields %s from %s possible grids",
+ field, len(self._grids))
+ if self._use_pbar: pbar = \
+ get_pbar('Searching grids for values ', len(self._grids))
+ field = fields_to_get[-1]
+ for i,grid in enumerate(self._grids):
+ if self._use_pbar: pbar.update(i)
+ self._get_data_from_grid(grid, fields_to_get)
+ if not na.any(self[field] == -999): break
+ if self._use_pbar: pbar.finish()
+ if na.any(self[field] == -999):
+ # and self.dx < self.hierarchy.grids[0].dx:
+ print "COVERING PROBLEM", na.where(self[field]==-999)[0].size
+ print na.where(self[fields_to_get[0]]==-999)
+ return
+ raise KeyError
def flush_data(self, field=None):
"""
Modified: trunk/yt/lagos/FieldInfoContainer.py
==============================================================================
--- trunk/yt/lagos/FieldInfoContainer.py (original)
+++ trunk/yt/lagos/FieldInfoContainer.py Wed Dec 3 19:03:06 2008
@@ -281,4 +281,4 @@
raise NeedsGridType(self.ghost_zones,self.fields)
if self.ghost_zones == data._num_ghost_zones:
return True
- raise NeedsGridType(self.ghost_zones)
+ raise NeedsGridType(self.ghost_zones,self.fields)
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Wed Dec 3 19:03:06 2008
@@ -671,7 +671,7 @@
tc_data = (PyObject *) PyList_GetItem(oc_data, n);
c_data[n] = (PyArrayObject *) PyArray_FromAny(tc_data,
PyArray_DescrFromType(NPY_FLOAT64), 3, 3,
- NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ NPY_UPDATEIFCOPY, NULL);
if((c_data[n]==NULL) || (c_data[n]->nd != 3)) {
PyErr_Format(_dataCubeError,
"CombineGrids: Three dimensions required for c_data[%i].",n);
@@ -691,10 +691,13 @@
malloc(sizeof(PyArrayObject*)*n_fields);
for (n=0;n<n_fields;n++)g_data[n]=NULL;
for (n=0;n<n_fields;n++){
+ /* Borrows a reference */
tg_data = (PyObject *) PyList_GetItem(og_data, n);
+ /* We set up an array so we only have to do this once
+ Note that this is an array in the C sense, not the NumPy sense */
g_data[n] = (PyArrayObject *) PyArray_FromAny(tg_data,
PyArray_DescrFromType(NPY_FLOAT64), 3, 3,
- NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ NPY_UPDATEIFCOPY, NULL);
if((g_data[n]==NULL) || (g_data[n]->nd != 3)) {
PyErr_Format(_dataCubeError,
"CombineGrids: Three dimensions required for g_data[%i].",n);
@@ -714,9 +717,11 @@
npy_float64 ac_le_p[3][3];
npy_float64 ac_re_p[3][3];
npy_float64 ag_re[3];
- for (i=0;i<3;i++){ag_re[i] = ag_le[i]+ag_dx[i]*(g_data[0]->dimensions[i]+1);}
+ /* This is for checking for periodic boundary conditions.
+ Manually set the right edge to be offset from the left. */
+ for(i=0;i<3;i++){ag_re[i] = ag_le[i]+ag_dx[i]*(g_data[0]->dimensions[i]+1);}
- for (i=0;i<3;i++){ac_le_p[i][0] = ac_le[i]; ac_re_p[i][0] = ac_re[i];}
+ for(i=0;i<3;i++){ac_le_p[i][0] = ac_le[i]; ac_re_p[i][0] = ac_re[i];}
for(i=0;i<3;i++) {
itc = 1;
if (ac_le[i] < adl_edge[i]) {
@@ -729,27 +734,49 @@
}
p_niter[i] = itc;
}
+ npy_intp nx, ny, nz;
+ /* This is easier than doing a lookup every loop */
+ nx = PyArray_DIM(c_data[0], 0);
+ ny = PyArray_DIM(c_data[0], 1);
+ nz = PyArray_DIM(c_data[0], 2);
+ npy_int64 xg_min, yg_min, zg_min;
+ npy_int64 xg_max, yg_max, zg_max;
- for (xg = 0; xg < g_data[0]->dimensions[0]; xg++) {
+ /* Periodic iterations, *if necessary* */
for (pxl = 0; pxl < p_niter[0]; pxl++) {
+ xg_min = max(floor((ac_le_p[0][pxl]-ag_le[0])/ag_dx[0]) - 1, 0);
+ xg_max = min(ceil((ac_re_p[0][pxl]-ag_le[0])/ag_dx[0]) + 1,
+ g_data[0]->dimensions[0]);
+ for (xg = xg_min; xg < xg_max; xg++) {
+ /* If we're off the destination cell boundary, skip */
if (ag_le[0]+ag_dx[0]*xg > ac_re_p[0][pxl]) continue;
if (ag_le[0]+ag_dx[0]*(xg+1) < ac_le_p[0][pxl]) continue;
+ /* Floor to the source edge */
cmin_x = max(floorl((ag_le[0]+ag_dx[0]*xg - ac_le_p[0][pxl])/ac_dx[0]),0);
- cmax_x = min( ceill((ag_le[0]+ag_dx[0]*(xg+1) - ac_le_p[0][pxl])/ac_dx[0]),PyArray_DIM(c_data[0],0));
- for (yg = 0; yg < g_data[0]->dimensions[1]; yg++) {
+ cmax_x = min( ceill((ag_le[0]+ag_dx[0]*(xg+1) - ac_le_p[0][pxl])/ac_dx[0]),nx);
+ if(cmin_x==cmax_x)continue;
for (pyl = 0; pyl < p_niter[1]; pyl++) {
+ yg_min = max(floor((ac_le_p[1][pyl]-ag_le[1])/ag_dx[1]) - 1, 0);
+ yg_max = min(ceil((ac_re_p[1][pyl]-ag_le[1])/ag_dx[1]),
+ g_data[0]->dimensions[1]);
+ for (yg = yg_min; yg < yg_max; yg++) {
if (ag_le[1]+ag_dx[1]*yg > ac_re_p[1][pyl]) continue;
if (ag_le[1]+ag_dx[1]*(yg+1) < ac_le_p[1][pyl]) continue;
cmin_y = max(floorl((ag_le[1]+ag_dx[1]*yg - ac_le_p[1][pyl])/ac_dx[1]),0);
- cmax_y = min( ceill((ag_le[1]+ag_dx[1]*(yg+1) - ac_le_p[1][pyl])/ac_dx[1]),PyArray_DIM(c_data[0],1));
- for (zg = 0; zg < g_data[0]->dimensions[2]; zg++) {
+ cmax_y = min( ceill((ag_le[1]+ag_dx[1]*(yg+1) - ac_le_p[1][pyl])/ac_dx[1]),ny);
+ if(cmin_y==cmax_y)continue;
for (pzl = 0; pzl < p_niter[2]; pzl++) {
- cm = *(npy_int *)PyArray_GETPTR3(g_cm,xg,yg,zg);
- if ((!ll) && (cm == 0)) continue;
+ zg_min = max(floor((ac_le_p[2][pzl]-ag_le[2])/ag_dx[2]) - 1, 0);
+ zg_max = min(ceil((ac_re_p[2][pzl]-ag_le[2])/ag_dx[2]),
+ g_data[0]->dimensions[2]);
+ for (zg = zg_min; zg < zg_max; zg++) {
+ cm = *(npy_int *)PyArray_GETPTR3(g_cm,xg,yg,zg);
+ if ((!ll) && (cm == 0)) continue;
if (ag_le[2]+ag_dx[2]*zg > ac_re_p[2][pzl]) continue;
if (ag_le[2]+ag_dx[2]*(zg+1) < ac_le_p[2][pzl]) continue;
cmin_z = max(floorl((ag_le[2]+ag_dx[2]*zg - ac_le_p[2][pzl])/ac_dx[2]),0);
- cmax_z = min( ceill((ag_le[2]+ag_dx[2]*(zg+1) - ac_le_p[2][pzl])/ac_dx[2]),PyArray_DIM(c_data[0],2));
+ cmax_z = min( ceill((ag_le[2]+ag_dx[2]*(zg+1) - ac_le_p[2][pzl])/ac_dx[2]),nz);
+ if(cmin_z==cmax_z)continue;
for (xc = cmin_x; xc < cmax_x ; xc++) {
for (yc = cmin_y; yc < cmax_y ; yc++) {
for (zc = cmin_z; zc < cmax_z ; zc++) {
More information about the yt-svn
mailing list