1. Matthew Turk
  2. yt-3.0

Commits

Sam Leitner  committed 27be21a

rewrite and clean up of grid_pos. first crack at working fill.

  • Participants
  • Parent commits ff52d81
  • Branches yt-3.0

Comments (0)

Files changed (6)

File yt/frontends/artio/_artio_caller.pyx

View file
 """
 
 """
-import sys #snl this is for debugging
+import numpy as np
+cimport numpy as np
+import sys 
 
 from libc.stdint cimport int32_t, int64_t
 from libc.stdlib cimport malloc, free
-#from  data_structures  import cell_pos_callback
 import  data_structures  
 
 cdef struct artio_context_struct: 
     ctypedef struct artio_file_struct "artio_file_struct" :
         pass
     ctypedef artio_file_struct *artio_file "artio_file"
+    ctypedef struct artio_grid_file_struct "artio_grid_file_struct" :
+        pass
+    ctypedef artio_grid_file_struct *artio_grid_file "artio_grid_file"
+
 
     ctypedef struct artio_context_struct "artio_context_struct" :
         pass
                                   float *variables, int32_t *num_oct_levels, int32_t *num_octs_per_level)
     int artio_grid_cache_sfc_range(artio_file handle, int64_t start, int64_t end)
 
-    ctypedef void (* GridCallBackPos)(float * variables, int level, int refined, int64_t sfc_index, double pos[3])
+    ctypedef void (* GridCallBackPos)(float * variables, int level, int refined, int64_t sfc_index, double pos[3], void *)
     int artio_grid_read_sfc_range_pos(artio_file handle,\
                 int64_t sfc1, int64_t sfc2,\
                 int min_level_to_read, int max_level_to_read,\
                 int options,\
-                GridCallBackPos callback)
+                GridCallBackPos callback, void *user_data)
+
+    ctypedef void (* GridCallBackBuffer)(float * variables, int level, int refined, int64_t sfc_index, void *)
+    int artio_grid_read_sfc_range_buffer(artio_file handle,\
+                int64_t sfc1, int64_t sfc2,\
+                int min_level_to_read, int max_level_to_read,\
+                int options,\
+                GridCallBackBuffer callback, void *user_data)
+
     ctypedef void (* GridCallBack)(float * variables, int level, int refined,int64_t sfc_index)
     int artio_grid_read_sfc_range(artio_file handle,\
                 int64_t sfc1, int64_t sfc2,\
 
             self.parameters[key] = parameter
 
-#        print self.parameters
-        
-#    def read_parameters(self) : 
-#        cdef char key[64]
-#        cdef int type
-#        cdef int length
-#        while artio_parameter_iterate( self.handle, key, &type, &length ) == ARTIO_SUCCESS :
-#            print 'hi!!'
-def check_artio_status(status, name):
-        if status!=ARTIO_SUCCESS :
-            print 'failure with status', status, 'in function ',name 
-            sys.exit(1)
+def check_artio_status(status, fname):
+    callername = sys._getframe().f_code.co_name
+    nline = sys._getframe().f_lineno
+    if status!=ARTIO_SUCCESS :
+        print 'failure with status', status, 'in function',fname,'from caller', callername, nline 
+        sys.exit(1)
 cdef class artio_fileset :
     cdef public object parameters 
     cdef artio_file handle
         print 'done reading header parameters'
 
 cdef class artio_fileset_grid :
-    cdef public object parameters #is self.parameters public or is parameters public?
+    cdef public object parameters
     cdef artio_file handle
+    cdef public file_prefix
 
     def __init__(self, char *file_prefix) :
         cdef int artio_type = ARTIO_OPEN_GRID
         d = read_parameters_artio(file_prefix, artio_type)
         self.parameters = {}
         self.parameters = d.parameters
-#        print self.parameters
+        self.file_prefix = file_prefix
         print 'done reading grid parameters'
 
-###### callback for positions #############
-def count_octs(char *file_prefix,
-               int64_t sfc1, int64_t sfc2,
-               int min_level_to_read, int max_level_to_read,
-               int num_grid_variables
-               ) :
-    cdef float * variables  
-    cdef int32_t * num_oct_levels 
-    cdef int32_t * num_octs_per_level 
+#snl: subclass some of this?
+class artio_grid_routines(object) : 
+    def __init__(self, param_handle) :
+        self.oct_handler=None
+        self.source = None
+        self.level_count = None
 
-    length = num_grid_variables * 8
-    variables = <float *>malloc(length*sizeof(float))
-    num_oct_levels = <int32_t *>malloc(1*sizeof(int32_t))
-    length = max_level_to_read
-    num_octs_per_level = <int32_t *>malloc(length*sizeof(int32_t))
+        self.min_level_to_read = 0
+        self.max_level_to_read = param_handle.parameters['grid_max_level'][0]
+        self.sfc1 = 0
+        self.sfc2 = param_handle.parameters['grid_file_sfc_index'][1]-1
+        self.num_grid_variables = param_handle.parameters['num_grid_variables'][0]
+        self.param_handle=param_handle
+        self.ng=1 
+        self.cpu=0
+        self.domain_id=0
 
-    cdef artio_file handle
-    handle = artio_fileset_open( file_prefix, ARTIO_OPEN_GRID, artio_context_global ) 
-    artio_fileset_open_grid( handle ) 
+        self.grid_variable_labels=param_handle.parameters['grid_variable_labels']
+ 
+        # dictionary from artio to yt ... should go elsewhere
+        self.label_artio_yt = {}
+        for label in self.grid_variable_labels :
+            if label == 'HVAR_GAS_DENSITY' :
+                self.label_artio_yt[label] = 'Density'
+            else :
+                self.label_artio_yt[label] = label
+        print self.label_artio_yt
 
-    num_total_octs =0
-    status = artio_grid_cache_sfc_range(handle, sfc1, sfc2)
-    check_artio_status(status, count_octs.__name__)
-    for sfc in xrange(sfc1,sfc2):
-        status = artio_grid_read_root_nocts(handle, sfc,
-                                            variables, num_oct_levels,
-                                            num_octs_per_level)
-        check_artio_status(status, count_octs.__name__)
-        noct_levels = num_oct_levels[0]
-        count_level_octs = {}          
-        count_level_octs = [ num_octs_per_level[i] for i in xrange(noct_levels) ]
-        num_sfc_octs = sum(count_level_octs)
-        num_total_octs += num_sfc_octs
 
-    status = artio_grid_read_root_cell_end(handle)
-    check_artio_status(status, count_octs.__name__)
-    artio_fileset_close_grid(handle)
+        self.grid_variable_labels        
+        self.label_index = {}
+        self.matched_fields = []
+        # not sure about file handling
+        # self.handle = <object> handle #<------seg faults
+        # and you never bother to close all of these handles
+        
+    def count_octs(self) :
+        cdef int min_level_to_read = self.min_level_to_read
+        cdef int max_level_to_read = self.max_level_to_read
+        cdef int64_t sfc1 = self.sfc1
+        cdef int64_t sfc2 = self.sfc2
+        cdef num_grid_variables = self.num_grid_variables
+        
+        cdef float * variables  
+        cdef int32_t * num_oct_levels 
+        cdef int32_t * num_octs_per_level 
+        
+        length = num_grid_variables * 8
+        variables = <float *>malloc(length*sizeof(float))
+        num_oct_levels = <int32_t *>malloc(1*sizeof(int32_t))
+        length = max_level_to_read
+        num_octs_per_level = <int32_t *>malloc(length*sizeof(int32_t))
+        
+        cdef artio_file handle
+        handle = artio_fileset_open( self.param_handle.file_prefix, 
+                                     ARTIO_OPEN_GRID, artio_context_global ) 
+        artio_fileset_open_grid( handle ) 
+        
+        cdef int64_t num_total_octs =0
+        n_levels = max_level_to_read - min_level_to_read + 1
+        level_count = np.zeros(n_levels, dtype='int64')
 
-    return num_total_octs
+        status = artio_grid_cache_sfc_range(handle, sfc1, sfc2)
+        check_artio_status(status, artio_grid_routines.__name__)
+        for sfc in xrange(sfc1,sfc2):
+            status = artio_grid_read_root_nocts(handle, sfc,
+                                                variables, num_oct_levels,
+                                                num_octs_per_level)
+            check_artio_status(status, artio_grid_routines.__name__)
+            noct_levels = num_oct_levels[0]
+            count_level_octs = {}          
+            count_level_octs = [ num_octs_per_level[i] for i in xrange(noct_levels) ]
+            for level in xrange(noct_levels) : 
+                level_count[level] += count_level_octs[level] 
+            num_sfc_octs = sum(count_level_octs)
+            num_total_octs += num_sfc_octs
+            status = artio_grid_read_root_cell_end(handle)
+            check_artio_status(status, artio_grid_routines.__name__)
+            
+            #    check_artio_status(-1, count_octs.__name__)
+            # dont close file until the end of the object... add __del__: artio_fileset_close_grid(handle)
+        self.level_count = level_count
+        return num_total_octs
 
-###### callback for positions #############
-def grid_pos_fill(char *file_prefix,\
-                int64_t sfc1, int64_t sfc2,\
-                int min_level_to_read, int max_level_to_read) :
-    cdef artio_file handle
-    handle = artio_fileset_open( file_prefix, ARTIO_OPEN_GRID, artio_context_global ) 
-    artio_fileset_open_grid( handle ) 
-    status = artio_grid_read_sfc_range_pos(handle,\
-                sfc1, sfc2,\
-                min_level_to_read, max_level_to_read,\
-                ARTIO_READ_LEAFS,\
-                wrap_cell_pos_callback)
-    check_artio_status(status, grid_pos_fill.__name__)
+    def grid_pos_fill(self, oct_handler) :
+        self.oct_handler = oct_handler
+        cdef artio_file handle
+        if self.oct_handler == None :
+            print 'oct_handler is not assigned!'
+            sys.exit(1)
+        handle = artio_fileset_open( self.param_handle.file_prefix, 
+                                     ARTIO_OPEN_GRID, artio_context_global ) 
+        status = artio_grid_read_sfc_range_pos(handle,\
+                    self.sfc1, self.sfc2,\
+                    self.min_level_to_read, self.max_level_to_read,\
+                    ARTIO_READ_REFINED,\
+                    wrap_cell_pos_callback, <void*>self) 
+        check_artio_status(status, artio_grid_routines.__name__)
+        print 'done filling oct positions'
+    def cell_pos_callback(self, level, refined, sfc_index, pos):
+        self.oct_handler.add(self.cpu + 1, level - self.min_level_to_read, self.ng, pos, self.domain_id)
 
-cdef void wrap_cell_pos_callback(float *variables, int level, int refined, int64_t sfc_index, double *pos):
-    position = {}
-    position = [ pos[i] for i in range(3) ]
-    data_structures.cell_pos_callback(level, refined, sfc_index, position)
+    def grid_var_fill(self, source, fields):
+        self.source = source
+        i=-1
+        for artlabel in self.grid_variable_labels :
+            label = self.label_artio_yt[artlabel]
+            i=i+1
+            for field in fields : 
+                if field == label :
+                    print 'match, in fields?', field,label, fields
+                    print '!!!!!!!!!!!!!!!'
+                    self.label_index[field]=i
+                    self.matched_fields.append(field)
+        print 'matched fields:',self.matched_fields
+        print 'art index of matched fields',self.label_index
+        self.count=0
+        cdef artio_file handle
+        if len(self.label_index) > 0 :
+            handle = artio_fileset_open( self.param_handle.file_prefix, 
+                                         ARTIO_OPEN_GRID, artio_context_global ) 
+            status = artio_grid_read_sfc_range_buffer(
+                handle, self.sfc1, self.sfc2,\
+                    self.min_level_to_read, self.max_level_to_read,\
+                    ARTIO_READ_REFINED,\
+                    wrap_cell_var_callback,\
+                    <void*>self
+                ) #only octs!
+            check_artio_status(status, artio_grid_routines.__name__)
+        print 'done buffering variables'
+    def cell_var_callback(self, level, refined, ichild, cell_var):
+        for field in self.matched_fields : 
+            self.source[field][self.count] = cell_var[self.label_index[field]] 
+        self.count=self.count+1
+ 
 
-###### callback for variables #############
-def grid_var_fill(char *file_prefix,\
-                int64_t sfc1, int64_t sfc2,\
-                int min_level_to_read, int max_level_to_read) :
-    cdef artio_file handle
-    handle = artio_fileset_open( file_prefix, ARTIO_OPEN_GRID, artio_context_global ) 
-    artio_fileset_open_grid( handle ) 
-    artio_grid_read_sfc_range(handle,\
-                sfc1, sfc2,\
-                min_level_to_read, max_level_to_read,\
-                ARTIO_READ_LEAFS,\
-                wrap_cell_var_callback)
-def cell_var_callback(level, refined, sfc_index, cell_var):
-    print "variable callback success! ",level, refined, sfc_index 
-cdef void wrap_cell_var_callback(float *variables, int level, int refined, int64_t sfc_index):
+###### callbacks #############
+cdef void wrap_cell_pos_callback(float *variables, int level, int refined, 
+                                 int64_t sfc_index, double *pos, void *user_data):
+    position = np.empty((1, 3), dtype='float64')
+    position[0,0] = pos[0] 
+    position[0,1] = pos[1] 
+    position[0,2] = pos[2] 
+    artioroutines = <object>user_data
+    artioroutines.cell_pos_callback(level, refined, sfc_index, position)
+
+cdef void wrap_cell_var_callback(float *variables, int level, int refined, 
+                                 int64_t ichild, void *user_data):
+    artioroutines = <object>user_data
     cell_var={}
-#     cdef int num_grid_variables=1 # really from read_parameters_artio(file_prefix, artio_type)
-#        for label in parameters.grid_variable_labels
-#            cell_var = [ variables[i] for i in range(num_grid_variables) ]
-#            if(label == 'density'): 
-#                cell_var_callback(cell_var, level, refined, sfc_index)
-    cell_var_callback(level, refined, sfc_index, cell_var)
- 
+    cell_var = [ variables[i] for i in range(artioroutines.num_grid_variables) ]
+    artioroutines.cell_var_callback(level, refined, ichild, cell_var)
+
+
 def artio_is_valid( char *file_prefix ) :
     cdef artio_file handle = artio_fileset_open( file_prefix, 
             ARTIO_OPEN_HEADER, artio_context_global )
     else :
         artio_fileset_close(handle) 
     return True
+
+
+

File yt/frontends/artio/artio_headers/artio.h

View file
 /* allocation strategy */
 #define ARTIO_ALLOC_EQUAL_SFC               0
 #define ARTIO_ALLOC_EQUAL_PROC              1
-#define ARTIO_ALLOC_MAX_FILE_SIZE  	        2
+#define ARTIO_ALLOC_MAX_FILE_SIZE  	    2
 
 #define ARTIO_TYPE_STRING                   0
 #define ARTIO_TYPE_CHAR                     1
 typedef void (* GridCallBack)(float * variables, int level, int refined,
 		int64_t sfc_index);
 typedef void (* GridCallBackPos)(double * variables, int level, int refined,
-                int64_t sfc_index, double pos[3]);
+                                 int64_t sfc_index, double pos[3], void * user_data);
+typedef void (* GridCallBackBuffer)(double * variables, int level, int refined,
+                                 int64_t sfc_index, void * user_data);
 
 /*
  * Description:	Add a grid component to a fileset open for writing
  *  option			1. refined nodes; 2 leaf nodes; 3 all nodes
  *  callback			callback function
  */
