[yt-svn] commit/yt: atmyers: Merged in xarthisius/yt (pull request #2117)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Fri Apr 29 12:39:20 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/2d7b93e59428/
Changeset: 2d7b93e59428
Branch: yt
User: atmyers
Date: 2016-04-29 19:39:05+00:00
Summary: Merged in xarthisius/yt (pull request #2117)
Optimize creation of masks and hdf5 reads during creation of ionization tables for OWLS. Closes #1206
Affected #: 2 files
diff -r 3f50767dec553868bf966a8e027f308f545e84c2 -r 2d7b93e59428f2dc6d4cde5030cacc4b113c6ea4 yt/frontends/owls/fields.py
--- a/yt/frontends/owls/fields.py
+++ b/yt/frontends/owls/fields.py
@@ -222,41 +222,40 @@
""" returns a function that calculates the ion density of a particle.
"""
- def _ion_density(field, data):
+ def get_owls_ion_density_field(ion, ftype, itab):
+ def _func(field, data):
- # get element symbol from ion string. ion string will
- # be a member of the tuple _ions (i.e. si13)
- #--------------------------------------------------------
- if ion[0:2].isalpha():
- symbol = ion[0:2].capitalize()
- else:
- symbol = ion[0:1].capitalize()
+ # get element symbol from ion string. ion string will
+ # be a member of the tuple _ions (i.e. si13)
+ #--------------------------------------------------------
+ if ion[0:2].isalpha():
+ symbol = ion[0:2].capitalize()
+ else:
+ symbol = ion[0:1].capitalize()
- # mass fraction for the element
- #--------------------------------------------------------
- m_frac = data[ftype, symbol+"_fraction"]
+ # mass fraction for the element
+ #--------------------------------------------------------
+ m_frac = data[ftype, symbol+"_fraction"]
- # get nH and T for lookup
- #--------------------------------------------------------
- log_nH = np.log10( data["PartType0", "H_number_density"] )
- log_T = np.log10( data["PartType0", "Temperature"] )
+ # get nH and T for lookup
+ #--------------------------------------------------------
+ log_nH = np.log10( data["PartType0", "H_number_density"] )
+ log_T = np.log10( data["PartType0", "Temperature"] )
- # get name of owls_ion_file for given ion
- #--------------------------------------------------------
- owls_ion_path = self._get_owls_ion_data_dir()
- fname = os.path.join( owls_ion_path, ion+".hdf5" )
+ # get name of owls_ion_file for given ion
+ #--------------------------------------------------------
+ itab.set_iz( data.ds.current_redshift )
- # create ionization table for this redshift
- #--------------------------------------------------------
- itab = oit.IonTableOWLS( fname )
- itab.set_iz( data.ds.current_redshift )
-
- # find ion balance using log nH and log T
- #--------------------------------------------------------
- i_frac = itab.interp( log_nH, log_T )
- return data[ftype,"Density"] * m_frac * i_frac
-
- return _ion_density
+ # find ion balance using log nH and log T
+ #--------------------------------------------------------
+ i_frac = itab.interp( log_nH, log_T )
+ return data[ftype,"Density"] * m_frac * i_frac
+ return _func
+
+ ion_path = self._get_owls_ion_data_dir()
+ fname = os.path.join( ion_path, ion+".hdf5" )
+ itab = oit.IonTableOWLS( fname )
+ return get_owls_ion_density_field(ion, ftype, itab)
diff -r 3f50767dec553868bf966a8e027f308f545e84c2 -r 2d7b93e59428f2dc6d4cde5030cacc4b113c6ea4 yt/frontends/owls/owls_ion_tables.py
--- a/yt/frontends/owls/owls_ion_tables.py
+++ b/yt/frontends/owls/owls_ion_tables.py
@@ -1,8 +1,8 @@
-"""
+"""
OWLS ion tables
A module to handle the HM01 UV background spectra and ionization data from the
-OWLS photoionization equilibrium lookup tables.
+OWLS photoionization equilibrium lookup tables.
@@ -17,27 +17,28 @@
#-----------------------------------------------------------------------------
from yt.utilities.on_demand_imports import _h5py as h5py
+import yt.extern.six as six
import numpy as np
-def h5rd( fname, path, dtype=None ):
+def h5rd(fname, path, dtype=None):
""" Read Data. Return a dataset located at <path> in file <fname> as
- a numpy array.
+ a numpy array.
e.g. rd( fname, '/PartType0/Coordinates' ). """
data = None
- with h5py.File( fname, 'r' ) as h5f:
- ds = h5f[path]
- if dtype is None:
- dtype = ds.dtype
- data = np.zeros( ds.shape, dtype=dtype )
- data = ds.value
+ fid = h5py.h5f.open(six.b(fname), h5py.h5f.ACC_RDONLY)
+ dg = h5py.h5d.open(fid, path.encode('ascii'))
+ if dtype is None:
+ dtype = dg.dtype
+ data = np.zeros(dg.shape, dtype=dtype)
+ dg.read(h5py.h5s.ALL, h5py.h5s.ALL, data)
+ fid.close()
return data
-
class IonTableSpectrum:
""" A class to handle the HM01 spectra in the OWLS ionization tables. """
@@ -45,17 +46,16 @@
def __init__(self, ion_file):
where = '/header/spectrum/gammahi'
- self.GH1 = h5rd( ion_file, where ) # GH1[1/s]
+ self.GH1 = h5rd(ion_file, where) # GH1[1/s]
where = '/header/spectrum/logenergy_ryd'
- self.logryd = h5rd( ion_file, where ) # E[ryd]
+ self.logryd = h5rd(ion_file, where) # E[ryd]
where = '/header/spectrum/logflux'
- self.logflux = h5rd( ion_file, where ) # J[ergs/s/Hz/Sr/cm^2]
+ self.logflux = h5rd(ion_file, where) # J[ergs/s/Hz/Sr/cm^2]
where = '/header/spectrum/redshift'
- self.z = h5rd( ion_file, where ) # z
-
+ self.z = h5rd(ion_file, where) # z
def return_table_GH1_at_z(self,z):
@@ -68,9 +68,9 @@
else:
i_zhi = i_zlo
i_zlo = i_zlo - 1
-
+
z_frac = (z - self.z[i_zlo]) / (self.z[i_zhi] - self.z[i_zlo])
-
+
# find GH1 from table
#-----------------------------------------------------------------
logGH1_all = np.log10( self.GH1 )
@@ -80,8 +80,6 @@
GH1_table = 10.0**logGH1_table
return GH1_table
-
-
class IonTableOWLS:
@@ -90,7 +88,7 @@
DELTA_nH = 0.25
DELTA_T = 0.1
-
+
def __init__(self, ion_file):
self.ion_file = ion_file
@@ -104,13 +102,13 @@
# read the ionization fractions
# linear values stored in file so take log here
- # ionbal is the ionization balance (i.e. fraction)
+ # ionbal is the ionization balance (i.e. fraction)
#---------------------------------------------------------------
- self.ionbal = h5rd( ion_file, '/ionbal' ).astype(np.float64)
+ self.ionbal = h5rd( ion_file, '/ionbal' ).astype(np.float64)
self.ionbal_orig = self.ionbal.copy()
- ipositive = np.where( self.ionbal > 0.0 )
- izero = np.where( self.ionbal <= 0.0 )
+ ipositive = self.ionbal > 0.0
+ izero = np.logical_not(ipositive)
self.ionbal[izero] = self.ionbal[ipositive].min()
self.ionbal = np.log10( self.ionbal )
@@ -118,7 +116,7 @@
# load in background spectrum
#---------------------------------------------------------------
- self.spectrum = IonTableSpectrum( ion_file )
+ self.spectrum = IonTableSpectrum( ion_file )
# calculate the spacing along each dimension
#---------------------------------------------------------------
@@ -129,9 +127,6 @@
self.order_str = '[log nH, log T, z]'
-
-
-
# sets iz and fz
#-----------------------------------------------------
def set_iz( self, z ):
@@ -149,11 +144,11 @@
self.fz = ( z - self.z[iz] ) / self.dz[iz]
break
-
+
# interpolate the table at a fixed redshift for the input
- # values of nH and T ( input should be log ). A simple
- # tri-linear interpolation is used.
+ # values of nH and T ( input should be log ). A simple
+ # tri-linear interpolation is used.
#-----------------------------------------------------
def interp( self, nH, T ):
@@ -162,7 +157,7 @@
if nH.size != T.size:
raise ValueError(' owls_ion_tables: array size mismatch !!! ')
-
+
# field discovery will have nH.size == 1 and T.size == 1
# in that case we simply return 1.0
@@ -185,14 +180,14 @@
x_T_clip = np.clip( x_T, 0.0, self.T.size-1.001 )
fT,iT = np.modf( x_T_clip )
iT = iT.astype( np.int32 )
-
+
# short names for previously calculated iz and fz
#-----------------------------------------------------
iz = self.iz
fz = self.fz
-
+
# calculate interpolated value
# use tri-linear interpolation on the log values
#-----------------------------------------------------
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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-svn-spacepope.org/attachments/20160429/3493d886/attachment.htm>
More information about the yt-svn
mailing list