[yt-svn] commit/yt: ngoldbaum: Merged in brittonsmith/yt (pull request #2345)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Sep 7 11:06:29 PDT 2016
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/4bc8ace93fda/
Changeset: 4bc8ace93fda
Branch: yt
User: ngoldbaum
Date: 2016-09-07 18:06:03+00:00
Summary: Merged in brittonsmith/yt (pull request #2345)
LightRay enhancements and bugfix (closes Issue #1258)
Affected #: 12 files
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b doc/source/analyzing/analysis_modules/cosmology_calculator.rst
--- a/doc/source/analyzing/analysis_modules/cosmology_calculator.rst
+++ b/doc/source/analyzing/analysis_modules/cosmology_calculator.rst
@@ -31,13 +31,13 @@
print("hubble distance", co.hubble_distance())
# distance from z = 0 to 0.5
- print("comoving radial distance", co.comoving_radial_distance(0, 0.5).in_units("Mpc/h"))
+ print("comoving radial distance", co.comoving_radial_distance(0, 0.5).in_units("Mpccm/h"))
# transverse distance
- print("transverse distance", co.comoving_transverse_distance(0, 0.5).in_units("Mpc/h"))
+ print("transverse distance", co.comoving_transverse_distance(0, 0.5).in_units("Mpccm/h"))
# comoving volume
- print("comoving volume", co.comoving_volume(0, 0.5).in_units("Gpc**3"))
+ print("comoving volume", co.comoving_volume(0, 0.5).in_units("Gpccm**3"))
# angulare diameter distance
print("angular diameter distance", co.angular_diameter_distance(0, 0.5).in_units("Mpc/h"))
@@ -67,7 +67,16 @@
# convert redshift to time after Big Bang (same as Hubble time)
print("t from z", co.t_from_z(0.5).in_units("Gyr"))
-Note, that all distances returned are comoving distances. All of the above
+.. warning::
+
+ Cosmological distance calculations return values that are either
+ in the comoving or proper frame, depending on the specific quantity. For
+ simplicity, the proper and comoving frames are set equal to each other
+ within the cosmology calculator. This means that for some distance value,
+ x, x.to("Mpc") and x.to("Mpccm") will be the same. The user should take
+ care to understand which reference frame is correct for the given calculation.
+
+All of the above
functions accept scalar values and arrays. The helper functions, `co.quan`
and `co.arr` exist to create unitful `YTQuantities` and `YTArray` with the
unit registry of the cosmology calculator. For more information on the usage
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b doc/source/analyzing/analysis_modules/light_ray_generator.rst
--- a/doc/source/analyzing/analysis_modules/light_ray_generator.rst
+++ b/doc/source/analyzing/analysis_modules/light_ray_generator.rst
@@ -49,13 +49,18 @@
* ``deltaz_min`` (*float*): Specifies the minimum Delta-z between
consecutive datasets in the returned list. Default: 0.0.
-* ``minimum_coherent_box_fraction`` (*float*): Used with
- ``use_minimum_datasets`` set to False, this parameter specifies the
- fraction of the total box size to be traversed before rerandomizing the
- projection axis and center. This was invented to allow light rays with
- thin slices to sample coherent large scale structure, but in practice
- does not work so well. Try setting this parameter to 1 and see what
- happens. Default: 0.0.
+* ``max_box_fraction`` (*float*): In terms of the size of the domain, the
+ maximum length a light ray segment can be in order to span the redshift interval
+ from one dataset to another. If using a zoom-in simulation, this parameter can
+ be set to the length of the high resolution region so as to limit ray segments
+ to that size. If the high resolution region is not cubical, the smallest side
+ should be used. Default: 1.0 (the size of the box)
+
+* ``minimum_coherent_box_fraction`` (*float*): Use to specify the minimum
+ length of a ray, in terms of the size of the domain, before the trajectory
+ is re-randomized. Set to 0 to have ray trajectory randomized for every
+ dataset. Set to np.inf (infinity) to use a single trajectory for the
+ entire ray. Default: 0.0.
* ``time_data`` (*bool*): Whether or not to include time outputs when
gathering datasets for time series. Default: True.
@@ -67,7 +72,7 @@
---------------------
Once the LightRay object has been instantiated, the
-:func:`~yt.analysis_modules.cosmological_observation.light_ray.light_ray.LightRay,make_light_ray`
+:func:`~yt.analysis_modules.cosmological_observation.light_ray.light_ray.LightRay.make_light_ray`
function will trace out the rays in each dataset and collect information for all the
fields requested. The output file will be an HDF5 file containing all the
cell field values for all the cells that were intersected by the ray. A
@@ -85,6 +90,21 @@
* ``seed`` (*int*): Seed for the random number generator. Default: None.
+* ``periodic`` (*bool*): If True, ray trajectories will make use of periodic
+ boundaries. If False, ray trajectories will not be periodic. Default : True.
+
+* ``left_edge`` (iterable of *floats* or *YTArray*): The left corner of the
+ region in which rays are to be generated. If None, the left edge will be
+ that of the domain. Default: None.
+
+* ``right_edge`` (iterable of *floats* or *YTArray*): The right corner of
+ the region in which rays are to be generated. If None, the right edge
+ will be that of the domain. Default: None.
+
+* ``min_level`` (*int*): The minimum refinement level of the spatial region in
+ which the ray passes. This can be used with zoom-in simulations where the
+ high resolution region does not keep a constant geometry. Default: None.
+
* ``start_position`` (*list* of floats): Used only if creating a light ray
from a single dataset. The coordinates of the starting position of the
ray. Default: None.
@@ -122,7 +142,82 @@
slice and 1 to have all processors work together on each projection.
Default: 1
-.. note:: As of :code:`yt-3.0`, the functionality for recording properties of the nearest halo to each element of the ray no longer exists. This is still available in :code:`yt-2.x`. If you would like to use this feature in :code:`yt-3.x`, help is needed to port it over. Contact the yt-users mailing list if you are interested in doing this.
+Useful Tips for Making LightRays
+--------------------------------
+
+Below are some tips that may come in handy for creating proper LightRays.
+
+How many snapshots do I need?
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The number of snapshots required to traverse some redshift interval depends
+on the simulation box size and cosmological parameters. Before running an
+expensive simulation only to find out that you don't have enough outputs
+to span the redshift interval you want, have a look at
+:ref:`planning-cosmology-simulations`. The functionality described there
+will allow you to calculate the precise number of snapshots and specific
+redshifts at which they should be written.
+
+My snapshots are too far apart!
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+The ``max_box_fraction`` keyword, provided when creating the `Lightray`,
+allows the user to control how long a ray segment can be for an
+individual dataset. Be default, the `LightRay` generator will try to
+make segments no longer than the size of the box to avoid sampling the
+same structures more than once. However, this can be increased in the
+case that the redshift interval between datasets is longer than the
+box size. Increasing this value should be done with caution as longer
+ray segments run a greater risk of coming back to somewhere near their
+original position.
+
+What if I have a zoom-in simulation?
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+A zoom-in simulation has a high resolution region embedded within a
+larger, low resolution volume. In this type of simulation, it is likely
+that you will want the ray segments to stay within the high resolution
+region. To do this, you must first specify the size of the high
+resolution region when creating the `LightRay` using the
+``max_box_fraction`` keyword. This will make sure that
+the calculation of the spacing of the segment datasets only takes into
+account the high resolution region and not the full box size. If your
+high resolution region is not a perfect cube, specify the smallest side.
+Then, in the call to
+:func:`~yt.analysis_modules.cosmological_observation.light_ray.light_ray.LightRay.make_light_ray`,
+use the ``left_edge`` and ``right_edge`` keyword arguments to specify the
+precise location of the high resolution region.
+
+Technically speaking, the ray segments should no longer be periodic
+since the high resolution region is only a sub-volume within the
+larger domain. To make the ray segments non-periodic, set the
+``periodic`` keyword to False. The LightRay generator will continue
+to generate randomly oriented segments until it finds one that fits
+entirely within the high resolution region. If you have a high
+resolution region that can move and change shape slightly as structure
+forms, use the `min_level` keyword to mandate that the ray segment only
+pass through cells that are refined to at least some minimum level.
+
+If the size of the high resolution region is not large enough to
+span the required redshift interval, the `LightRay` generator can
+be configured to treat the high resolution region as if it were
+periodic simply by setting the ``periodic`` keyword to True. This
+option should be used with caution as it will lead to the creation
+of disconnected ray segments within a single dataset.
+
+I want a continous trajectory over the entire ray.
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+Set the ``minimum_coherent_box_fraction`` keyword argument to a very
+large number, like infinity (`numpy.inf`).
+
+.. note::
+
+ As of :code:`yt-3.0`, the functionality for recording properties of
+ the nearest halo to each element of the ray no longer exists. This
+ is still available in :code:`yt-2.x`. If you would like to use this
+ feature in :code:`yt-3.x`, help is needed to port it over. Contact
+ the yt-users mailing list if you are interested in doing this.
What Can I do with this?
------------------------
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b doc/source/analyzing/analysis_modules/planning_cosmology_simulations.rst
--- a/doc/source/analyzing/analysis_modules/planning_cosmology_simulations.rst
+++ b/doc/source/analyzing/analysis_modules/planning_cosmology_simulations.rst
@@ -4,7 +4,7 @@
===================================================
If you want to run a cosmological simulation that will have just enough data
-outputs to create a cosmology splice, the
+outputs to create a light cone or light ray, the
:meth:`~yt.analysis_modules.cosmological_observation.cosmology_splice.CosmologySplice.plan_cosmology_splice`
function will calculate a list of redshifts outputs that will minimally
connect a redshift interval.
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -67,7 +67,7 @@
local_ytdata_000:
- yt/frontends/ytdata
- local_absorption_spectrum_004:
+ local_absorption_spectrum_005:
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_non_cosmo
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_non_cosmo_novpec
- yt/analysis_modules/absorption_spectrum/tests/test_absorption_spectrum.py:test_absorption_spectrum_cosmo
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/analysis_modules/cosmological_observation/cosmology_splice.py
--- a/yt/analysis_modules/cosmological_observation/cosmology_splice.py
+++ b/yt/analysis_modules/cosmological_observation/cosmology_splice.py
@@ -21,6 +21,8 @@
from yt.funcs import mylog
from yt.utilities.cosmology import \
Cosmology
+from yt.utilities.physical_constants import \
+ c
class CosmologySplice(object):
"""
@@ -67,7 +69,11 @@
max_box_fraction : float
In terms of the size of the domain, the maximum length a light
ray segment can be in order to span the redshift interval from
- one dataset to another.
+ one dataset to another. If using a zoom-in simulation, this
+ parameter can be set to the length of the high resolution
+ region so as to limit ray segments to that size. If the
+ high resolution region is not cubical, the smallest side
+ should be used.
Default: 1.0 (the size of the box)
deltaz_min : float
Specifies the minimum delta z between consecutive datasets
@@ -115,6 +121,7 @@
output['next'] = self.splice_outputs[i + 1]
# Calculate maximum delta z for each data dump.
+ self.max_box_fraction = max_box_fraction
self._calculate_deltaz_max()
# Calculate minimum delta z for each data dump.
@@ -144,7 +151,7 @@
self.splice_outputs.sort(key=lambda obj:np.fabs(z - obj['redshift']))
cosmology_splice.append(self.splice_outputs[0])
z = cosmology_splice[-1]["redshift"]
- z_target = z - max_box_fraction * cosmology_splice[-1]["dz_max"]
+ z_target = z - cosmology_splice[-1]["dz_max"]
# fill redshift space with datasets
while ((z_target > near_redshift) and
@@ -172,7 +179,7 @@
cosmology_splice.append(current_slice)
z = current_slice["redshift"]
- z_target = z - max_box_fraction * current_slice["dz_max"]
+ z_target = z - current_slice["dz_max"]
# Make light ray using maximum number of datasets (minimum spacing).
else:
@@ -199,8 +206,8 @@
mylog.info("create_cosmology_splice: Used %d data dumps to get from z = %f to %f." %
(len(cosmology_splice), far_redshift, near_redshift))
- # change the 'next' and 'previous' pointers to point to the correct outputs for the created
- # splice
+ # change the 'next' and 'previous' pointers to point to the correct outputs
+ # for the created splice
for i, output in enumerate(cosmology_splice):
if len(cosmology_splice) == 1:
output['previous'] = None
@@ -264,7 +271,8 @@
rounded += np.power(10.0, (-1.0*decimals))
z = rounded
- deltaz_max = self._deltaz_forward(z, self.simulation.box_size)
+ deltaz_max = self._deltaz_forward(z, self.simulation.box_size *
+ self.max_box_fraction)
outputs.append({'redshift': z, 'dz_max': deltaz_max})
z -= deltaz_max
@@ -282,72 +290,23 @@
from z to (z - delta z).
"""
- d_Tolerance = 1e-4
- max_Iterations = 100
+ target_distance = self.simulation.box_size * \
+ self.max_box_fraction
+ for output in self.splice_outputs:
+ output['dz_max'] = self._deltaz_forward(output['redshift'],
+ target_distance)
- target_distance = self.simulation.box_size
-
- for output in self.splice_outputs:
- z = output['redshift']
-
- # Calculate delta z that corresponds to the length of the box
- # at a given redshift using Newton's method.
- z1 = z
- z2 = z1 - 0.1 # just an initial guess
- distance1 = self.simulation.quan(0.0, "Mpccm / h")
- distance2 = self.cosmology.comoving_radial_distance(z2, z)
- iteration = 1
-
- while ((np.abs(distance2-target_distance)/distance2) > d_Tolerance):
- m = (distance2 - distance1) / (z2 - z1)
- z1 = z2
- distance1 = distance2
- z2 = ((target_distance - distance2) / m.in_units("Mpccm / h")) + z2
- distance2 = self.cosmology.comoving_radial_distance(z2, z)
- iteration += 1
- if (iteration > max_Iterations):
- mylog.error("calculate_deltaz_max: Warning - max iterations " +
- "exceeded for z = %f (delta z = %f)." %
- (z, np.abs(z2 - z)))
- break
- output['dz_max'] = np.abs(z2 - z)
-
def _calculate_deltaz_min(self, deltaz_min=0.0):
r"""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
-
target_distance = self.simulation.box_size / \
self.simulation.domain_dimensions[0]
-
for output in self.splice_outputs:
- z = output['redshift']
-
- # Calculate delta z that corresponds to the length of a
- # top grid pixel at a given redshift using Newton's method.
- z1 = z
- z2 = z1 - 0.01 # just an initial guess
- distance1 = self.simulation.quan(0.0, "Mpccm / h")
- distance2 = self.cosmology.comoving_radial_distance(z2, z)
- iteration = 1
-
- while ((np.abs(distance2 - target_distance) / distance2) > d_Tolerance):
- m = (distance2 - distance1) / (z2 - z1)
- z1 = z2
- distance1 = distance2
- z2 = ((target_distance - distance2) / m.in_units("Mpccm / h")) + z2
- distance2 = self.cosmology.comoving_radial_distance(z2, z)
- iteration += 1
- if (iteration > max_Iterations):
- mylog.error("calculate_deltaz_max: Warning - max iterations " +
- "exceeded for z = %f (delta z = %f)." %
- (z, np.abs(z2 - z)))
- break
- # Use this calculation or the absolute minimum specified by the user.
- output['dz_min'] = max(np.abs(z2 - z), deltaz_min)
+ zf = self._deltaz_forward(output['redshift'],
+ target_distance)
+ output['dz_min'] = max(zf, deltaz_min)
def _deltaz_forward(self, z, target_distance):
r"""Calculate deltaz corresponding to moving a comoving distance
@@ -357,10 +316,13 @@
d_Tolerance = 1e-4
max_Iterations = 100
- # Calculate delta z that corresponds to the length of the
- # box at a given redshift.
z1 = z
- z2 = z1 - 0.1 # just an initial guess
+ # Use Hubble's law for initial guess
+ target_distance = self.cosmology.quan(target_distance.to("Mpccm / h"))
+ v = self.cosmology.hubble_parameter(z) * target_distance
+ v = min(v, 0.9 * c)
+ dz = np.sqrt((1. + v/c) / (1. - v/c)) - 1.
+ z2 = z1 - dz
distance1 = self.cosmology.quan(0.0, "Mpccm / h")
distance2 = self.cosmology.comoving_radial_distance(z2, z)
iteration = 1
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
--- a/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
+++ b/yt/analysis_modules/cosmological_observation/light_ray/light_ray.py
@@ -79,21 +79,23 @@
max_box_fraction : optional, float
In terms of the size of the domain, the maximum length a light
ray segment can be in order to span the redshift interval from
- one dataset to another.
+ one dataset to another. If using a zoom-in simulation, this
+ parameter can be set to the length of the high resolution
+ region so as to limit ray segments to that size. If the
+ high resolution region is not cubical, the smallest side
+ should be used.
Default: 1.0 (the size of the box)
deltaz_min : optional, float
Specifies the minimum :math:`\Delta z` between consecutive
datasets in the returned list. Do not use for simple rays.
Default: 0.0.
minimum_coherent_box_fraction : optional, float
- Used with use_minimum_datasets set to False, this parameter
- specifies the fraction of the total box size to be traversed
- before rerandomizing the projection axis and center. This
- was invented to allow light rays with thin slices to sample
- coherent large scale structure, but in practice does not work
- so well. Try setting this parameter to 1 and see what happens.
- Do not use for simple rays.
- Default: 0.0.
+ Use to specify the minimum length of a ray, in terms of the
+ size of the domain, before the trajectory is re-randomized.
+ Set to 0 to have ray trajectory randomized for every dataset.
+ Set to np.inf (infinity) to use a single trajectory for the
+ entire ray.
+ Default: 0.
time_data : optional, bool
Whether or not to include time outputs when gathering
datasets for time series. Do not use for simple rays.
@@ -123,6 +125,11 @@
time_data=True, redshift_data=True,
find_outputs=False, load_kwargs=None):
+ if near_redshift is not None and far_redshift is not None and \
+ near_redshift >= far_redshift:
+ raise RuntimeError(
+ "near_redshift must be less than far_redshift.")
+
self.near_redshift = near_redshift
self.far_redshift = far_redshift
self.use_minimum_datasets = use_minimum_datasets
@@ -154,8 +161,7 @@
self.cosmology = Cosmology(
hubble_constant=self.ds.hubble_constant,
omega_matter=self.ds.omega_matter,
- omega_lambda=self.ds.omega_lambda,
- unit_registry=self.ds.unit_registry)
+ omega_lambda=self.ds.omega_lambda)
else:
redshift = 0.
self.light_ray_solution.append({"filename": self.parameter_filename,
@@ -169,20 +175,23 @@
CosmologySplice.__init__(self, self.parameter_filename, simulation_type,
find_outputs=find_outputs)
self.light_ray_solution = \
- self.create_cosmology_splice(self.near_redshift, self.far_redshift,
- minimal=self.use_minimum_datasets,
- max_box_fraction=max_box_fraction,
- deltaz_min=self.deltaz_min,
- time_data=time_data,
- redshift_data=redshift_data)
+ self.create_cosmology_splice(
+ self.near_redshift, self.far_redshift,
+ minimal=self.use_minimum_datasets,
+ max_box_fraction=max_box_fraction,
+ deltaz_min=self.deltaz_min,
+ time_data=time_data,
+ redshift_data=redshift_data)
def _calculate_light_ray_solution(self, seed=None,
+ left_edge=None, right_edge=None,
+ min_level=None, periodic=True,
start_position=None, end_position=None,
trajectory=None, filename=None):
"Create list of datasets to be added together to make the light ray."
# Calculate dataset sizes, and get random dataset axes and centers.
- np.random.seed(seed)
+ my_random = np.random.RandomState(seed)
# If using only one dataset, set start and stop manually.
if start_position is not None:
@@ -192,9 +201,9 @@
if not ((end_position is None) ^ (trajectory is None)):
raise RuntimeError("LightRay Error: must specify either end_position " + \
"or trajectory, but not both.")
- self.light_ray_solution[0]['start'] = np.asarray(start_position)
+ self.light_ray_solution[0]['start'] = start_position
if end_position is not None:
- self.light_ray_solution[0]['end'] = np.asarray(end_position)
+ self.light_ray_solution[0]['end'] = end_position
else:
# assume trajectory given as r, theta, phi
if len(trajectory) != 3:
@@ -228,29 +237,40 @@
# Get dataset axis and center.
# If using box coherence, only get start point and vector if
- # enough of the box has been used,
- # or if box_fraction_used will be greater than 1 after this slice.
- if (q == 0) or (self.minimum_coherent_box_fraction == 0) or \
- (box_fraction_used >
- self.minimum_coherent_box_fraction) or \
- (box_fraction_used +
- self.light_ray_solution[q]['traversal_box_fraction'] > 1.0):
- # Random start point
- self.light_ray_solution[q]['start'] = np.random.random(3)
- theta = np.pi * np.random.random()
- phi = 2 * np.pi * np.random.random()
- box_fraction_used = 0.0
+ # enough of the box has been used.
+ if (q == 0) or (box_fraction_used >=
+ self.minimum_coherent_box_fraction):
+ if periodic:
+ self.light_ray_solution[q]['start'] = left_edge + \
+ (right_edge - left_edge) * my_random.random_sample(3)
+ theta = np.pi * my_random.random_sample()
+ phi = 2 * np.pi * my_random.random_sample()
+ box_fraction_used = 0.0
+ else:
+ ds = load(self.light_ray_solution[q]["filename"])
+ ray_length = \
+ ds.quan(self.light_ray_solution[q]['traversal_box_fraction'],
+ "unitary")
+ self.light_ray_solution[q]['start'], \
+ self.light_ray_solution[q]['end'] = \
+ non_periodic_ray(ds, left_edge, right_edge, ray_length,
+ my_random=my_random, min_level=min_level)
+ del ds
else:
- # Use end point of previous segment and same theta and phi.
+ # Use end point of previous segment, adjusted for periodicity,
+ # and the same trajectory.
self.light_ray_solution[q]['start'] = \
- self.light_ray_solution[q-1]['end'][:]
+ periodic_adjust(self.light_ray_solution[q-1]['end'][:],
+ left=left_edge, right=right_edge)
- self.light_ray_solution[q]['end'] = \
- self.light_ray_solution[q]['start'] + \
- self.light_ray_solution[q]['traversal_box_fraction'] * \
- np.array([np.cos(phi) * np.sin(theta),
- np.sin(phi) * np.sin(theta),
- np.cos(theta)])
+ if "end" not in self.light_ray_solution[q]:
+ self.light_ray_solution[q]['end'] = \
+ self.light_ray_solution[q]['start'] + \
+ self.light_ray_solution[q]['traversal_box_fraction'] * \
+ self.simulation.box_size * \
+ np.array([np.cos(phi) * np.sin(theta),
+ np.sin(phi) * np.sin(theta),
+ np.cos(theta)])
box_fraction_used += \
self.light_ray_solution[q]['traversal_box_fraction']
@@ -261,15 +281,18 @@
'far_redshift':self.far_redshift,
'near_redshift':self.near_redshift})
- def make_light_ray(self, seed=None,
+ def make_light_ray(self, seed=None, periodic=True,
+ left_edge=None, right_edge=None, min_level=None,
start_position=None, end_position=None,
trajectory=None,
fields=None, setup_function=None,
solution_filename=None, data_filename=None,
- get_los_velocity=None, use_peculiar_velocity=True,
- redshift=None, njobs=-1):
+ get_los_velocity=None, use_peculiar_velocity=True,
+ redshift=None, field_parameters=None, njobs=-1):
"""
- make_light_ray(seed=None, start_position=None, end_position=None,
+ make_light_ray(seed=None, periodic=True,
+ left_edge=None, right_edge=None, min_level=None,
+ start_position=None, end_position=None,
trajectory=None, fields=None, setup_function=None,
solution_filename=None, data_filename=None,
use_peculiar_velocity=True, redshift=None,
@@ -285,6 +308,29 @@
seed : optional, int
Seed for the random number generator.
Default: None.
+ periodic : optional, bool
+ If True, ray trajectories will make use of periodic
+ boundaries. If False, ray trajectories will not be
+ periodic.
+ Default : True.
+ left_edge : optional, iterable of floats or YTArray
+ The left corner of the region in which rays are to be
+ generated. If None, the left edge will be that of the
+ domain. If specified without units, it is assumed to
+ be in code units.
+ Default: None.
+ right_edge : optional, iterable of floats or YTArray
+ The right corner of the region in which rays are to be
+ generated. If None, the right edge will be that of the
+ domain. If specified without units, it is assumed to
+ be in code units.
+ Default: None.
+ min_level : optional, int
+ The minimum refinement level of the spatial region in which
+ the ray passes. This can be used with zoom-in simulations
+ where the high resolution region does not keep a constant
+ geometry.
+ Default: None.
start_position : optional, iterable of floats or YTArray.
Used only if creating a light ray from a single dataset.
The coordinates of the starting position of the ray.
@@ -363,30 +409,56 @@
... use_peculiar_velocity=True)
"""
+ if self.simulation_type is None:
+ domain = self.ds
+ else:
+ domain = self.simulation
- if start_position is not None and hasattr(start_position, 'units'):
- start_position = start_position.to('unitary')
- elif start_position is not None :
- start_position = self.ds.arr(
- start_position, 'code_length').to('unitary')
+ assumed_units = "code_length"
+ if left_edge is None:
+ left_edge = domain.domain_left_edge
+ elif not hasattr(left_edge, 'units'):
+ left_edge = domain.arr(left_edge, assumed_units)
+ left_edge.convert_to_units('unitary')
- if end_position is not None and hasattr(end_position, 'units'):
- end_position = end_position.to('unitary')
- elif end_position is not None :
- end_position = self.ds.arr(
- end_position, 'code_length').to('unitary')
+ if right_edge is None:
+ right_edge = domain.domain_right_edge
+ elif not hasattr(right_edge, 'units'):
+ right_edge = domain.arr(right_edge, assumed_units)
+ right_edge.convert_to_units('unitary')
+
+ if start_position is not None:
+ if hasattr(start_position, 'units'):
+ start_position = start_position
+ else:
+ start_position = self.ds.arr(start_position, assumed_units)
+ start_position.convert_to_units('unitary')
+
+ if end_position is not None:
+ if hasattr(end_position, 'units'):
+ end_position = end_position
+ else:
+ end_position = self.ds.arr(end_position, assumed_units)
+ end_position.convert_to_units('unitary')
if get_los_velocity is not None:
use_peculiar_velocity = get_los_velocity
- mylog.warn("'get_los_velocity' kwarg is deprecated. Use 'use_peculiar_velocity' instead.")
+ mylog.warn("'get_los_velocity' kwarg is deprecated. " + \
+ "Use 'use_peculiar_velocity' instead.")
# Calculate solution.
self._calculate_light_ray_solution(seed=seed,
+ left_edge=left_edge,
+ right_edge=right_edge,
+ min_level=min_level, periodic=periodic,
start_position=start_position,
end_position=end_position,
trajectory=trajectory,
filename=solution_filename)
+ if field_parameters is None:
+ field_parameters = {}
+
# Initialize data structures.
self._data = {}
# temperature field is automatically added to fields
@@ -427,19 +499,11 @@
if setup_function is not None:
setup_function(ds)
- if start_position is not None:
- my_segment["start"] = ds.arr(my_segment["start"], "unitary")
- my_segment["end"] = ds.arr(my_segment["end"], "unitary")
- else:
- my_segment["start"] = ds.domain_width * my_segment["start"] + \
- ds.domain_left_edge
- my_segment["end"] = ds.domain_width * my_segment["end"] + \
- ds.domain_left_edge
-
if not ds.cosmological_simulation:
next_redshift = my_segment["redshift"]
elif self.near_redshift == self.far_redshift:
- if isinstance(my_segment["traversal_box_fraction"], YTArray):
+ if isinstance(my_segment["traversal_box_fraction"], YTArray) and \
+ not my_segment["traversal_box_fraction"].units.is_dimensionless:
segment_length = \
my_segment["traversal_box_fraction"].in_units("Mpccm / h")
else:
@@ -453,18 +517,18 @@
else:
next_redshift = my_segment['next']['redshift']
+ # Make sure start, end, left, right
+ # are using the dataset's unit system.
+ my_start = ds.arr(my_segment['start'])
+ my_end = ds.arr(my_segment['end'])
+ my_left = ds.arr(left_edge)
+ my_right = ds.arr(right_edge)
mylog.info("Getting segment at z = %s: %s to %s." %
- (my_segment['redshift'], my_segment['start'],
- my_segment['end']))
-
- # Convert segment units from unitary to code length for sub_ray
- my_segment['start'] = my_segment['start'].to('code_length')
- my_segment['end'] = my_segment['end'].to('code_length')
+ (my_segment['redshift'], my_start, my_end))
# Break periodic ray into non-periodic segments.
- sub_segments = periodic_ray(my_segment['start'], my_segment['end'],
- left=ds.domain_left_edge,
- right=ds.domain_right_edge)
+ sub_segments = periodic_ray(my_start, my_end,
+ left=my_left, right=my_right)
# Prepare data structure for subsegment.
sub_data = {}
@@ -477,6 +541,8 @@
mylog.info("Getting subsegment: %s to %s." %
(list(sub_segment[0]), list(sub_segment[1])))
sub_ray = ds.ray(sub_segment[0], sub_segment[1])
+ for key, val in field_parameters.items():
+ sub_ray.set_field_parameter(key, val)
asort = np.argsort(sub_ray["t"])
sub_data['dl'].extend(sub_ray['dts'][asort] *
vector_length(sub_ray.start_point,
@@ -515,7 +581,7 @@
# sight) and the velocity vectors: a dot b = ab cos(theta)
sub_vel_mag = sub_ray['velocity_magnitude']
- cos_theta = np.dot(line_of_sight, sub_vel) / sub_vel_mag
+ cos_theta = line_of_sight.dot(sub_vel) / sub_vel_mag
# Protect against stituations where velocity mag is exactly
# zero, in which case zero / zero = NaN.
cos_theta = np.nan_to_num(cos_theta)
@@ -535,8 +601,7 @@
# Get redshift for each lixel. Assume linear relation between l
# and z.
sub_data['dredshift'] = (my_segment['redshift'] - next_redshift) * \
- (sub_data['dl'] / vector_length(my_segment['start'],
- my_segment['end']).in_cgs())
+ (sub_data['dl'] / vector_length(my_start, my_end).in_cgs())
sub_data['redshift'] = my_segment['redshift'] - \
sub_data['dredshift'].cumsum() + sub_data['dredshift']
@@ -672,6 +737,22 @@
return np.sqrt(np.power((end - start), 2).sum())
+def periodic_adjust(p, left=None, right=None):
+ """
+ Return the point p adjusted for periodic boundaries.
+
+ """
+ if isinstance(p, YTArray):
+ p.convert_to_units("unitary")
+ if left is None:
+ left = np.zeros_like(p)
+ if right is None:
+ right = np.ones_like(p)
+
+ w = right - left
+ p -= left
+ return np.mod(p, w)
+
def periodic_distance(coord1, coord2):
"""
periodic_distance(coord1, coord2)
@@ -713,7 +794,7 @@
dim = right - left
vector = end - start
- wall = np.zeros(start.shape)
+ wall = np.zeros_like(start)
close = np.zeros(start.shape, dtype=object)
left_bound = vector < 0
@@ -733,7 +814,6 @@
this_end = end.copy()
t = 0.0
tolerance = 1e-6
-
while t < 1.0 - tolerance:
hit_left = (this_start <= left) & (vector < 0)
if (hit_left).any():
@@ -751,8 +831,44 @@
now = this_start + vector * dt
close_enough = np.abs(now - nearest) / np.abs(vector.max()) < 1e-10
now[close_enough] = nearest[close_enough]
- segments.append([np.copy(this_start), np.copy(now)])
+ segments.append([this_start.copy(), now.copy()])
this_start = now.copy()
t += dt
return segments
+
+def non_periodic_ray(ds, left_edge, right_edge, ray_length, max_iter=5000,
+ min_level=None, my_random=None):
+
+ max_length = vector_length(left_edge, right_edge)
+ if ray_length > max_length:
+ raise RuntimeError(
+ ("The maximum segment length in the region %s to %s is %s, " +
+ "but the ray length requested is %s. Decrease ray length.") %
+ (left_edge, right_edge, max_length, ray_length))
+
+ if my_random is None:
+ my_random = np.random.RandomState()
+ i = 0
+ while True:
+ start = my_random.random_sample(3) * \
+ (right_edge - left_edge) + left_edge
+ theta = np.pi * my_random.random_sample()
+ phi = 2 * np.pi * my_random.random_sample()
+ end = start + ray_length * \
+ np.array([np.cos(phi) * np.sin(theta),
+ np.sin(phi) * np.sin(theta),
+ np.cos(theta)])
+ i += 1
+ test_ray = ds.ray(start, end)
+ if (end >= left_edge).all() and (end <= right_edge).all() and \
+ (min_level is None or min_level <= 0 or
+ (test_ray["grid_level"] >= min_level).all()):
+ mylog.info("Found ray after %d attempts." % i)
+ del test_ray
+ return start, end
+ del test_ray
+ if i > max_iter:
+ raise RuntimeError(
+ ("Failed to create segment in %d attempts. " +
+ "Decreasing ray length is recommended") % i)
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/analysis_modules/cosmological_observation/light_ray/tests/test_light_ray.py
--- a/yt/analysis_modules/cosmological_observation/light_ray/tests/test_light_ray.py
+++ b/yt/analysis_modules/cosmological_observation/light_ray/tests/test_light_ray.py
@@ -10,6 +10,8 @@
# The full license is in the file COPYING.txt, distributed with this software.
#-----------------------------------------------------------------------------
+import numpy as np
+
from yt.testing import \
requires_file
from yt.analysis_modules.cosmological_observation.api import LightRay
@@ -41,6 +43,48 @@
os.chdir(curdir)
shutil.rmtree(tmpdir)
+ at requires_file(COSMO_PLUS)
+def test_light_ray_cosmo_nested():
+ """
+ This test generates a cosmological light ray confing the ray to a subvolume
+ """
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ left = np.ones(3) * 0.25
+ right = np.ones(3) * 0.75
+
+ lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.03)
+
+ lr.make_light_ray(seed=1234567, left_edge=left, right_edge=right,
+ fields=['temperature', 'density', 'H_number_density'],
+ data_filename='lightray.h5')
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
+
+ at requires_file(COSMO_PLUS)
+def test_light_ray_cosmo_nonperiodic():
+ """
+ This test generates a cosmological light ray using non-periodic segments
+ """
+ # Set up in a temp dir
+ tmpdir = tempfile.mkdtemp()
+ curdir = os.getcwd()
+ os.chdir(tmpdir)
+
+ lr = LightRay(COSMO_PLUS, 'Enzo', 0.0, 0.03)
+
+ lr.make_light_ray(seed=1234567, periodic=False,
+ fields=['temperature', 'density', 'H_number_density'],
+ data_filename='lightray.h5')
+
+ # clean up
+ os.chdir(curdir)
+ shutil.rmtree(tmpdir)
@requires_file(COSMO_PLUS_SINGLE)
def test_light_ray_non_cosmo():
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -197,10 +197,20 @@
def __init__(self, start_point, end_point, ds=None,
field_parameters=None, data_source=None):
super(YTRay, self).__init__(ds, field_parameters, data_source)
- self.start_point = self.ds.arr(start_point,
- 'code_length', dtype='float64')
- self.end_point = self.ds.arr(end_point,
- 'code_length', dtype='float64')
+ if isinstance(start_point, YTArray):
+ self.start_point = \
+ self.ds.arr(start_point).to("code_length")
+ else:
+ self.start_point = \
+ self.ds.arr(start_point, 'code_length',
+ dtype='float64')
+ if isinstance(end_point, YTArray):
+ self.end_point = \
+ self.ds.arr(end_point).to("code_length")
+ else:
+ self.end_point = \
+ self.ds.arr(end_point, 'code_length',
+ dtype='float64')
self.vec = self.end_point - self.start_point
self._set_center(self.start_point)
self.set_field_parameter('center', self.start_point)
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/data_objects/static_output.py
--- a/yt/data_objects/static_output.py
+++ b/yt/data_objects/static_output.py
@@ -868,8 +868,7 @@
self.cosmology = \
Cosmology(hubble_constant=self.hubble_constant,
omega_matter=self.omega_matter,
- omega_lambda=self.omega_lambda,
- unit_registry=self.unit_registry)
+ omega_lambda=self.omega_lambda)
self.critical_density = \
self.cosmology.critical_density(self.current_redshift)
self.scale_factor = 1.0 / (1.0 + self.current_redshift)
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/frontends/enzo/simulation_handling.py
--- a/yt/frontends/enzo/simulation_handling.py
+++ b/yt/frontends/enzo/simulation_handling.py
@@ -110,6 +110,8 @@
self.domain_right_edge = self.domain_right_edge * self.length_unit
self.unit_registry.modify("code_time", self.time_unit)
self.unit_registry.modify("code_length", self.length_unit)
+ self.unit_registry.add("unitary", float(self.box_size.in_base()),
+ self.length_unit.units.dimensions)
def get_time_series(self, time_data=True, redshift_data=True,
initial_time=None, final_time=None,
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/frontends/gadget/simulation_handling.py
--- a/yt/frontends/gadget/simulation_handling.py
+++ b/yt/frontends/gadget/simulation_handling.py
@@ -102,6 +102,8 @@
self.box_size = self.box_size * self.length_unit
self.domain_left_edge = self.domain_left_edge * self.length_unit
self.domain_right_edge = self.domain_right_edge * self.length_unit
+ self.unit_registry.add("unitary", float(self.box_size.in_base()),
+ self.length_unit.units.dimensions)
else:
# Read time from file for non-cosmological sim
self.time_unit = self.quan(
diff -r 9ccda6b503c68b4944523f4dc8afb5dbbfb94b1e -r 4bc8ace93fda81e800a2e90c052c431c2d48a01b yt/utilities/cosmology.py
--- a/yt/utilities/cosmology.py
+++ b/yt/utilities/cosmology.py
@@ -33,7 +33,14 @@
For an explanation of the various cosmological measures, see, for example
Hogg (1999, http://xxx.lanl.gov/abs/astro-ph/9905116).
-
+
+ WARNING: Cosmological distance calculations return values that are either
+ in the comoving or proper frame, depending on the specific quantity. For
+ simplicity, the proper and comoving frames are set equal to each other
+ within the cosmology calculator. This means that for some distance value,
+ x, x.to("Mpc") and x.to("Mpccm") will be the same. The user should take
+ care to understand which reference frame is correct for the given calculation.
+
Parameters
----------
hubble_constant : float
@@ -58,7 +65,7 @@
>>> from yt.utilities.cosmology import Cosmology
>>> co = Cosmology()
>>> print(co.hubble_time(0.0).in_units("Gyr"))
-
+
"""
def __init__(self, hubble_constant = 0.71,
omega_matter = 0.27,
@@ -66,9 +73,9 @@
omega_curvature = 0.0,
unit_registry = None,
unit_system = "cgs"):
- self.omega_matter = omega_matter
- self.omega_lambda = omega_lambda
- self.omega_curvature = omega_curvature
+ self.omega_matter = float(omega_matter)
+ self.omega_lambda = float(omega_lambda)
+ self.omega_curvature = float(omega_curvature)
if unit_registry is None:
unit_registry = UnitRegistry()
unit_registry.modify("h", hubble_constant)
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.
More information about the yt-svn
mailing list