-int artio_grid_read_sfc_range(artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read,int max_level_to_read, int options, GridCallBack callback);
-int artio_grid_read_sfc_range_pos(artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read,int max_level_to_read, int options, GridCallBackPos callback);
+int artio_grid_read_sfc_range(       artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read, int max_level_to_read, int options, GridCallBack callback);
+int artio_grid_read_sfc_range_pos(artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read,int max_level_to_read, int options, GridCallBackPos callback, void *user_data);
+int artio_grid_read_sfc_range_buffer(artio_file handle, int64_t sfc1, int64_t sfc2, int min_level_to_read, int max_level_to_read, int options, GridCallBackBuffer callback, void *user_data);
 		
 
 typedef void (* ParticleCallBack)(int64_t pid, 

File yt/frontends/artio/artio_headers/artio_grid.c

View file
 
 	return ARTIO_SUCCESS;
 }
+int artio_grid_read_sfc_range_buffer(artio_file handle, 
+		int64_t sfc1, int64_t sfc2,
+		int min_level_to_read, int max_level_to_read, 
+		int options,
+                GridCallBackBuffer callback, 
+                void *userdata) {
+	int64_t sfc;
+	int oct, level, j;
+	int ret;
+	int *octs_per_level = NULL;
+	int refined;
+	int oct_refined[8];
+	int root_tree_levels;
+	float * variables = NULL;
+
+	artio_grid_file ghandle;
+
+	if ( handle == NULL ) {
+		return ARTIO_ERR_INVALID_HANDLE;
+	}
+
+	if (handle->open_mode != ARTIO_FILESET_READ ||
+			!(handle->open_type & ARTIO_OPEN_GRID) ) {
+		return ARTIO_ERR_INVALID_FILESET_MODE;
+	}
+
+	ghandle = handle->grid;
+
+	if ((min_level_to_read < 0) || (min_level_to_read > max_level_to_read)) {
+		return ARTIO_ERR_INVALID_LEVEL;
+	}
+
+	octs_per_level = (int *)malloc(ghandle->file_max_level * sizeof(int));
+	variables = (float *)malloc(8*ghandle->num_grid_variables * sizeof(float));
+
+	if ( octs_per_level == NULL || variables == NULL ) {
+		if ( octs_per_level != NULL ) free(octs_per_level);
+		if ( variables != NULL ) free(variables);
+		return ARTIO_ERR_MEMORY_ALLOCATION;
+	}
+
+	ret = artio_grid_cache_sfc_range(handle, sfc1, sfc2);
+	if ( ret != ARTIO_SUCCESS ) {
+		free(octs_per_level);
+		free(variables);
+		return ret;
+	}
+		
+	for (sfc = sfc1; sfc <= sfc2; sfc++) {
+		ret = artio_grid_read_root_cell_begin(handle, sfc, variables,
+				&root_tree_levels, octs_per_level);
+		if ( ret != ARTIO_SUCCESS ) {
+			free(octs_per_level);
+			free(variables);
+			return ret;
+		}
+
+		if (min_level_to_read == 0 && (options == ARTIO_READ_ALL || 
+				(options == ARTIO_READ_REFINED && root_tree_levels > 0) || 
+				(options == ARTIO_READ_LEAFS && root_tree_levels == 0)) ) {
+			refined = (root_tree_levels > 0) ? 1 : 0;
+                        j=-1;
+			callback(variables, 0, refined, j, userdata);
+		}
+
+		for (level = MAX(min_level_to_read,1); 
+				level <= MIN(root_tree_levels,max_level_to_read); level++) {
+			ret = artio_grid_read_level_begin(handle, level);
+			if ( ret != ARTIO_SUCCESS ) {
+				free(octs_per_level);
+				free(variables);
+				return ret;
+			}
+
+			for (oct = 0; oct < octs_per_level[level - 1]; oct++) {
+				ret = artio_grid_read_oct(handle, variables, oct_refined);
+				if ( ret != ARTIO_SUCCESS ) {
+					free(octs_per_level);
+					free(variables);
+					return ret;
+				}
+
+				for (j = 0; j < 8; j++) {
+					if (options == ARTIO_READ_ALL || 
+							(options == ARTIO_READ_REFINED && oct_refined[j]) ||
+							(options == ARTIO_READ_LEAFS && !oct_refined[j]) ) {
+						callback(&variables[j * ghandle->num_grid_variables],
+                                                         level, oct_refined[j], j, userdata);
+					}
+				}
+			}
+			artio_grid_read_level_end(handle);
+		}
+		artio_grid_read_root_cell_end(handle);
+	}
+
+	free(variables);
+	free(octs_per_level);
+
+	return ARTIO_SUCCESS;
+}
 
 
 /* array which describes how child cells are offset from 
 #endif
 };
 int artio_grid_read_sfc_range_pos(artio_file handle, 
-		int64_t sfc1, int64_t sfc2,
-		int min_level_to_read, int max_level_to_read, 
-		int options,
-		GridCallBackPos callback) {
+                                  int64_t sfc1, int64_t sfc2,
+                                  int min_level_to_read, int max_level_to_read, 
+                                  int options,
+                                  GridCallBackPos callback, 
+                                  void *user_data
+    ) {
 	int64_t sfc;
 	int oct, level, j, i;
 	int ret;
                 //////////////////////////////////////
                         
 
-		if (min_level_to_read == 0 && (options == ARTIO_READ_ALL || 
-				(options == ARTIO_READ_REFINED && root_tree_levels > 0) || 
-				(options == ARTIO_READ_LEAFS && root_tree_levels == 0)) ) {
+		if (min_level_to_read == 0 && 
+                    (options == ARTIO_READ_ALL || 
+                     (options == ARTIO_READ_REFINED && root_tree_levels > 0) || 
+                     (options == ARTIO_READ_LEAFS && root_tree_levels == 0)) ) {
 			refined = (root_tree_levels > 0) ? 1 : 0;
-			callback(variables, min_level, refined, sfc, pos);
+			callback(variables, min_level, refined, sfc, pos, user_data);
 		}
-// level is the cell_level; current octs live at level-1.
+                // level is the cell_level; current octs live at level-1.
 		for (level = min_level+1; level <= MIN(root_tree_levels,max_level_to_read); level++) { 
 			ret = artio_grid_read_level_begin(handle, level);
 			if ( ret != ARTIO_SUCCESS ) {
                                                 num_next_level_octs++ ;
                                         }
 					if (options == ARTIO_READ_ALL || 
-							(options == ARTIO_READ_REFINED && oct_refined[j]) ||
-							(options == ARTIO_READ_LEAFS && !oct_refined[j]) ) {
-                                            callback(variables, level, oct_refined[j], sfc, pos);
+                                            (options == ARTIO_READ_REFINED && oct_refined[j]) ||
+                                            (options == ARTIO_READ_LEAFS && !oct_refined[j]) ) {
+                                            callback(variables, level, oct_refined[j], sfc, pos, user_data);
 								
 					}
 				}

File yt/frontends/artio/data_structures.py

View file
   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
-
 import numpy as np
 import stat
 import weakref
 
 from _artio_caller import \
     artio_is_valid, artio_fileset , artio_fileset_grid, \
-    grid_pos_fill, grid_var_fill, count_octs#, read_header 
+    artio_grid_routines#, read_header 
 from yt.utilities.definitions import \
     mpc_conversion, sec_conversion 
 from .fields import ARTIOFieldInfo, KnownARTIOFields
 
 
 
-def cell_pos_callback( level, refined, sfc_index, pos):
-    print pos[0],pos[1],pos[2], 'cell position from callback'
-    ng=1 #until you write something to buffer
-    cpu=0
-    domain_id=0
-    oct_handler.add(cpu + 1, level - min_level, ng, pos, domain_id) 
-
-# the following classes are old RAMSES oct handling until otherwise noted
 class ARTIODomainFile(object):
-    _handle = None
+#    _last_mask = None
+#    _last_selector_id = None
+#    nvar = 6
+#    _handle = None
 
     def __init__(self, pf, domain_id):
         self.pf = pf
         self.domain_id = domain_id
         self.part_fn = "%s.p%03i" % (pf.parameter_filename[:-4],domain_id)
         self.grid_fn = "%s.g%03i" % (pf.parameter_filename[:-4],domain_id)
-        #self.parameter_file = weakref.proxy(pf)
         self._fileset_prefix = pf.parameter_filename[:-4]
-        self._handle = artio_fileset(self._fileset_prefix) 
-
+        self._handle = artio_fileset_grid(self._fileset_prefix) 
+        
+        self.artiogrid = artio_grid_routines(self._handle)
         self._read_grid_header()
         self._read_particle_header()
+        self.level_count
 
-#     @property
-#     def level_count(self):
-#         if self._level_count is not None: return self._level_count
-#         self.hydro_offset
-#         return self._level_count
-
-#     @property
-#     def hydro_offset(self):
-#         if self._hydro_offset is not None: return self._hydro_offset
-#         # We now have to open the file and calculate it
-#         f = open(self.hydro_fn, "rb")
-#         fpu.skip(f, 6)
-#         # It goes: level, CPU, 8-variable
-#         min_level = self.pf.min_level
-#         n_levels = self.grid_header['nlevelmax'] - min_level
-#         hydro_offset = np.zeros(n_levels, dtype='int64')
-#         hydro_offset -= 1
-#         level_count = np.zeros(n_levels, dtype='int64')
-#         for level in range(self.grid_header['nlevelmax']):
-#             for cpu in range(self.grid_header['nboundary'] +
-#                              self.grid_header['ncpu']):
-#                 header = ( ('file_ilevel', 1, 'I'),
-#                            ('file_ncache', 1, 'I') )
-#                 hvals = fpu.read_attrs(f, header)
-#                 if hvals['file_ncache'] == 0: continue
-#                 assert(hvals['file_ilevel'] == level+1)
-#                 if cpu + 1 == self.domain_id and level >= min_level:
-#                     hydro_offset[level - min_level] = f.tell()
-#                     level_count[level - min_level] = hvals['file_ncache']
-#                 fpu.skip(f, 8 * self.nvar)
-#         self._hydro_offset = hydro_offset
-#         self._level_count = level_count
-#         return self._hydro_offset
-
+    ###########################################
     def _read_particle_header(self):
         if not os.path.exists(self.part_fn):
             self.local_particle_count = 0
             self.particle_field_offsets = {}
             return
-#        f = open(self.part_fn, "rb")
-        hvals = {}
-        attrs = ( ('ncpu', 1, 'I'),
-                  ('ndim', 1, 'I'),
-                  ('npart', 1, 'I') )
-#        hvals.update(fpu.read_attrs(f, attrs))
-#        fpu.read_vector(f, 'I')
-
-        attrs = ( ('nstar_tot', 1, 'I'),
-                  ('mstar_tot', 1, 'd'),
-                  ('mstar_lost', 1, 'd'),
-                  ('nsink', 1, 'I') )
-#        hvals.update(fpu.read_attrs(f, attrs))
-#        self.particle_header = hvals
-#        self.local_particle_count = hvals['npart']
-        particle_fields = [
-                ("particle_position_x", "d"),
-                ("particle_position_y", "d"),
-                ("particle_position_z", "d"),
-                ("particle_velocity_x", "d"),
-                ("particle_velocity_y", "d"),
-                ("particle_velocity_z", "d"),
-                ("particle_mass", "d"),
-                ("particle_identifier", "I"),
-                ("particle_refinement_level", "I")]
-#        if hvals["nstar_tot"] > 0:
-#            particle_fields += [("particle_age", "d"),
-#                                ("particle_metallicity", "d")]
-#        field_offsets = {particle_fields[0][0]: f.tell()}
-#        for field, vtype in particle_fields[1:]:
-#            fpu.skip(f, 1)
-#            field_offsets[field] = f.tell()
-#        self.particle_field_offsets = field_offsets
 
     def _read_grid_header(self):
-#        hvals = {}
-#        f = open(self.grid_fn, "rb")
-#        for header in ramses_header(hvals):
-#            hvals.update(fpu.read_attrs(f, header))
-        # That's the header, now we skip a few.
-#        hvals['numbl'] = np.array(hvals['numbl']).reshape(
-#            (hvals['nlevelmax'], hvals['ncpu']))
-#        fpu.skip(f)
-#        if hvals['nboundary'] > 0:
-#            fpu.skip(f, 2)
-#            self.ngridbound = fpu.read_vector(f, 'i').astype("int64")
-#        else:
-#            self.ngridbound = np.zeros(hvals['nlevelmax'], dtype='int64')
-#        free_mem = fpu.read_attrs(f, (('free_mem', 5, 'i'), ) )
-#        ordering = fpu.read_vector(f, 'c')
-#        fpu.skip(f, 4)
-        # Now we're at the tree itself
-        # Now we iterate over each level and each CPU.
-#        self.grid_header = hvals
-#        self.grid_offset = f.tell()
-#        self.local_oct_count = hvals['numbl'][self.pf.min_level:, self.domain_id - 1].sum()
+        sfc_max = self._handle.parameters['grid_file_sfc_index'][1]-1
+        self.local_oct_count = self.artiogrid.count_octs() #count_octs(self._handle)
+        self.level_count = self.artiogrid.level_count
+        print 'snl data_structures local_octs count',self.local_oct_count
+    ###########################################
         
-        sfc_max = self._handle.parameters['grid_file_sfc_index'][1]-1
-        self.local_oct_count = \
-            count_octs(self._fileset_prefix,
-                       0,sfc_max, 
-                       self.pf.min_level,
-                       self._handle.parameters['grid_max_level'][0], 
-                       self._handle.parameters['num_grid_variables'][0])
-        print self.local_oct_count
-        
-
-
     def _read_grid(self, oct_handler):
         """Open the oct file, read in octs level-by-level.
            For each oct, only the position, index, level and domain 
            The most important is finding all the information to feed
            oct_handler.add
         """
-        sfc_max = self._handle.parameters['grid_file_sfc_index'][1]-1
-        sfc_max = 5   # for testing
-        grid_pos_fill(self._fileset_prefix, 0,sfc_max, 0, 10)##########
-
-    #snl this is where I would LIKE the callback to go:
-#    def cell_pos_callback(self, level, refined, sfc_index, pos):
-#        print pos[0],pos[1],pos[2], 'snlpos'
-#        ng=1 #until you write something to buffer
-#        cpu=0
-#        self.domain_id=0
-#        oct_handler.add(cpu + 1, level - min_level, ng, pos, self.domain_id) 
+        #oct_handler may not be attributed yet
+        self.artiogrid=artio_grid_routines(self._handle) 
+        self.artiogrid.grid_pos_fill(oct_handler)
+        print 'segfault debug'
         
-    
-
     def select(self, selector):
         if id(selector) == self._last_selector_id:
             return self._last_mask
         return self.count(selector)
 
 class ARTIODomainSubset(object):
+
     def __init__(self, domain, mask, cell_count):
         self.mask = mask
         self.domain = domain
         level_counts[1:] = level_counts[:-1]
         level_counts[0] = 0
         self.level_counts = np.add.accumulate(level_counts)
-
+        print 'cumulative masked level counts',self.level_counts
+        
+        self.artiogrid = self.domain.artiogrid
+        #self.artiogrid=artio_grid_routines(self._handle) 
+        
     def icoords(self, dobj):
         return self.oct_handler.icoords(self.domain.domain_id, self.mask,
                                         self.cell_count,
                                      self.cell_count,
                                      self.level_counts.copy())
 
-    #unchanged...
-    def fill(self, content, fields):
+    def fill(self, fields):
+        oct_handler = self.oct_handler
+        all_fields = self.domain.pf.h.fluid_field_list
+        min_level = self.domain.pf.min_level
+        fields = [f for ft, f in fields]
+        print 'all_fields:', all_fields 
+        print 'fields:', fields
+        
+
+
+        tr = {}
+        for field in fields: 
+            tr[field] = np.zeros(self.cell_count, 'float64')
+        temp = {}
+        nc = self.domain.level_count.sum()
+        print 'ncells total',self.cell_count, 'ncells masked', nc
+        for field in fields:
+            #temp[field] = np.empty((no,8), dtype="float64") 
+            temp[field] = np.empty(nc, dtype="float32") 
+
+        #buffer variables 
+        self.artiogrid.grid_var_fill(temp, fields)
+        
+        #mask unused cells 
+        oct_handler = self.oct_handler
+        level_offset = 0
+        level_offset += oct_handler.fill_mask(
+             self.domain.domain_id,
+            tr, temp, self.mask, level_offset) #[oct_container.pyx]
+
+        return tr
+    
+    
+
+    def fill_old(self, content, fields):
         # Here we get a copy of the file, which we skip through and read the
         # bits we want.
         oct_handler = self.oct_handler
         tr = {}
         filled = pos = level_offset = 0
         min_level = self.domain.pf.min_level
+        print 'hydro offset', self.domain.hydro_offset
+        
         for field in fields:
             tr[field] = np.zeros(self.cell_count, 'float64')
         for level, offset in enumerate(self.domain.hydro_offset):
+            print 'level', level, 'offset', offset
             if offset == -1: continue
             content.seek(offset)
             nc = self.domain.level_count[level]
                         #print "Reading %s in %s : %s" % (field, level,
                         #        self.domain.domain_id)
                         temp[field][:,i] = fpu.read_vector(content, 'd') # cell 1
+                        #print field, temp[field][:,i] 
             level_offset += oct_handler.fill_level(self.domain.domain_id, level,
                                    tr, temp, self.mask, level_offset)
             #print "FILL (%s : %s) %s" % (self.domain.domain_id, level, level_offset)
         #print "DONE (%s) %s of %s" % (self.domain.domain_id, level_offset, self.cell_count)
+        raise NotImplementedError #snl
         return tr
 
 class ARTIOGeometryHandler(OctreeGeometryHandler):
         super(ARTIOGeometryHandler, self).__init__(pf, data_style)
 
     def _initialize_oct_handler(self):
-        #domains are the class object ... ncpu == 1 for parallelization currently 
+        #domains are the class object ... ncpu == 1 currently 
         #only one file/"domain" for ART
         self.domains = [ARTIODomainFile(self.parameter_file, i + 1)
                         for i in range(self.parameter_file['ncpu'])]
-        #this merely allocates space for the oct tree
-        #and nothing else
+        #this allocates space for the oct tree
         self.oct_handler = ARTIOOctreeContainer(
-            self.parameter_file.domain_dimensions/2,
+            self.parameter_file.domain_dimensions, 
             self.parameter_file.domain_left_edge,
-            self.parameter_file.domain_right_edge)
+            self.parameter_file.domain_right_edge) 
         mylog.debug("Allocating octs")
         self.oct_handler.allocate_domains(
             [dom.local_oct_count for dom in self.domains])
         for dom in self.domains:
             dom._read_grid(self.oct_handler)
+        print 'seg fault debugging in _initialize_oct_handler' 
 
     def _detect_fields(self):
-        # TODO: Add additional fields #snl how are field names dealt with in ART?
+        # snl Add additional fields and translator from artio <-> yt
         self.fluid_field_list = [ "Density", "x-velocity", "y-velocity",
 	                        "z-velocity", "Pressure", "Metallicity" ]
         pfl = set([])
         raise NotImplementedError
 
     def _chunk_io(self, dobj):
-        raise NotImplementedError
-#        oobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
-#        for subset in oobjs:
-#            yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
+        # _current_chunk is made from identify_base_chunk 
+        #object = dobj._current_chunk.objs or dobj._current_chunk.${dobj._chunk_info}
+        oobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+        for subset in oobjs:
+            yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
 
 
-
-#snl: the following is for complete and understood header reading and unit setting
 class ARTIOStaticOutput(StaticOutput):
     _handle = None
     _hierarchy_class = ARTIOGeometryHandler
         self._filename = filename
         self._fileset_prefix = filename[:-4]
         self._handle = artio_fileset(self._fileset_prefix) 
-#        grid_pos_fill(self._fileset_prefix, 0,sfc_max, 0, 10)##########
-#        grid_var_fill(self._fileset_prefix, 0,100, 0, 10)##########
-        
         # Here we want to initiate a traceback, if the reader is not built.
         StaticOutput.__init__(self, filename, data_style)
         self.storage_filename = storage_filename 
         self.refine_by = 2
         self.parameters["HydroMethod"] = 'artio'
         self.parameters["Time"] = 1. # default unit is 1...
-        self.min_level = 0  #self._handle.parameters['grid_min_level']
-        # ^hard-coded
+        self.min_level = 0  # ART has min_level=0. No self._handle.parameters['grid_min_level']
 
         # read header
         self.unique_identifier = \
             int(os.stat(self.parameter_filename)[stat.ST_CTIME])
 
-        self.domain_dimensions = np.ones(3,dtype='int64') * \
-            self._handle.parameters["num_root_cells"][0]**(1./3.)
-
+        num_grid = (int)(self._handle.parameters["num_root_cells"][0]**(1./3.)+0.5)
+        self.domain_dimensions = np.ones(3,dtype='int32') * num_grid
         self.domain_left_edge = np.zeros(3, dtype="float64")
-        self.domain_right_edge = self.domain_dimensions
+        self.domain_right_edge = np.ones(3, dtype='float64')*num_grid
 
         self.min_level = 0
         self.max_level = self._handle.parameters["max_refinement_level"][0]
         self.parameters['unit_t'] = self._handle.parameters["time_unit"][0]
         self.parameters['unit_m'] = self._handle.parameters["mass_unit"][0]
         self.parameters['unit_d'] = self.parameters['unit_m']/self.parameters['unit_l']**(3.0)
-#snl hard coded number of domains in ART = 1 ... that may change for parallelization 
+
+        # hard coded number of domains in ART = 1 ... that may change for parallelization 
         self.parameters['ncpu'] = 1
 
     @classmethod

File yt/frontends/artio/fields.py

View file
 
 #need all of these to be defined
 known_artio_fields = [
+    "VAR_POTENTIAL",
     "Density",
     "x-velocity",
     "y-velocity",
     f.take_log = False
 
 known_artio_particle_fields = [
-    "particle_position_x",
-    "particle_position_y",
-    "particle_position_z",
-    "particle_velocity_x",
-    "particle_velocity_y",
-    "particle_velocity_z",
+    "POSITION_X",
+    "POSITION_Y",
+    "POSITION_Z",
+    "VELOCITY_X",
+    "VELOCITY_Y",
+    "VELOCITY_Z",
+#    "particle_position_x",
+#    "particle_position_y",
+#    "particle_position_z",
+#    "particle_velocity_x",
+#    "particle_velocity_y",
+#    "particle_velocity_z",
     "particle_mass",
     "particle_identifier",
-    "particle_refinement_level",
+    "particle_refinement_level"
 ]
 
 for f in known_artio_particle_fields:

File yt/frontends/artio/io.py

View file
   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 """
-
 from collections import defaultdict
 import numpy as np
 
     _data_style = "artio"
 
     def _read_fluid_selection(self, chunks, selector, fields, size):
-#chunks: list of YTDataChunk objects
-#selector: SelectorObjects (you can see the definition in yt/geometry/selection_routines.pxd) which can tell you which cells or particles to include in an object
-#fields: list of fields of the form (ftype, fname) where ftype is a string (for the case of multi-fluid simulations) and fname is the name of the fluid type.
-#size: the total (expected) number of cells that intersect, so that fields can be pre-allocated (to avoid memory fragmentation)
-
-        # Chunks in this case will have affiliated domain subset objects
-        # Each domain subset will contain a hydro_offset array, which gives
-        # pointers to level-by-level hydro information
         tr = dict((f, np.empty(size, dtype='float64')) for f in fields)
         cp = 0
-        for chunk in chunks:     #from _chunk_io in grid_geometry_handler.py yielding YTDataChunk 
-            for subset in chunk.objs:        #chunk.objs =  gs = files  [=gobjs appended to gfiles]
-                # Now we read the entire thing
-                f = open(subset.domain.grid_fn, "rb") 
-                # This contains the boundary information, so we skim through
-                # and pick off the right vectors
-                content = cStringIO.StringIO(f.read())
-                rv = subset.fill(content, fields)
+        for chunk in chunks:
+            for subset in chunk.objs:
+                print 'reading from ',fields, subset.domain.grid_fn
+                rv = subset.fill(fields) 
                 for ft, f in fields:
                     mylog.debug("Filling %s with %s (%0.3e %0.3e) (%s:%s)",
                         f, subset.cell_count, rv[f].min(), rv[f].max(),
         return tr
 
     def _read_particle_selection(self, chunks, selector, fields):
+        raise NotImplementedError #snl 
+
         size = 0
         masks = {}
         for chunk in chunks: