[Yt-svn] yt: Initial attempt at converting comoving units in usecosmology...
hg at spacepope.org
hg at spacepope.org
Tue Dec 7 13:28:09 PST 2010
hg Repository: yt
details: yt/rev/bdc4b6935e3f
changeset: 3578:bdc4b6935e3f
user: Matthew Turk <matthewturk at gmail.com>
date:
Tue Dec 07 13:28:05 2010 -0800
description:
Initial attempt at converting comoving units in usecosmology=1 FLASH simulations
diffstat:
yt/frontends/flash/data_structures.py | 53 ++++++++++++++++++++++++++++++++-
yt/frontends/flash/fields.py | 55 +++++++++++++++++++++++++++++++++-
2 files changed, 103 insertions(+), 5 deletions(-)
diffs (167 lines):
diff -r 672c309bd892 -r bdc4b6935e3f yt/frontends/flash/data_structures.py
--- a/yt/frontends/flash/data_structures.py Mon Dec 06 14:54:04 2010 -0800
+++ b/yt/frontends/flash/data_structures.py Tue Dec 07 13:28:05 2010 -0800
@@ -83,7 +83,7 @@
def _detect_fields(self):
ncomp = self._handle["/unknown names"].shape[0]
- self.field_list = [s.strip() for s in self._handle["/unknown names"][:].flat]
+ self.field_list = [s for s in self._handle["/unknown names"][:].flat]
if ("/particle names" in self._handle) :
self.field_list += ["particle_" + s[0].strip() for s
in self._handle["/particle names"][:]]
@@ -163,6 +163,17 @@
def _setup_derived_fields(self):
self.derived_field_list = []
+ for field in self.parameter_file.field_info:
+ try:
+ fd = self.parameter_file.field_info[field].get_dependencies(
+ pf = self.parameter_file)
+ except:
+ continue
+ available = na.all([f in self.field_list for f in fd.requested])
+ if available: self.derived_field_list.append(field)
+ for field in self.field_list:
+ if field not in self.derived_field_list:
+ self.derived_field_list.append(field)
def _setup_data_io(self):
self.io = io_registry[self.data_style](self.parameter_file)
@@ -200,8 +211,11 @@
self.time_units = {}
if len(self.parameters) == 0:
self._parse_parameter_file()
- self._setup_nounits_units()
self.conversion_factors = defaultdict(lambda: 1.0)
+ if self.cosmological_simulation == 1:
+ self._setup_comoving_units()
+ else:
+ self._setup_nounits_units()
self.time_units['1'] = 1
self.units['1'] = 1.0
self.units['unitary'] = 1.0 / \
@@ -212,7 +226,26 @@
for p, v in self._conversion_override.items():
self.conversion_factors[p] = v
+ def _setup_comoving_units(self):
+ self.conversion_factors['dens'] = (1.0 + self.current_redshift)**3.0
+ self.conversion_factors['pres'] = (1.0 + self.current_redshift)**1.0
+ self.conversion_factors['temp'] = (1.0 + self.current_redshift)**-2.0
+ self.conversion_factors['velx'] = (1.0 + self.current_redshift)
+ self.conversion_factors['vely'] = self.conversion_factors['velx']
+ self.conversion_factors['velz'] = self.conversion_factors['velx']
+ 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 _setup_nounits_units(self):
+ self.conversion_factors['dens'] = 1.0
+ self.conversion_factors['pres'] = 1.0
+ self.conversion_factors['temp'] = 1.0
+ self.conversion_factors['velx'] = 1.0
+ self.conversion_factors['vely'] = 1.0
+ self.conversion_factors['velz'] = 1.0
z = 0
mylog.warning("Setting 1.0 in code units to be 1.0 cm")
if not self.has_key("TimeUnits"):
@@ -255,6 +288,22 @@
else:
self.current_time = \
float(self._find_parameter("real", "time", scalar=True))
+
+ try:
+ use_cosmo = self._find_parameter("logical", "usecosmology")
+ except KeyError:
+ use_cosmo = 0
+
+ if use_cosmo == 1:
+ self.cosmological_simulation = 1
+ self.current_redshift = self._find_parameter("real", "redshift",
+ scalar = True)
+ self.omega_lambda = self._find_parameter("real", "cosmologicalconstant")
+ self.omega_matter = self._find_parameter("real", "omegamatter")
+ self.hubble_constant = self._find_parameter("real", "hubbleconstant")
+ else:
+ self.current_redshift = self.omega_lambda = self.omega_matter = \
+ self.hubble_constant = self.cosmological_simulation = 0.0
self._handle.close()
@classmethod
diff -r 672c309bd892 -r bdc4b6935e3f yt/frontends/flash/fields.py
--- a/yt/frontends/flash/fields.py Mon Dec 06 14:54:04 2010 -0800
+++ b/yt/frontends/flash/fields.py Tue Dec 07 13:28:05 2010 -0800
@@ -67,13 +67,63 @@
"particle_position_x" : "particle_posx",
"particle_position_y" : "particle_posy",
"particle_position_z" : "particle_posz",
- }
+ "Electron_Fraction" : "elec",
+ "HI_Fraction" : "h ",
+ "HD_Fraction" : "hd ",
+ "HeI_Fraction": "hel ",
+ "HeII_Fraction": "hep ",
+ "HeIII_Fraction": "hepp",
+ "HM_Fraction": "hmin",
+ "HII_Fraction": "hp ",
+ "H2I_Fraction": "htwo",
+ "H2II_Fraction": "htwp",
+ "DI_Fraction": "deut",
+ "DII_Fraction": "dplu"}
+
+def _get_density(fname):
+ def _dens(field, data):
+ return data[fname] * data['Density']
+ return _dens
+
+for fn1, fn2 in translation_dict.items():
+ if fn1.endswith("_Fraction"):
+ add_field(fn1.split("_")[0] + "_Density",
+ function=_get_density(fn1), take_log=True)
+
+def _get_alias(alias):
+ def _alias(field, data):
+ return data[alias]
+ return _alias
def _generate_translation(mine, theirs):
pfield = theirs.startswith("particle")
- add_field(theirs, function=lambda a, b: b[mine], take_log=True,
+ add_field(theirs, function=_get_alias(mine), take_log=True,
particle_type = pfield)
+def _get_convert(fname):
+ def _conv(data):
+ return data.convert(fname)
+ return _conv
+
+add_field("dens", function=lambda a,b: None, take_log=True,
+ convert_function=_get_convert("dens"),
+ units=r"\rm{g}/\rm{cm}^3")
+add_field("xvel", function=lambda a,b: None, take_log=False,
+ convert_function=_get_convert("xvel"),
+ units=r"\rm{cm}/\rm{s}")
+add_field("yvel", function=lambda a,b: None, take_log=False,
+ convert_function=_get_convert("yvel"),
+ units=r"\rm{cm}/\rm{s}")
+add_field("zvel", function=lambda a,b: None, take_log=False,
+ convert_function=_get_convert("zvel"),
+ units=r"\rm{cm}/\rm{s}")
+add_field("temp", function=lambda a,b: None, take_log=True,
+ convert_function=_get_convert("temp"),
+ units=r"\rm{K}")
+add_field("pres", function=lambda a,b: None, take_log=True,
+ convert_function=_get_convert("pres"),
+ units=r"\rm{unknown}")
+
for f,v in translation_dict.items():
if v not in FLASHFieldInfo:
pfield = v.startswith("particle")
@@ -130,4 +180,3 @@
add_field("divb", function=lambda a,b: None, take_log=False,
validators = [ValidateDataField("divb")],
units = r"\rm{G}\/\rm{cm}")
-
More information about the yt-svn
mailing list