[Yt-svn] yt: Added ART IO methods. One is in python, in frontends/art/io....
hg at spacepope.org
hg at spacepope.org
Thu Oct 21 14:53:05 PDT 2010
hg Repository: yt
details: yt/rev/8f25a0593545
changeset: 3457:8f25a0593545
user: Christopher Erick Moody <cemoody at ucsc.edu>
date:
Thu Oct 21 14:52:55 2010 -0700
description:
Added ART IO methods. One is in python, in frontends/art/io.py and the other is functionally identical but in cython and in utilities/_amr_utils/fortran_reader.pyx
diffstat:
yt/frontends/art/data_structures.py | 21 +-
yt/frontends/art/io.py | 66 +-
yt/utilities/_amr_utils/fortran_reader.pyx | 106 +-
yt/utilities/amr_utils.c | 5404 +++++++++++++++++++--------
4 files changed, 3856 insertions(+), 1741 deletions(-)
diffs (truncated from 15993 to 300 lines):
diff -r 76ed5c3f466c -r 8f25a0593545 yt/frontends/art/data_structures.py
--- a/yt/frontends/art/data_structures.py Thu Oct 21 12:55:58 2010 -0700
+++ b/yt/frontends/art/data_structures.py Thu Oct 21 14:52:55 2010 -0700
@@ -128,10 +128,14 @@
def _count_grids(self):
# We have to do all the patch-coalescing here.
- level_info = [self.pf.ncell] # skip root grid for now
+ #level_info is used by the IO so promoting it to the static
+ # output class
+ self.pf.level_info = [self.pf.ncell] # skip root grid for now
+ #leve_info = []
amr_utils.count_art_octs(
self.pf.parameter_filename, self.pf.child_grid_offset,
- self.pf.min_level, self.pf.max_level, level_info)
+ self.pf.min_level, self.pf.max_level, self.pf.nhydro_vars,
+ self.pf.level_info)
num_ogrids = sum(level_info) + self.pf.iOctFree
ogrid_left_indices = na.zeros((num_ogrids,3), dtype='int64') - 999
ogrid_levels = na.zeros(num_ogrids, dtype='int64')
@@ -144,7 +148,7 @@
ogrid_left_indices, ogrid_levels,
ogrid_parents, ochild_masks)
ochild_masks.reshape((num_ogrids, 8), order="F")
- ogrid_levels[ogrid_left_indices[:,0] == -999] = -1
+ ogrid_levels[ogrid_left_indices[:,0] == -999] = 0
# This bit of code comes from Chris, and I'm still not sure I have a
# handle on what it does.
final_indices = ogrid_left_indices[na.where(ogrid_levels==self.pf.max_level)[0]]
@@ -344,14 +348,15 @@
storage_filename = None):
StaticOutput.__init__(self, filename, data_style)
self.storage_filename = storage_filename
-
+
self.field_info = self._fieldinfo_class()
self.current_time = 0.0
self.dimensionality = 3
self.refine_by = 2
self.parameters["HydroMethod"] = 'ramses'
self.parameters["Time"] = 1. # default unit is 1...
-
+ self.fh = open(storage_filename,'rb') #used by the io
+
def __repr__(self):
return self.basename.rsplit(".", 1)[0]
@@ -448,11 +453,13 @@
self.hubble_constant = header_vals['hubble']
self.min_level = header_vals['min_level']
self.max_level = header_vals['max_level']
-
+ self.nhydro_vars = 10 #this gets updated later, but we'll default to this
+ #nchem is nhydrovars-8, so we typically have 2 extra chem species
+
for to_skip in ['tl','dtl','tlold','dtlold','iSO']:
_skip_record(f)
- self.ncell = struct.unpack('>l', _read_record(f))
+ (self.ncell,) = struct.unpack('>l', _read_record(f))
# Try to figure out the root grid dimensions
est = na.log2(self.ncell) / 3
if int(est) != est: raise RuntimeError
diff -r 76ed5c3f466c -r 8f25a0593545 yt/frontends/art/io.py
--- a/yt/frontends/art/io.py Thu Oct 21 12:55:58 2010 -0700
+++ b/yt/frontends/art/io.py Thu Oct 21 14:52:55 2010 -0700
@@ -30,12 +30,68 @@
class IOHandlerART(BaseIOHandler):
_data_style = "art"
-
+ #at which position in the child record does the field occur
+ field_dict = {'gas density':0}
+
def __init__(self, art_tree, *args, **kwargs):
- self.art_tree = ramses_tree
BaseIOHandler.__init__(self, *args, **kwargs)
- def _read_data_set(self, grid, field):
- fullfieldname = 'grid_fluid_'+field
- return self.hierarchy.pf.art[fullfieldname][grid.id]
+ def _read_data_set(self,grid,field):
+ fn = grid.pf.storage_filename
+ # because of reuse promote this handler to pf class?
+ min_level = grid.pf.min_level
+ max_level = grid.pf.max_level
+ nhydro_vars = grid.pf.nhydro_vars
+ grid_level = grid.Level
+ #grid_idc = (grid.id-1)*8+grid.pf.ncell+1 (verbatim analysis_ART)
+ grid_idc = (grid.id)*8+grid.pf.ncell #this gives only the first cell of an oct
+ header_offset = grid.pf.child_grid_offset
+
+ record_size = 2+1+nhydro_vars+3 #two pads, idc, hydro vars, 3 vars
+ #go past the previous levels
+ # skip the header
+ offset = (header_offset
+ #per level:
+ # 4 bytes per integer, float or pad
+ # first section has three integers + 2 pads
+ + 4*5*(grid_level)
+ # second section has 13 floats +2 pads per oct
+ # (level_info) is the number of octs per level
+ + 4*15*sum(level_info[:grid_level])
+ # after the oct section is the child section.
+ # there are 2 pads, 1 integer child ID (idc)
+ # then #nhydro_vars of floats + 3 vars
+ #there 8 times as many children as octs
+ + 4*(2+1+nhydro_vars+3)*8*sum(level_info[:grid_level-1]))
+ pdb.set_trace()
+ fh = open(fn,'rb')
+ fh.seek(offset)
+ dtype='>i,'+'>i,'+'>%df,'%(nhydro_vars+3)+'>i'
+ first_idc = na.fromfile(fh,dtype=dtype,count=1)[0][1]
+ #this is the first idc of our section. we can guess
+ # how many records away our idc is since idc increments
+ # by 1 for almost every record. because idc can increase
+ # by more than 1 sometimes, we will always overestimate
+ # the number of records we should skip ahead. so work
+ # backwards from our guess
+ seek_guess = (grid_idc - first_idc)*record_size*4
+ fh.seek(offset+seek_guess)
+ while True:
+ record = na.fromfile(fh,dtype=dtype,count=1)[0]
+ # make sure that the pad bytes line up
+ # otherwise we may have a phase offset or have stepped out
+ # our section
+ assert record[0]==record[-1]
+ # we better have overestimated the idc
+ assert record[1]>= grid_idc
+ if record[1] == grid_idc:
+ fh.close()
+ return record[2][field_dict[field]]
+ else:
+ # in the next iteration we'll read the previous record
+ # so rewind the last read, then rewind one further record
+ fh.seek(fh.tell()-2*record_size)
+ fh.close()
+
+
\ No newline at end of file
diff -r 76ed5c3f466c -r 8f25a0593545 yt/utilities/_amr_utils/fortran_reader.pyx
--- a/yt/utilities/_amr_utils/fortran_reader.pyx Thu Oct 21 12:55:58 2010 -0700
+++ b/yt/utilities/_amr_utils/fortran_reader.pyx Thu Oct 21 14:52:55 2010 -0700
@@ -83,6 +83,7 @@
def count_art_octs(char *fn, long offset,
int min_level, int max_level,
+ int nhydro_vars,
level_info):
cdef int nchild = 8
cdef int i, Lev, next_record, nLevel
@@ -93,6 +94,7 @@
for Lev in range(min_level + 1, max_level + 1):
fread(dummy_records, sizeof(int), 2, f);
fread(&nLevel, sizeof(int), 1, f); FIX_LONG(nLevel)
+ print level_info
level_info.append(nLevel)
fread(dummy_records, sizeof(int), 2, f);
fread(&next_record, sizeof(int), 1, f); FIX_LONG(next_record)
@@ -102,17 +104,23 @@
fseek(f, next_record, SEEK_CUR)
# Now we skip the second section
fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
+ nhydro_vars = next_record/4-2-3 #nhvar in daniel's code
+ #record length is normally 2 pad bytes, 8 + 2 hvars (the 2 is nchem)
+ # and then 3 vars, but we can find nhvars only here and not in other
+ # file headers
next_record = (2*sizeof(int) + readin) * (nLevel * nchild)
next_record -= sizeof(int)
fseek(f, next_record, SEEK_CUR)
+ print "nhvars",nhydro_vars
fclose(f)
def read_art_tree(char *fn, long offset,
- int min_level, int max_level,
+ int min_level, int max_level, int nhydro_vars,
np.ndarray[np.int64_t, ndim=2] oct_indices,
np.ndarray[np.int64_t, ndim=1] oct_levels,
np.ndarray[np.int64_t, ndim=1] oct_parents,
- np.ndarray[np.int64_t, ndim=1] oct_mask):
+ np.ndarray[np.int64_t, ndim=1] oct_mask
+ ):
# This accepts the filename of the ART header and an integer offset that
# points to the start of the record *following* the reading of iOctFree and
# nOct. For those following along at home, we only need to read:
@@ -143,7 +151,7 @@
print "Reading Hierarchy for Level", Lev, Level, nLevel, iOct
total_cells += nLevel
for ic1 in range(nLevel):
- iOctMax = imax(iOctMax, iOct)
+ iOctMax = max(iOctMax, iOct)
#print readin, iOct, nLevel, sizeof(int)
next_record = ftell(f)
fread(&readin, sizeof(int), 1, f); FIX_LONG(readin)
@@ -172,3 +180,95 @@
print "Masked cells", total_masked
fclose(f)
print "Read this many cells", total_cells, iOctMax
+
+def read_art_vars(char *fn,
+ int min_level, int max_level, int nhydro_vars,
+ int grid_level,long grid_idc,long header_offset,
+ np.ndarray[np.int64_t, ndim=1] level_info,
+ np.ndarray[np.float64_t, ndim=1] hvars,
+ np.ndarray[np.float64_t, ndim=1] var):
+ # nhydro_vars is the number of columns- 3 (adjusting for vars)
+ # this is normally 10=(8+2chem species)
+ cdef FILE *f = fopen(fn, "rb")
+ cdef long offset
+ cdef long nocts,nc #number of octs/children in previous levels
+ cdef int j,lev
+
+ #parameters
+ cdef int record_size = 2+1+nhydro_vars+3
+
+ #record values
+ cdef int pada,padb
+ cdef int idc,idc_readin
+ cdef float temp
+ offset = 0
+
+ #total number of octs in previous levels
+ # including current level
+ nocts=0
+ for lev in range(min_level-1,grid_level):
+ #print lev
+ nocts += level_info[lev]
+ #total number of children in prev levels
+ # not including the current level
+ # there are 8 children for every oct
+ nc=0
+ for lev in range(min_level-1,grid_level-1):
+ #print lev
+ nc += 8*level_info[lev]
+
+
+ #skip the header
+ offset += header_offset
+ #per level:
+ # 4 bytes per integer, float or pad
+ # first section has three integers + 2 pads
+ offset += 4*5*grid_level
+ # second section has 13 floats +2 pads per oct
+ offset += 4*15*nocts
+ # after the oct section is the child section.
+ # there are 2 pads, 1 integer child ID (idc)
+ # then #nhydro_vars of floats + 3 vars
+ offset += + 4*(record_size)*nc
+
+ #now we read in the first idc, then make our
+ #seek guess
+ #print 'entry',offset
+ fseek(f, offset, SEEK_SET)
+ fread(&pada, sizeof(int), 1, f); FIX_LONG(pada)
+ fread(&idc_readin, sizeof(int), 1, f); FIX_LONG(idc_readin)
+
+ offset += (grid_idc - idc_readin)*record_size*4
+ fseek(f, offset, SEEK_SET)
+ j=0
+
+ while j<100:
+ fread(&pada, sizeof(int), 1, f); FIX_LONG(pada)
+ fread(&idc_readin, sizeof(int), 1, f); FIX_LONG(idc_readin)
+ if grid_idc != idc_readin:
+ # in the next iteration we'll read the previous record
+ # so rewind one further record
+ j+=1
+ offset = offset - record_size*4
+ fseek(f, offset, SEEK_SET)
+ else: break
+ for j in range(nhydro_vars):
+ fread(&temp, sizeof(float), 1, f);
+ FIX_FLOAT(temp)
+ hvars[j] = temp
+ for j in range(3):
+ fread(&temp, sizeof(float), 1, f);
+ FIX_FLOAT(temp)
+ var[j] = temp
+ fclose(f)
+
+
+
+
+
+
+
+
+
+
+
diff -r 76ed5c3f466c -r 8f25a0593545 yt/utilities/amr_utils.c
--- a/yt/utilities/amr_utils.c Thu Oct 21 12:55:58 2010 -0700
+++ b/yt/utilities/amr_utils.c Thu Oct 21 14:52:55 2010 -0700
@@ -1,4 +1,4 @@
-/* Generated by Cython 0.13.beta0 on Fri Sep 17 12:25:22 2010 */
+/* Generated by Cython 0.13 on Thu Oct 21 14:40:23 2010 */
#define PY_SSIZE_T_CLEAN
#include "Python.h"
More information about the yt-svn
mailing list