[Yt-svn] yt: Fixed domain issues in the Sunrise exporter and a few small ...
hg at spacepope.org
hg at spacepope.org
Tue Nov 9 21:54:32 PST 2010
hg Repository: yt
details: yt/rev/a198f8022927
changeset: 3518:a198f8022927
user: Christopher Erick Moody <cemoody at ucsc.edu>
date:
Tue Nov 09 21:54:19 2010 -0800
description:
Fixed domain issues in the Sunrise exporter and a few small bugs. Updated the ART current time calculation, and renamed the metal fields to reasonable and intutive things.
diffstat:
yt/analysis_modules/sunrise_export/sunrise_exporter.py | 28 +++++++----
yt/frontends/art/data_structures.py | 19 +++++--
yt/frontends/art/fields.py | 43 +++++++++--------
3 files changed, 54 insertions(+), 36 deletions(-)
diffs (194 lines):
diff -r 4b7fca0be218 -r a198f8022927 yt/analysis_modules/sunrise_export/sunrise_exporter.py
--- a/yt/analysis_modules/sunrise_export/sunrise_exporter.py Tue Nov 09 17:07:24 2010 -0700
+++ b/yt/analysis_modules/sunrise_export/sunrise_exporter.py Tue Nov 09 21:54:19 2010 -0800
@@ -92,7 +92,12 @@
# output->addKey("logflux", true, "Column L_lambda values are log (L_lambda)");
col_list = []
- DLE, DRE = pf.domain_left_edge, pf.domain_right_edge
+ if subregion_bounds == None:
+ DLE, DRE = pf.domain_left_edge, pf.domain_right_edge
+ else:
+ DLE, DX = zip(*subregion_bounds)
+ DLE, DX = na.array(DLE), na.array(DX)
+ DRE = DLE + DX
reg = pf.h.region((DRE+DLE)/2.0, DLE, DRE)
if write_particles is True:
@@ -132,9 +137,9 @@
write_particles = True
def _MetalMass(field, data):
- return data["Metal_Density"] * data["CellVolume"]
+ return data["Metallicity"] * data["CellMassMsun"]
def _convMetalMass(data):
- return 1.0/1.989e33
+ return 1.0
add_field("MetalMass", function=_MetalMass,
convert_function=_convMetalMass)
@@ -167,9 +172,9 @@
st_table.header.update("hierarch lengthunit", "kpc", comment="Length unit for grid")
for i,a in enumerate('xyz'):
- st_table.header.update("min%s" % a, pf.domain_left_edge[i] * pf['kpc'])
- st_table.header.update("max%s" % a, pf.domain_right_edge[i] * pf['kpc'])
- st_table.header.update("n%s" % a, pf.domain_dimensions[i])
+ st_table.header.update("min%s" % a, DLE[i] * pf['kpc'])
+ st_table.header.update("max%s" % a, DRE[i] * pf['kpc'])
+ st_table.header.update("n%s" % a, DX[i])
st_table.header.update("subdiv%s" % a, 2)
st_table.header.update("subdivtp", "UNIFORM", "Type of grid subdivision")
@@ -194,12 +199,13 @@
col_list.append(pyfits.Column("mass_gas", format='D',
array=output.pop('CellMassMsun'), unit="Msun"))
col_list.append(pyfits.Column("mass_metals", format='D',
- array=output.pop('MetalMass')*cvcgs/1.989e33, unit="Msun"))
- # Unit is set to K but it's really K*Msun
+ array=output.pop('MetalMass'), unit="Msun"))
+ # The units for gas_temp are really K*Msun. For older Sunrise versions
+ # you must set the unit to just K
col_list.append(pyfits.Column("gas_temp_m", format='D',
- array=output['TemperatureTimesCellMassMsun'], unit="K"))
+ array=output['TemperatureTimesCellMassMsun'], unit="K*Msun"))
col_list.append(pyfits.Column("gas_teff_m", format='D',
- array=output.pop('TemperatureTimesCellMassMsun'), unit="K"))
+ array=output.pop('TemperatureTimesCellMassMsun'), unit="K*Msun"))
col_list.append(pyfits.Column("cell_volume", format='D',
array=output.pop('CellVolumeCode').astype('float64')*pf['kpc']**3.0,
unit="kpc^3"))
@@ -213,13 +219,13 @@
mg_table.name = "GRIDDATA"
# Add a dummy Primary; might be a better way to do this!
- hls = [pyfits.PrimaryHDU(), st_table, mg_table]
col_list = [pyfits.Column("dummy", format="F", array=na.zeros(1, dtype='float32'))]
cols = pyfits.ColDefs(col_list)
md_table = pyfits.new_table(cols)
md_table.header.update("snaptime", pf.current_time*pf["years"])
md_table.name = "YT"
+ hls = [pyfits.PrimaryHDU(), st_table, mg_table,md_table]
if write_particles: hls.append(pd_table)
hdus = pyfits.HDUList(hls)
hdus.writeto(fn, clobber=True)
diff -r 4b7fca0be218 -r a198f8022927 yt/frontends/art/data_structures.py
--- a/yt/frontends/art/data_structures.py Tue Nov 09 17:07:24 2010 -0700
+++ b/yt/frontends/art/data_structures.py Tue Nov 09 21:54:19 2010 -0800
@@ -127,7 +127,7 @@
self.field_list = [ 'Density','Total_Energy',
'x-momentum','y-momentum','z-momentum',
'Pressure','Gamma','Gas_Energy',
- 'Metal_Density1', 'Metal_Density2',
+ 'Metal_DensitySNII', 'Metal_DensitySNIa',
'Potential_New','Potential_Old']
def _setup_classes(self):
@@ -377,11 +377,11 @@
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"] = 'art'
self.parameters["Time"] = 1. # default unit is 1...
+ self.parameters["InitialTime"]=self.current_time
def __repr__(self):
return self.basename.rsplit(".", 1)[0]
@@ -438,9 +438,9 @@
self.rho0*self.v0**2*(aexpn**-5.0)
tr = self.tr
self.conversion_factors["Temperature"] = tr
- self.conversion_factors["Metal_Density"] = 1
- self.conversion_factors["Metal_Density1"] = 1
- self.conversion_factors["Metal_Density2"] = 1
+ self.conversion_factors["Metallicity"] = 1
+ self.conversion_factors["MetallicitySNII"] = 1
+ self.conversion_factors["MetallicitySNIa"] = 1
# Now our conversion factors
for ax in 'xyz':
@@ -527,6 +527,15 @@
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
+ self.hubble_time = 1.0/(self.hubble_constant*100/3.08568025e19)
+ #self.hubble_time /= 3.168876e7 #Gyr in s
+ def integrand(x,oml=self.omega_lambda,omb=self.omega_matter):
+ return 1./(x*na.sqrt(oml+omb*x**-3.0))
+ spacings = na.logspace(-5,na.log10(self.parameters['aexpn']),1e5)
+ integrand_arr = integrand(spacings)
+ self.current_time = na.trapz(integrand_arr,dx=na.diff(spacings))
+ self.current_time *= self.hubble_time
+
for to_skip in ['tl','dtl','tlold','dtlold','iSO']:
_skip_record(f)
diff -r 4b7fca0be218 -r a198f8022927 yt/frontends/art/fields.py
--- a/yt/frontends/art/fields.py Tue Nov 09 17:07:24 2010 -0700
+++ b/yt/frontends/art/fields.py Tue Nov 09 21:54:19 2010 -0800
@@ -93,36 +93,39 @@
ARTFieldInfo["Temperature"]._units = r"\mathrm{K}"
ARTFieldInfo["Temperature"]._convert_function=_convertTemperature
-def _Metal_DensitySNII(field, data):
- tr = data["Metal_Density1"] / data["Density"]
+def _MetallicitySNII(field, data):
+ #get the dimensionaless mass fraction
+ tr = data["Metal_DensitySNII"] / data["Density"]
tr *= data.pf.conversion_factors["Density"]
return tr
-def _convertMetal_DensitySNII(data):
- return data.convert("Metal_Density1")
+def _convert_MetallicitySNII(data):
+ return data.convert("MetallicitySNII")
-add_field("Metal_DensitySNII", function=_Temperature, units = r"\mathrm{K}")
-ARTFieldInfo["Metal_DensitySNII"]._units = r"\mathrm{K}"
-ARTFieldInfo["Metal_DensitySNII"]._convert_function=_convertMetal_DensitySNII
+add_field("MetallicitySNII", function=_MetallicitySNII, units = r"\mathrm{K}")
+ARTFieldInfo["MetallicitySNII"]._units = r"\mathrm{K}"
+ARTFieldInfo["MetallicitySNII"]._convert_function=_convert_MetallicitySNII
-def _Metal_DensitySNIa(field, data):
- tr = data["Metal_Density2"] / data["Density"]
+def _MetallicitySNIa(field, data):
+ #get the dimensionaless mass fraction
+ tr = data["Metal_DensitySNIa"] / data["Density"]
tr *= data.pf.conversion_factors["Density"]
return tr
-def _convertMetal_DensitySNIa(data):
- return data.convert("Metal_Density2")
+def _convert_MetallicitySNIa(data):
+ return data.convert("MetallicitySNIa")
-add_field("Metal_DensitySNIa", function=_Temperature, units = r"\mathrm{K}")
-ARTFieldInfo["Metal_DensitySNIa"]._units = r"\mathrm{K}"
-ARTFieldInfo["Metal_DensitySNIa"]._convert_function=_convertMetal_DensitySNIa
+add_field("MetallicitySNIa", function=_MetallicitySNIa, units = r"\mathrm{K}")
+ARTFieldInfo["MetallicitySNIa"]._units = r"\mathrm{K}"
+ARTFieldInfo["MetallicitySNIa"]._convert_function=_convert_MetallicitySNIa
-def _Metal_Density(field, data):
+def _Metallicity(field, data):
+ #get the dimensionaless mass fraction of the total metals
tr = data["Metal_Density2"] / data["Density"]
tr += data["Metal_Density1"] / data["Density"]
tr *= data.pf.conversion_factors["Density"]
return tr
-def _convertMetal_Density(data):
- return data.convert("Metal_Density")
+def _convert_Metallicity(data):
+ return data.convert("Metallicity")
-add_field("Metal_Density", function=_Temperature, units = r"\mathrm{K}")
-ARTFieldInfo["Metal_Density"]._units = r"\mathrm{K}"
-ARTFieldInfo["Metal_Density"]._convert_function=_convertMetal_Density
\ No newline at end of file
+add_field("Metallicity", function=_Metallicity, units = r"\mathrm{K}")
+ARTFieldInfo["Metallicity"]._units = r"\mathrm{K}"
+ARTFieldInfo["Metallicity"]._convert_function=_convert_Metallicity
\ No newline at end of file
More information about the yt-svn
mailing list