[Yt-svn] yt: 2 new changesets
hg at spacepope.org
hg at spacepope.org
Thu Oct 28 18:15:01 PDT 2010
hg Repository: yt
details: yt/rev/57f2fabaebf2
changeset: 3492:57f2fabaebf2
user: Matthew Turk <matthewturk at gmail.com>
date:
Thu Oct 28 10:09:47 2010 -0700
description:
Adding subregions for exporting data to Sunrise
hg Repository: yt
details: yt/rev/782bda7e39ad
changeset: 3493:782bda7e39ad
user: Matthew Turk <matthewturk at gmail.com>
date:
Thu Oct 28 18:14:55 2010 -0700
description:
Attempting first pass at units. Temperature is calculated incorrectly, but I
believe Density and Gas_Energy are correct.
diffstat:
yt/analysis_modules/sunrise_export/sunrise_exporter.py | 37 ++++++++--
yt/frontends/art/data_structures.py | 60 +++++++++++++----
yt/frontends/art/fields.py | 19 +++++
yt/frontends/art/io.py | 2 +-
4 files changed, 92 insertions(+), 26 deletions(-)
diffs (235 lines):
diff -r a4fb580b735c -r 782bda7e39ad yt/analysis_modules/sunrise_export/sunrise_exporter.py
--- a/yt/analysis_modules/sunrise_export/sunrise_exporter.py Wed Oct 27 23:34:28 2010 -0700
+++ b/yt/analysis_modules/sunrise_export/sunrise_exporter.py Thu Oct 28 18:14:55 2010 -0700
@@ -30,11 +30,13 @@
# We silently fail here
pass
+import time
+import numpy as na
+
from yt.funcs import *
-import numpy as na
import yt.utilities.amr_utils as amr_utils
-def export_to_sunrise(pf, fn, write_particles = True):
+def export_to_sunrise(pf, fn, write_particles = True, subregion_bounds = None):
r"""Convert the contents of a dataset to a FITS file format that Sunrise
understands.
@@ -53,6 +55,12 @@
The filename of the FITS file.
write_particles : bool, default is True
Whether to write out the star particles or not
+ subregion_bounds : list of tuples
+ This is a list of tuples describing the subregion of the top grid to
+ export. This will only work when only *one* root grid exists.
+ It is of the format:
+ [ (start_index_x, nx), (start_index_y, ny), (start_index_z, nz) ]
+ where nx, ny, nz are the number of cells to extract.
Notes
-----
@@ -109,7 +117,7 @@
output, refined = generate_flat_octree(pf,
["CellMassMsun","Temperature", "Metal_Density",
- "CellVolumeCode"])
+ "CellVolumeCode"], subregion_bounds = subregion_bounds)
cvcgs = output["CellVolumeCode"].astype('float64') * pf['cm']**3.0
# First the structure
@@ -204,7 +212,7 @@
ogl = amr_utils.OctreeGridList(grids)
return ogl, levels_finest, levels_all
-def generate_flat_octree(pf, fields):
+def generate_flat_octree(pf, fields, subregion_bounds = None):
"""
Generates two arrays, one of the actual values in a depth-first flat
octree array, and the other of the values describing the refinement.
@@ -218,14 +226,25 @@
output = na.zeros((o_length,len(fields)), dtype='float64')
refined = na.zeros(r_length, dtype='int32')
position = amr_utils.position()
- amr_utils.RecurseOctreeDepthFirst(0, 0, 0,
- ogl[0].dimensions[0],
- ogl[0].dimensions[1],
- ogl[0].dimensions[2],
+ if subregion_bounds is None:
+ sx, sy, sz = 0, 0, 0
+ nx, ny, nz = ogl[0].dimensions
+ else:
+ ss, ns = zip(*subregion_bounds)
+ sx, sy, sz = ss
+ nx, ny, nz = ns
+ print "Running from %s for %s cells" % (
+ (sx,sy,sz), (nx,ny,nz))
+ t1 = time.time()
+ amr_utils.RecurseOctreeDepthFirst(
+ sx, sy, sz, nx, ny, nz,
position, ogl[0].offset,
output, refined, ogl)
+ t2 = time.time()
+ print "Finished. Took %0.3e seconds." % (t2-t1)
dd = {}
- for i,field in enumerate(fields): dd[field] = output[:,i]
+ for i, field in enumerate(fields):
+ dd[field] = output[:position.output_pos,i]
return dd, refined
def generate_levels_octree(pf, fields):
diff -r a4fb580b735c -r 782bda7e39ad yt/frontends/art/data_structures.py
--- a/yt/frontends/art/data_structures.py Wed Oct 27 23:34:28 2010 -0700
+++ b/yt/frontends/art/data_structures.py Thu Oct 28 18:14:55 2010 -0700
@@ -47,6 +47,9 @@
import yt.frontends.ramses._ramses_reader as _ramses_reader
+from yt.utilities.physical_constants import \
+ mass_hydrogen_cgs
+
def num_deep_inc(f):
def wrap(self, *args, **kwargs):
self.num_deep += 1
@@ -120,10 +123,11 @@
pass
def _detect_fields(self):
- self.field_list = [ 'Density','Gas_Energy',
+ # This will need to be generalized to be used elsewhere.
+ self.field_list = [ 'Density','Total_Energy',
'x-momentum','y-momentum','z-momentum',
- 'Pressure','Gamma','Total_Energy','Potential_New'
- 'Potential_Old']
+ 'Pressure','Gamma','Gas_Energy',
+ 'Metal_Density1', 'Metal_Density2']
def _setup_classes(self):
dd = self._get_data_reader_dict()
@@ -389,24 +393,41 @@
self.time_units = {}
if len(self.parameters) == 0:
self._parse_parameter_file()
- self._setup_nounits_units()
self.conversion_factors = defaultdict(lambda: 1.0)
- self.time_units['1'] = 1
self.units['1'] = 1.0
self.units['unitary'] = 1.0 / (self.domain_right_edge - self.domain_left_edge).max()
- seconds = 1 #self["Time"]
+
+ z = self.current_redshift
+ h = self.hubble_constant
+ boxcm_cal = self["boxh"]
+ boxcm_uncal = boxcm_cal / h
+ box_proper = boxcm_uncal/(1+z)
+ for unit in mpc_conversion:
+ self.units[unit] = mpc_conversion[unit] * box_proper
+ self.units[unit+'h'] = mpc_conversion[unit] * box_proper * h
+ self.units[unit+'cm'] = mpc_conversion[unit] * boxcm_uncal
+ self.units[unit+'hcm'] = mpc_conversion[unit] * boxcm_cal
+ # Variable names have been chosen to reflect primary reference
+ Om0 = self["Om0"]
+ boxh = self["boxh"]
+ aexpn = self["aexpn"]
+ wmu = self["wmu"]
+ ng = self.domain_dimensions[0]
+ r0 = self["cmh"]/ng # comoving cm h^-1
+ t0 = 6.17e17/(self.hubble_constant + na.sqrt(self.omega_matter))
+ v0 = r0 / t0
+ rho0 = 1.8791e-29 * self.hubble_constant**2.0 * self.omega_matter
+ e0 = v0**2.0
+ self.conversion_factors["Density"] = rho0 * aexpn**-3
+ self.conversion_factors["Gas_Energy"] = e0/aexpn**2
+ # Now our conversion factors
+ for ax in 'xyz':
+ # Add on the 1e5 to get to cm/s
+ self.conversion_factors["%s-velocity" % ax] = v0/aexpn
+ seconds = t0
self.time_units['years'] = seconds / (365*3600*24.0)
self.time_units['days'] = seconds / (3600*24.0)
- def _setup_nounits_units(self):
- z = 0
- mylog.warning("Setting 1.0 in code units to be 1.0 cm")
- if not self.has_key("TimeUnits"):
- mylog.warning("No time units. Setting 1.0 = 1 second.")
- self.conversion_factors["Time"] = 1.0
- for unit in mpc_conversion.keys():
- self.units[unit] = mpc_conversion[unit] / mpc_conversion["cm"]
-
def _parse_parameter_file(self):
# We set our domain to run from 0 .. 1 since we are otherwise
# unconstrained.
@@ -467,10 +488,17 @@
header_vals[name] = output
self.dimensionality = 3 # We only support three
self.refine_by = 2 # Octree
+ # Update our parameters with the header and with some compile-time
+ # constants we will set permanently.
+ self.parameters.update(header_vals)
+ self.parameters["Y_p"] = 0.245
+ self.parameters["wmu"] = 4.0/(8.0-5.0*self.parameters["Y_p"])
+ self.parameters["gamma"] = 5./3.
+ self.current_redshift = self.parameters["aexpn"]**-1.0 - 1.0
self.data_comment = header_vals['jname']
self.current_time = header_vals['t']
self.omega_lambda = header_vals['Oml0']
- self.omega_matter = header_vals['Om0'] - header_vals['Oml0']
+ self.omega_matter = header_vals['Om0']
self.hubble_constant = header_vals['hubble']
self.min_level = header_vals['min_level']
self.max_level = header_vals['max_level']
diff -r a4fb580b735c -r 782bda7e39ad yt/frontends/art/fields.py
--- a/yt/frontends/art/fields.py Wed Oct 27 23:34:28 2010 -0700
+++ b/yt/frontends/art/fields.py Thu Oct 28 18:14:55 2010 -0700
@@ -31,6 +31,8 @@
ValidateSpatial, \
ValidateGridType
import yt.data_objects.universal_fields
+from yt.utilities.physical_constants import \
+ boltzmann_constant_cgs, mass_hydrogen_cgs
class ARTFieldContainer(CodeFieldInfoContainer):
_shared_state = {}
@@ -46,6 +48,8 @@
"z-velocity":"velocity_z",
"Pressure":"pressure",
"Metallicity":"metallicity",
+ "Gas_Energy":"Gas_Energy",
+ "Total_Energy":"Total_Energy",
}
def _generate_translation(mine, theirs):
@@ -58,4 +62,19 @@
#print "Setting up translator from %s to %s" % (v, f)
_generate_translation(v, f)
+def _convertDensity(data):
+ return data.convert("Density")
+ARTFieldInfo["Density"]._units = r"\rm{g}/\rm{cm}^3"
+ARTFieldInfo["Density"]._projected_units = r"\rm{g}/\rm{cm}^2"
+ARTFieldInfo["Density"]._convert_function=_convertDensity
+def _convertEnergy(data):
+ return data.convert("Gas_Energy")
+ARTFieldInfo["Gas_Energy"]._units = r"\rm{ergs}/\rm{g}"
+ARTFieldInfo["Gas_Energy"]._convert_function=_convertEnergy
+
+def _Temperature(field, data):
+ tr = (data["Gamma"] - 1.0) * data["Gas_Energy"]
+ tr *= mass_hydrogen_cgs * data.pf["wmu"] / boltzmann_constant_cgs
+ return tr
+add_field("Temperature", function=_Temperature, units = r"\mathrm{K}")
diff -r a4fb580b735c -r 782bda7e39ad yt/frontends/art/io.py
--- a/yt/frontends/art/io.py Wed Oct 27 23:34:28 2010 -0700
+++ b/yt/frontends/art/io.py Thu Oct 28 18:14:55 2010 -0700
@@ -78,7 +78,7 @@
if grid.Level == 0: # We only have one root grid
self.preload_level(0)
tr = self.level_data[0][field_id,:].reshape(
- pf.domain_dimensions, order="F")
+ pf.domain_dimensions, order="F").copy()
return tr.swapaxes(0, 2)
tr = na.zeros(grid.ActiveDimensions, dtype='float64')
filled = na.zeros(grid.ActiveDimensions, dtype='int32')
More information about the yt-svn
mailing list