[Yt-svn] yt-commit r1421 - trunk/yt/extensions
britton at wrangler.dreamhost.com
britton at wrangler.dreamhost.com
Wed Sep 2 10:06:24 PDT 2009
Author: britton
Date: Wed Sep 2 10:06:23 2009
New Revision: 1421
URL: http://yt.spacepope.org/changeset/1421
Log:
Changed method names to be closer to PEP 8 style and added a Cosmology
calculator attribute. Also added _create_cosmology_splice method as the
first part of calculating light cone solutions that simply involves connecting
dataset together to get from one redshift to another.
Modified:
trunk/yt/extensions/EnzoSimulation.py
Modified: trunk/yt/extensions/EnzoSimulation.py
==============================================================================
--- trunk/yt/extensions/EnzoSimulation.py (original)
+++ trunk/yt/extensions/EnzoSimulation.py Wed Sep 2 10:06:23 2009
@@ -53,20 +53,29 @@
# Convert initial/final redshifts to times.
if self.enzoParameters['ComovingCoordinates']:
- self.cosmology = lagos.EnzoCosmology(HubbleConstantNow =
+ # Instantiate a cosmology calculator.
+ self.cosmology = lagos.Cosmology(HubbleConstantNow =
+ (100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
+ OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
+ OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'])
+
+ # Instantiate EnzoCosmology object for units and time conversions.
+ self.enzo_cosmology = lagos.EnzoCosmology(HubbleConstantNow =
(100.0 * self.enzoParameters['CosmologyHubbleConstantNow']),
OmegaMatterNow = self.enzoParameters['CosmologyOmegaMatterNow'],
OmegaLambdaNow = self.enzoParameters['CosmologyOmegaLambdaNow'],
InitialRedshift = self.enzoParameters['CosmologyInitialRedshift'])
if self.InitialRedshift is not None:
- self.InitialTime = self.cosmology.ComputeTimeFromRedshift(self.InitialRedshift) / self.cosmology.TimeUnits
+ self.InitialTime = self.enzo_cosmology.ComputeTimeFromRedshift(self.InitialRedshift) / \
+ self.enzo_cosmology.TimeUnits
if self.FinalRedshift is not None:
- self.FinalTime = self.cosmology.ComputeTimeFromRedshift(self.FinalRedshift) / self.cosmology.TimeUnits
+ self.FinalTime = self.enzo_cosmology.ComputeTimeFromRedshift(self.FinalRedshift) / \
+ self.enzo_cosmology.TimeUnits
# Get initial time of simulation.
if self.enzoParameters['ComovingCoordinates'] and \
self.enzoParameters.has_key('CosmologyInitialRedshift'):
- self.SimulationInitialTime = self.cosmology.InitialTime / self.cosmology.TimeUnits
+ self.SimulationInitialTime = self.enzo_cosmology.InitialTime / self.enzo_cosmology.TimeUnits
elif self.enzoParameters.has_key('InitialTime'):
self.SimulationInitialTime = self.enzoParameters['InitialTime']
else:
@@ -87,7 +96,8 @@
"Calculates time from redshift of redshift dumps."
for output in self.redshiftOutputs:
- output['time'] = self.cosmology.ComputeTimeFromRedshift(output['redshift']) / self.cosmology.TimeUnits
+ output['time'] = self.enzo_cosmology.ComputeTimeFromRedshift(output['redshift']) / \
+ self.enzo_cosmology.TimeUnits
def _CalculateTimeDumps(self):
"Calculates time dumps and their redshifts if cosmological."
@@ -102,8 +112,8 @@
self.timeOutputs.append({'index':index,'filename':filename,'time':current_time})
if self.enzoParameters['ComovingCoordinates']:
- t = self.cosmology.InitialTime + (self.enzoParameters['dtDataDump'] * self.cosmology.TimeUnits * index)
- self.timeOutputs[-1]['redshift'] = self.cosmology.ComputeRedshiftFromTime(t)
+ t = self.enzo_cosmology.InitialTime + (self.enzoParameters['dtDataDump'] * self.enzo_cosmology.TimeUnits * index)
+ self.timeOutputs[-1]['redshift'] = self.enzo_cosmology.ComputeRedshiftFromTime(t)
current_time += self.enzoParameters['dtDataDump']
index += 1
@@ -191,6 +201,139 @@
self.enzoParameters['DataDumpDir'] = "DD"
self.enzoParameters['ComovingCoordinates'] = 0
+ def _create_cosmology_splice(self, minimal=True, deltaz_min=0.0):
+ "Create list of datasets to be used for LightCones or LightRays."
+
+ # Calculate maximum delta z for each data dump.
+ self._calculate_deltaz_max()
+
+ # Calculate minimum delta z for each data dump.
+ self._calculate_deltaz_min(deltaz_min=deltaz_min)
+
+ cosmology_splice = []
+
+ # Use minimum number of datasets to go from z_i to z_f.
+ if minimal:
+
+ z_Tolerance = 1e-4
+ z = self.InitialRedshift
+
+ # fill redshift space with datasets
+ while ((z > self.FinalRedshift) and
+ (na.fabs(z - self.FinalRedshift) > z_Tolerance)):
+ # Sort data outputs by proximity to current redsfhit.
+ self.allOutputs.sort(key=lambda obj:na.fabs(z - obj['redshift']))
+ # For first data dump, choose closest to desired redshift.
+ if (len(cosmology_splice) == 0):
+ cosmology_splice.append(self.allOutputs[0])
+ # Start with data dump closest to desired redshift and move backward
+ # until one is within max delta z of last output in solution list.
+ else:
+ output = self.allOutputs[0]
+ while (z > output['redshift']):
+ output = output['previous']
+ if (output is None):
+ mylog.error("CalculateLightRaySolution: search for data output went off the end of the stack.")
+ mylog.error("Could not calculate light ray solution.")
+ return
+ if (output['redshift'] == cosmology_splice[-1]['redshift']):
+ mylog.error("CalculateLightRaySolution: No data dump between z = %f and %f." % \
+ ((cosmology_splice[-1]['redshift'] - cosmology_splice[-1]['deltazMax']),
+ cosmology_splice[-1]['redshift']))
+ mylog.error("Could not calculate light ray solution.")
+ return
+ cosmology_splice.append(output)
+ z = cosmology_splice[-1]['redshift'] - cosmology_splice[-1]['deltazMax']
+
+ # Make light ray using maximum number of datasets (minimum spacing).
+ else:
+ # Sort data outputs by proximity to current redsfhit.
+ self.allOutputs.sort(key=lambda obj:na.fabs(self.InitialRedshift - obj['redshift']))
+ # For first data dump, choose closest to desired redshift.
+ cosmology_splice.append(self.allOutputs[0])
+
+ nextOutput = cosmology_splice[-1]['next']
+ while (nextOutput is not None):
+ if (nextOutput['redshift'] <= self.FinalRedshift):
+ break
+ if ((cosmology_splice[-1]['redshift'] - nextOutput['redshift']) > cosmology_splice[-1]['deltazMin']):
+ cosmology_splice.append(nextOutput)
+ nextOutput = nextOutput['next']
+
+ mylog.info("create_cosmology_splice: Used %d data dumps to get from z = %f to %f." %
+ (len(cosmology_splice),self.InitialRedshift,self.FinalRedshift))
+
+ return cosmology_splice
+
+ def _calculate_deltaz_max(self):
+ "Calculate delta z that corresponds to full box length going from z to (z - delta z)."
+
+ d_Tolerance = 1e-4
+ max_Iterations = 100
+
+ targetDistance = self.enzoParameters['CosmologyComovingBoxSize']
+
+ for output in self.allOutputs:
+ z = output['redshift']
+
+ # Calculate delta z that corresponds to the length of the box at a given redshift.
+ # Use Newton's method to calculate solution.
+ z1 = z
+ z2 = z1 - 0.1 # just an initial guess
+ distance1 = 0.0
+ iteration = 1
+
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+
+ while ((na.fabs(distance2-targetDistance)/distance2) > d_Tolerance):
+ m = (distance2 - distance1) / (z2 - z1)
+ z1 = z2
+ distance1 = distance2
+ z2 = ((targetDistance - distance2) / m) + z2
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+ iteration += 1
+ if (iteration > max_Iterations):
+ if self.verbose: mylog.error("calculate_deltaz_max: Warning - max iterations exceeded for z = %f (delta z = %f)." %
+ (z,na.fabs(z2-z)))
+ break
+ output['deltazMax'] = na.fabs(z2-z)
+
+ def _calculate_deltaz_min(self, deltaz_min=0.0):
+ "Calculate delta z that corresponds to a single top grid pixel going from z to (z - delta z)."
+
+ d_Tolerance = 1e-4
+ max_Iterations = 100
+
+ targetDistance = self.enzoParameters['CosmologyComovingBoxSize'] / self.enzoParameters['TopGridDimensions'][0]
+
+ for output in self.allOutputs:
+ z = output['redshift']
+
+ # Calculate delta z that corresponds to the length of a top grid pixel at a given redshift.
+ # Use Newton's method to calculate solution.
+ z1 = z
+ z2 = z1 - 0.01 # just an initial guess
+ distance1 = 0.0
+ iteration = 1
+
+ # Convert comoving radial distance into Mpc / h, since that's how box size is stored.
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+
+ while ((na.fabs(distance2 - targetDistance) / distance2) > d_Tolerance):
+ m = (distance2 - distance1) / (z2 - z1)
+ z1 = z2
+ distance1 = distance2
+ z2 = ((targetDistance - distance2) / m) + z2
+ distance2 = self.cosmology.ComovingRadialDistance(z2,z) * self.enzoParameters['CosmologyHubbleConstantNow']
+ iteration += 1
+ if (iteration > max_Iterations):
+ if self.verbose: mylog.error("calculate_deltaz_max: Warning - max iterations exceeded for z = %f (delta z = %f)." %
+ (z,na.fabs(z2-z)))
+ break
+ # Use this calculation or the absolute minimum specified by the user.
+ output['deltazMin'] = max(na.fabs(z2-z),deltaz_min)
+
EnzoParameterDict = {"CosmologyOmegaMatterNow": float,
"CosmologyOmegaLambdaNow": float,
"CosmologyHubbleConstantNow": float,
More information about the yt-svn
mailing list