[Yt-svn] yt-commit r589 - trunk/yt/lagos
mturk at wrangler.dreamhost.com
mturk at wrangler.dreamhost.com
Sat Jun 21 08:52:14 PDT 2008
Author: mturk
Date: Sat Jun 21 08:52:13 2008
New Revision: 589
URL: http://yt.spacepope.org/changeset/589
Log:
Added periodic boundary conditions for covering grids. (Actually, now, all
covering grids will assume the domain is periodic. This might be unwanted, so
perhaps it will be revisited in the future.)
Modified:
trunk/yt/lagos/BaseDataTypes.py
trunk/yt/lagos/HierarchyType.py
trunk/yt/lagos/PointCombine.c
Modified: trunk/yt/lagos/BaseDataTypes.py
==============================================================================
--- trunk/yt/lagos/BaseDataTypes.py (original)
+++ trunk/yt/lagos/BaseDataTypes.py Sat Jun 21 08:52:13 2008
@@ -1457,7 +1457,13 @@
self._refresh_data()
def _get_list_of_grids(self):
- grids, ind = self.pf.hierarchy.get_box_grids(self.left_edge, self.right_edge)
+ if na.any(self.left_edge < self.pf["DomainLeftEdge"]) or \
+ na.any(self.right_edge > self.pf["DomainRightEdge"]):
+ grids,ind = self.pf.hierarchy.get_periodic_box_grids(
+ self.left_edge.copy(), self.right_edge.copy())
+ else:
+ grids,ind = self.pf.hierarchy.get_box_grids(
+ 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,)]
@@ -1493,6 +1499,8 @@
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
def flush_data(self, field=None):
@@ -1523,7 +1531,7 @@
PointCombine.DataCubeRefine(
grid.LeftEdge, g_dx, g_fields, grid.child_mask,
self.left_edge, self.right_edge, c_dx, c_fields,
- ll)
+ ll, self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"])
def _flush_data_to_grid(self, grid, fields):
ll = int(grid.Level == self.level)
@@ -1538,7 +1546,7 @@
PointCombine.DataCubeReplace(
grid.LeftEdge, g_dx, g_fields, grid.child_mask,
self.left_edge, self.right_edge, c_dx, c_fields,
- ll)
+ ll, self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"])
class EnzoSmoothedCoveringGrid(EnzoCoveringGrid):
def __init__(self, *args, **kwargs):
@@ -1633,7 +1641,7 @@
grid.LeftEdge, g_dx, g_fields, grid.child_mask,
self.left_edge-c_dx, self.right_edge+c_dx,
c_dx, c_fields,
- 1)
+ 1, self.pf["DomainLeftEdge"], self.pf["DomainRightEdge"])
def flush_data(self, *args, **kwargs):
raise KeyError("Can't do this")
Modified: trunk/yt/lagos/HierarchyType.py
==============================================================================
--- trunk/yt/lagos/HierarchyType.py (original)
+++ trunk/yt/lagos/HierarchyType.py Sat Jun 21 08:52:13 2008
@@ -615,13 +615,31 @@
gridI = na.where(na.logical_and((self.gridDxs<=radius)[:,0],(dist < (radius + long_axis))) == 1)
return self.grids[gridI], gridI
- def get_box_grids(self, leftEdge, rightEdge):
+ def get_box_grids(self, left_edge, right_edge):
"""
Gets back all the grids between a left edge and right edge
"""
- gridI = na.where((na.all(self.gridRightEdge > leftEdge, axis=1)
- & na.all(self.gridLeftEdge < rightEdge, axis=1)) == True)
- return self.grids[gridI], gridI
+ grid_i = na.where((na.all(self.gridRightEdge > left_edge, axis=1)
+ & na.all(self.gridLeftEdge < right_edge, axis=1)) == True)
+ return self.grids[grid_i], grid_i
+
+ def get_periodic_box_grids(self, left_edge, right_edge):
+ mask = na.zeros(self.grids.shape, dtype='bool')
+ dl = self.parameters["DomainLeftEdge"]
+ dr = self.parameters["DomainRightEdge"]
+ db = right_edge - left_edge
+ for off_x in [-1, 0, 1]:
+ nle = left_edge[:]
+ nre = left_edge[:]
+ nle[0] = dl[0] + (dr[0]-dl[0])*off_x
+ for off_y in [-1, 0, 1]:
+ nle[1] = dl[1] + (dr[1]-dl[1])*off_y
+ for off_z in [-1, 0, 1]:
+ nle[2] = dl[2] + (dr[2]-dl[2])*off_z
+ nre = nle + db
+ g, gi = self.get_box_grids(nle, nre)
+ mask[gi] = True
+ return self.grids[mask], na.where(mask)
@time_execution
def find_max(self, field, finestLevels = True):
Modified: trunk/yt/lagos/PointCombine.c
==============================================================================
--- trunk/yt/lagos/PointCombine.c (original)
+++ trunk/yt/lagos/PointCombine.c Sat Jun 21 08:52:13 2008
@@ -531,20 +531,21 @@
int ll, i, n;
PyObject *og_le, *og_dx, *og_data, *og_cm,
- *oc_le, *oc_re, *oc_dx, *oc_data;
+ *oc_le, *oc_re, *oc_dx, *oc_data, *odr_edge, *odl_edge;
PyArrayObject *g_le, *g_dx, *g_cm,
- *c_le, *c_re, *c_dx;
+ *c_le, *c_re, *c_dx, *dr_edge, *dl_edge;
PyArrayObject **g_data, **c_data;
g_data = c_data = NULL;
npy_int *ag_cm;
npy_float64 ag_le[3], ag_dx[3],
- ac_le[3], ac_re[3], ac_dx[3];
+ ac_le[3], ac_re[3], ac_dx[3],
+ adl_edge[3], adr_edge[3];
Py_ssize_t n_fields = 0;
- if (!PyArg_ParseTuple(args, "OOOOOOOOi",
+ if (!PyArg_ParseTuple(args, "OOOOOOOOiOO",
&og_le, &og_dx, &og_data, &og_cm,
&oc_le, &oc_re, &oc_dx, &oc_data,
- &ll))
+ &ll, &odl_edge, &odr_edge))
return PyErr_Format(_dataCubeError,
"DataCubeGeneric: Invalid parameters.");
@@ -618,6 +619,26 @@
goto _fail;
}
+ dl_edge = (PyArrayObject *) PyArray_FromAny(odl_edge,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((dl_edge==NULL) || (PyArray_SIZE(dl_edge) != 3)) {
+ PyErr_Format(_dataCubeError,
+ "CombineGrids: Three values, one dimension required for dl_edge.");
+ goto _fail;
+ }
+ for(i=0;i<3;i++)adl_edge[i]=*(npy_float64*)PyArray_GETPTR1(dl_edge,i);
+
+ dr_edge = (PyArrayObject *) PyArray_FromAny(odr_edge,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((dr_edge==NULL) || (PyArray_SIZE(dr_edge) != 3)) {
+ PyErr_Format(_dataCubeError,
+ "CombineGrids: Three values, one dimension required for dr_edge.");
+ goto _fail;
+ }
+ for(i=0;i<3;i++)adr_edge[i]=*(npy_float64*)PyArray_GETPTR1(dr_edge,i);
+
n_fields = PyList_Size(oc_data);
if(n_fields == 0) {
PyErr_Format(_dataCubeError,
@@ -667,27 +688,63 @@
/* And let's begin */
npy_int64 xg, yg, zg, xc, yc, zc, cmax_x, cmax_y, cmax_z,
- cmin_x, cmin_y, cmin_z, cm;
+ cmin_x, cmin_y, cmin_z, cm, pxl, pyl, pzl;
npy_float64 *val1, *val2;
long int total=0;
+ npy_int64 p_niter[3] = {1,1,1};
+ npy_float64 ac_le_p[3][2];
+ npy_float64 ac_re_p[3][2];
+ 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);}
+
+ for (i=0;i<3;i++){ac_le_p[i][0] = ac_le[i]; ac_re_p[i][0] = ac_re[i];}
+ /*fprintf(stderr, "%0.3f %0.3f %0.3f %0.3f %0.3f %0.3f\n",
+ ac_le[0], ac_le[1], ac_le[2],
+ adl_edge[0], adl_edge[1], adl_edge[2]);
+ fprintf(stderr, "%0.3f %0.3f %0.3f %0.3f %0.3f %0.3f\n",
+ ac_re[0], ac_re[1], ac_re[2],
+ adr_edge[0], adr_edge[1], adr_edge[2]);*/
+ for(i=0;i<3;i++) {
+ n = 1;
+ if (ac_le[i] < adl_edge[i]){// && ag_le[i] > ac_re[i]) {
+ ac_le_p[i][n ] = adr_edge[i] + ac_le[i];
+ ac_re_p[i][n++] = adr_edge[i] + ac_re[i];
+ //fprintf(stderr, "le axis: %d le: %0.5f re: %0.5f\n",
+ // i, ac_le_p[i][n-1], ac_re_p[i][n-1]);
+ }
+ if (ac_re[i] > adr_edge[i]){// && ag_re[i] < ac_le[i]) {
+ ac_le_p[i][n ] = adl_edge[i] - (adr_edge[i]-adl_edge[i])
+ + ac_le[i];
+ ac_re_p[i][n++] = adl_edge[i] - (adr_edge[i]-adl_edge[i])
+ + ac_re[i];
+ //fprintf(stderr, "re axis: %d le: %0.5f re: %0.5f\n",
+ // i, ac_le_p[i][n-1], ac_re_p[i][n-1]);
+ }
+ p_niter[i] = n;
+ //fprintf(stderr, "Iterating on %d %d times\n", i, n);
+ }
+
for (xg = 0; xg < g_data[0]->dimensions[0]; xg++) {
- if (ag_le[0]+ag_dx[0]*xg > ac_re[0]) continue;
- if (ag_le[0]+ag_dx[0]*(xg+1) < ac_le[0]) continue;
- cmin_x = max(floorl((ag_le[0]+ag_dx[0]*xg - ac_le[0])/ac_dx[0]),0);
- cmax_x = min( ceill((ag_le[0]+ag_dx[0]*(xg+1) - ac_le[0])/ac_dx[0]),PyArray_DIM(c_data[0],0));
+ for (pxl = 0; pxl < p_niter[0]; pxl++) {
+ 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;
+ 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++) {
- if (ag_le[1]+ag_dx[1]*yg > ac_re[1]) continue;
- if (ag_le[1]+ag_dx[1]*(yg+1) < ac_le[1]) continue;
- cmin_y = max(floorl((ag_le[1]+ag_dx[1]*yg - ac_le[1])/ac_dx[1]),0);
- cmax_y = min( ceill((ag_le[1]+ag_dx[1]*(yg+1) - ac_le[1])/ac_dx[1]),PyArray_DIM(c_data[0],1));
+ for (pyl = 0; pyl < p_niter[1]; pyl++) {
+ 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++) {
+ for (pzl = 0; pzl < p_niter[2]; pzl++) {
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[2]) continue;
- if (ag_le[2]+ag_dx[2]*(zg+1) < ac_le[2]) continue;
- cmin_z = max(floorl((ag_le[2]+ag_dx[2]*zg - ac_le[2])/ac_dx[2]),0);
- cmax_z = min( ceill((ag_le[2]+ag_dx[2]*(zg+1) - ac_le[2])/ac_dx[2]),PyArray_DIM(c_data[0],2));
+ 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));
for (xc = cmin_x; xc < cmax_x ; xc++) {
for (yc = cmin_y; yc < cmax_y ; yc++) {
for (zc = cmin_z; zc < cmax_z ; zc++) {
@@ -699,9 +756,9 @@
}
}
}
- }
- }
- }
+ }}
+ }}
+ }}
/* Cleanup time */
@@ -711,6 +768,8 @@
Py_DECREF(c_le);
Py_DECREF(c_re);
Py_DECREF(c_dx);
+ Py_DECREF(dl_edge);
+ Py_DECREF(dr_edge);
for(n=0;n<n_fields;n++) {
Py_DECREF(g_data[n]);
Py_DECREF(c_data[n]);
More information about the yt-svn
mailing list