[yt-svn] commit/yt: xarthisius: Merged in jzuhone/yt (pull request #2149)
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Jan 17 08:12:40 PST 2017
1 new commit in yt:
https://bitbucket.org/yt_analysis/yt/commits/e03273b808d6/
Changeset: e03273b808d6
Branch: yt
User: xarthisius
Date: 2017-01-17 16:12:12+00:00
Summary: Merged in jzuhone/yt (pull request #2149)
Athena++ frontend
Affected #: 21 files
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -98,13 +98,13 @@
ds = load("./A11QR1/s11Qzm1h2_a1.0000.art")
-.. _loading_athena_data:
+.. _loading-athena-data:
Athena Data
-----------
-Athena 4.x VTK data is *mostly* supported and cared for by John
-ZuHone. Both uniform grid and SMR datasets are supported.
+Athena 4.x VTK data is supported and cared for by John ZuHone. Both uniform grid
+and SMR datasets are supported.
.. note:
yt also recognizes Fargo3D data written to VTK files as
@@ -138,10 +138,14 @@
which will pick up all of the files in the different ``id*`` directories for
the entire dataset.
-yt works in cgs ("Gaussian") units by default, but Athena data is not
-normally stored in these units. If you would like to convert data to
-cgs units, you may supply conversions for length, time, and mass to ``load`` using
-the ``units_override`` functionality:
+The default unit system in yt is cgs ("Gaussian") units, but Athena data is not
+normally stored in these units, so the code unit system is the default unit
+system for Athena data. This means that answers to field queries from data
+objects and plots of data will be expressed in code units. Note that the default
+conversions from these units will still be in terms of cgs units, e.g. 1
+``code_length`` equals 1 cm, and so on. If you would like to provided different
+conversions, you may supply conversions for length, time, and mass to ``load``
+using the ``units_override`` functionality:
.. code-block:: python
@@ -153,14 +157,16 @@
ds = yt.load("id0/cluster_merger.0250.vtk", units_override=units_override)
-This means that the yt fields, e.g. ``("gas","density")``, ``("gas","x-velocity")``,
-``("gas","magnetic_field_x")``, will be in cgs units, but the Athena fields, e.g.,
-``("athena","density")``, ``("athena","velocity_x")``, ``("athena","cell_centered_B_x")``, will be
-in code units.
+This means that the yt fields, e.g. ``("gas","density")``,
+``("gas","velocity_x")``, ``("gas","magnetic_field_x")``, will be in cgs units
+(or whatever unit system was specified), but the Athena fields, e.g.,
+``("athena","density")``, ``("athena","velocity_x")``,
+``("athena","cell_centered_B_x")``, will be in code units.
-Some 3D Athena outputs may have large grids (especially parallel datasets subsequently joined with
-the `join_vtk` script), and may benefit from being subdivided into "virtual grids". For this purpose,
-one can pass in the `nprocs` parameter:
+Some 3D Athena outputs may have large grids (especially parallel datasets
+subsequently joined with the ``join_vtk`` script), and may benefit from being
+subdivided into "virtual grids". For this purpose, one can pass in the
+``nprocs`` parameter:
.. code-block:: python
@@ -168,43 +174,100 @@
ds = yt.load("sloshing.0000.vtk", nprocs=8)
-which will subdivide each original grid into `nprocs` grids.
+which will subdivide each original grid into ``nprocs`` grids.
.. note::
Virtual grids are only supported (and really only necessary) for 3D data.
-Alternative values for the following simulation parameters may be specified using a ``parameters``
-dict, accepting the following keys:
+Alternative values for the following simulation parameters may be specified
+using a ``parameters`` dict, accepting the following keys:
* ``Gamma``: ratio of specific heats, Type: Float
-* ``geometry``: Geometry type, currently accepts ``"cartesian"`` or ``"cylindrical"``
-* ``periodicity``: Is the domain periodic? Type: Tuple of boolean values corresponding to each dimension
+* ``geometry``: Geometry type, currently accepts ``"cartesian"`` or
+ ``"cylindrical"``
+* ``periodicity``: Is the domain periodic? Type: Tuple of boolean values
+ corresponding to each dimension
.. code-block:: python
import yt
- parameters = {"gamma":4./3., "geometry":"cylindrical", "periodicity":(False,False,False)}
+ parameters = {"gamma":4./3., "geometry":"cylindrical",
+ "periodicity":(False,False,False)}
ds = yt.load("relativistic_jet_0000.vtk", parameters=parameters)
.. rubric:: Caveats
-* yt primarily works with primitive variables. If the Athena
- dataset contains conservative variables, the yt primitive fields will be generated from the
+* yt primarily works with primitive variables. If the Athena dataset contains
+ conservative variables, the yt primitive fields will be generated from the
conserved variables on disk.
-* Special relativistic datasets may be loaded, but are not fully supported. In particular, the relationships between
- quantities such as pressure and thermal energy will be incorrect, as it is currently assumed that their relationship
- is that of an ideal a :math:`\gamma`-law equation of state.
+* Special relativistic datasets may be loaded, but at this time not all of
+ their fields are fully supported. In particular, the relationships between
+ quantities such as pressure and thermal energy will be incorrect, as it is
+ currently assumed that their relationship is that of an ideal a
+ :math:`\gamma`-law equation of state. This will be rectified in a future
+ release.
* Domains may be visualized assuming periodicity.
* Particle list data is currently unsupported.
.. note::
The old behavior of supplying unit conversions using a ``parameters``
- dict supplied to ``load`` for Athena datasets is still supported, but is being deprecated in
- favor of ``units_override``, which provides the same functionality.
+ dict supplied to ``load`` for Athena datasets is still supported, but is
+ being deprecated in favor of ``units_override``, which provides the same
+ functionality.
+
+.. _loading-athena-pp-data:
+
+Athena++ Data
+-------------
+
+Athena++ HDF5 data is supported and cared for by John ZuHone. Uniform-grid, SMR,
+and AMR datasets in cartesian coordinates are fully supported. Support for
+curvilinear coordinates and logarithmic cell sizes exists, but is preliminary.
+For the latter type of dataset, the data will be loaded in as a semi-structured
+mesh dataset. See :ref:`loading-semi-structured-mesh-data` for more details on
+how this works in yt.
+
+The default unit system in yt is cgs ("Gaussian") units, but Athena++ data is
+not normally stored in these units, so the code unit system is the default unit
+system for Athena++ data. This means that answers to field queries from data
+objects and plots of data will be expressed in code units. Note that the default
+conversions from these units will still be in terms of cgs units, e.g. 1
+``code_length`` equals 1 cm, and so on. If you would like to provided different
+conversions, you may supply conversions for length, time, and mass to ``load``
+using the ``units_override`` functionality:
+
+.. code-block:: python
+
+ import yt
+
+ units_override = {"length_unit":(1.0,"Mpc"),
+ "time_unit"(1.0,"Myr"),
+ "mass_unit":(1.0e14,"Msun")}
+
+ ds = yt.load("AM06/AM06.out1.00400.athdf", units_override=units_override)
+
+This means that the yt fields, e.g. ``("gas","density")``,
+``("gas","velocity_x")``, ``("gas","magnetic_field_x")``, will be in cgs units
+(or whatever unit system was specified), but the Athena fields, e.g.,
+``("athena_pp","density")``, ``("athena_pp","vel1")``, ``("athena_pp","Bcc1")``,
+will be in code units.
+
+.. rubric:: Caveats
+
+* yt primarily works with primitive variables. If the Athena++ dataset contains
+ conservative variables, the yt primitive fields will be generated from the
+ conserved variables on disk.
+* Special relativistic datasets may be loaded, but at this time not all of their
+ fields are fully supported. In particular, the relationships between
+ quantities such as pressure and thermal energy will be incorrect, as it is
+ currently assumed that their relationship is that of an ideal
+ :math:`\gamma`-law equation of state. This will be rectified in a future
+ release.
+* Domains may be visualized assuming periodicity.
.. _loading-orion-data:
@@ -1182,6 +1245,8 @@
* Particles may be difficult to integrate.
* Data must already reside in memory.
+.. _loading-semi-structured-mesh-data:
+
Semi-Structured Grid Data
-------------------------
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 tests/tests.yaml
--- a/tests/tests.yaml
+++ b/tests/tests.yaml
@@ -2,9 +2,12 @@
local_artio_001:
- yt/frontends/artio/tests/test_outputs.py
- local_athena_002:
+ local_athena_003:
- yt/frontends/athena
+ local_athena_pp_000:
+ - yt/frontends/athena_pp
+
local_chombo_002:
- yt/frontends/chombo/tests/test_outputs.py
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/fields/magnetic_field.py
--- a/yt/fields/magnetic_field.py
+++ b/yt/fields/magnetic_field.py
@@ -36,15 +36,17 @@
def setup_magnetic_field_fields(registry, ftype = "gas", slice_info = None):
unit_system = registry.ds.unit_system
- if (ftype,"magnetic_field_x") not in registry:
+ axis_names = registry.ds.coordinates.axis_order
+
+ if (ftype,"magnetic_field_%s" % axis_names[0]) not in registry:
return
- u = registry[ftype,"magnetic_field_x"].units
+ u = registry[ftype,"magnetic_field_%s" % axis_names[0]].units
def _magnetic_field_strength(field,data):
- B2 = (data[ftype,"magnetic_field_x"]**2 +
- data[ftype,"magnetic_field_y"]**2 +
- data[ftype,"magnetic_field_z"]**2)
+ B2 = (data[ftype,"magnetic_field_%s" % axis_names[0]]**2 +
+ data[ftype,"magnetic_field_%s" % axis_names[1]]**2 +
+ data[ftype,"magnetic_field_%s" % axis_names[2]]**2)
return np.sqrt(B2)
registry.add_field((ftype,"magnetic_field_strength"), sampling_type="cell",
function=_magnetic_field_strength,
@@ -69,37 +71,63 @@
function=_magnetic_pressure,
units=unit_system["pressure"])
- def _magnetic_field_poloidal(field,data):
- normal = data.get_field_parameter("normal")
- d = data[ftype,'magnetic_field_x']
- Bfields = data.ds.arr(
- [data[ftype,'magnetic_field_x'],
- data[ftype,'magnetic_field_y'],
- data[ftype,'magnetic_field_z']],
- d.units)
+ if registry.ds.geometry == "cartesian":
+ def _magnetic_field_poloidal(field,data):
+ normal = data.get_field_parameter("normal")
+ d = data[ftype,'magnetic_field_x']
+ Bfields = data.ds.arr(
+ [data[ftype,'magnetic_field_x'],
+ data[ftype,'magnetic_field_y'],
+ data[ftype,'magnetic_field_z']],
+ d.units)
- theta = data["index", 'spherical_theta']
- phi = data["index", 'spherical_phi']
+ theta = data["index", 'spherical_theta']
+ phi = data["index", 'spherical_phi']
- return get_sph_theta_component(Bfields, theta, phi, normal)
+ return get_sph_theta_component(Bfields, theta, phi, normal)
+
+ def _magnetic_field_toroidal(field,data):
+ normal = data.get_field_parameter("normal")
+ d = data[ftype,'magnetic_field_x']
+ Bfields = data.ds.arr(
+ [data[ftype,'magnetic_field_x'],
+ data[ftype,'magnetic_field_y'],
+ data[ftype,'magnetic_field_z']],
+ d.units)
+
+ phi = data["index", 'spherical_phi']
+ return get_sph_phi_component(Bfields, phi, normal)
+
+ elif registry.ds.geometry == "cylindrical":
+ def _magnetic_field_poloidal(field, data):
+ r = data["index", "r"]
+ z = data["index", "z"]
+ d = np.sqrt(r*r+z*z)
+ return (data[ftype, "magnetic_field_r"]*(r/d) +
+ data[ftype, "magnetic_field_z"]*(z/d))
+
+ def _magnetic_field_toroidal(field, data):
+ return data[ftype,"magnetic_field_theta"]
+
+ elif registry.ds.geometry == "spherical":
+ def _magnetic_field_poloidal(field, data):
+ return data[ftype,"magnetic_field_theta"]
+
+ def _magnetic_field_toroidal(field, data):
+ return data[ftype,"magnetic_field_phi"]
+
+ else:
+
+ # Unidentified geometry--set to None
+
+ _magnetic_field_toroidal = None
+ _magnetic_field_poloidal = None
registry.add_field((ftype, "magnetic_field_poloidal"), sampling_type="cell",
function=_magnetic_field_poloidal,
units=u, validators=[ValidateParameter("normal")])
- def _magnetic_field_toroidal(field,data):
- normal = data.get_field_parameter("normal")
- d = data[ftype,'magnetic_field_x']
- Bfields = data.ds.arr(
- [data[ftype,'magnetic_field_x'],
- data[ftype,'magnetic_field_y'],
- data[ftype,'magnetic_field_z']],
- d.units)
-
- phi = data["index", 'spherical_phi']
- return get_sph_phi_component(Bfields, phi, normal)
-
- registry.add_field((ftype, "magnetic_field_toroidal"), sampling_type="cell",
+ registry.add_field((ftype, "magnetic_field_toroidal"), sampling_type="cell",
function=_magnetic_field_toroidal,
units=u, validators=[ValidateParameter("normal")])
@@ -160,7 +188,7 @@
def _mag_field(field, data):
return convert(data[fd])
return _mag_field
- for ax, fd in zip("xyz", ds_fields):
+ for ax, fd in zip(registry.ds.coordinates.axis_order, ds_fields):
registry.add_field((ftype,"magnetic_field_%s" % ax), sampling_type="cell",
function=mag_field(fd),
units=unit_system[to_units.dimensions])
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/api.py
--- a/yt/frontends/api.py
+++ b/yt/frontends/api.py
@@ -20,6 +20,7 @@
'art',
'artio',
'athena',
+ 'athena_pp',
'boxlib',
'chombo',
'eagle',
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena/fields.py
--- a/yt/frontends/athena/fields.py
+++ b/yt/frontends/athena/fields.py
@@ -60,29 +60,17 @@
self.add_field(("gas","velocity_%s" % comp), sampling_type="cell",
function=velocity_field(comp), units = unit_system["velocity"])
# Add pressure, energy, and temperature fields
- def ekin1(data):
- return 0.5*(data["athena","momentum_x"]**2 +
- data["athena","momentum_y"]**2 +
- data["athena","momentum_z"]**2)/data["athena","density"]
- def ekin2(data):
- return 0.5*(data["athena","velocity_x"]**2 +
- data["athena","velocity_y"]**2 +
- data["athena","velocity_z"]**2)*data["athena","density"]
- def emag(data):
- return 0.5*(data["cell_centered_B_x"]**2 +
- data["cell_centered_B_y"]**2 +
- data["cell_centered_B_z"]**2)
def eint_from_etot(data):
eint = data["athena","total_energy"]
- eint -= ekin1(data)
+ eint -= data["gas","kinetic_energy"]
if ("athena","cell_centered_B_x") in self.field_list:
- eint -= emag(data)
+ eint -= data["gas","magnetic_energy"]
return eint
def etot_from_pres(data):
etot = data["athena","pressure"]/(data.ds.gamma-1.)
- etot += ekin2(data)
+ etot += data["gas", "kinetic_energy"]
if ("athena","cell_centered_B_x") in self.field_list:
- etot += emag(data)
+ etot += data["gas", "magnetic_energy"]
return etot
if ("athena","pressure") in self.field_list:
self.add_output_field(("athena","pressure"), sampling_type="cell",
@@ -103,10 +91,6 @@
elif ("athena","total_energy") in self.field_list:
self.add_output_field(("athena","total_energy"), sampling_type="cell",
units=pres_units)
- def _pressure(field, data):
- return eint_from_etot(data)*(data.ds.gamma-1.0)
- self.add_field(("gas","pressure"), sampling_type="cell", function=_pressure,
- units=unit_system["pressure"])
def _thermal_energy(field, data):
return eint_from_etot(data)/data["athena","density"]
self.add_field(("gas","thermal_energy"), sampling_type="cell",
@@ -117,7 +101,7 @@
self.add_field(("gas","total_energy"), sampling_type="cell",
function=_total_energy,
units=unit_system["specific_energy"])
-
+ # Add temperature field
def _temperature(field, data):
if data.has_field_parameter("mu"):
mu = data.get_field_parameter("mu")
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena_pp/api.py
--- /dev/null
+++ b/yt/frontends/athena_pp/api.py
@@ -0,0 +1,26 @@
+"""
+API for yt.frontends.athena++
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+from .data_structures import \
+ AthenaPPGrid, \
+ AthenaPPHierarchy, \
+ AthenaPPDataset
+
+from .fields import \
+ AthenaPPFieldInfo
+
+from .io import \
+ IOHandlerAthenaPP
+
+from . import tests
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena_pp/data_structures.py
--- /dev/null
+++ b/yt/frontends/athena_pp/data_structures.py
@@ -0,0 +1,360 @@
+"""
+Data structures for Athena.
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+import os
+import weakref
+
+from yt.funcs import \
+ mylog, get_pbar, \
+ ensure_tuple
+from yt.data_objects.grid_patch import \
+ AMRGridPatch
+from yt.geometry.grid_geometry_handler import \
+ GridIndex
+from yt.data_objects.static_output import \
+ Dataset
+from yt.geometry.geometry_handler import \
+ YTDataChunk
+from yt.utilities.file_handler import \
+ HDF5FileHandler
+from yt.geometry.unstructured_mesh_handler import \
+ UnstructuredIndex
+from yt.data_objects.unstructured_mesh import \
+ SemiStructuredMesh
+from itertools import chain, product
+from .fields import AthenaPPFieldInfo
+
+geom_map = {"cartesian": "cartesian",
+ "cylindrical": "cylindrical",
+ "spherical_polar": "spherical",
+ "minkowski": "cartesian",
+ "tilted": "cartesian",
+ "sinusoidal": "cartesian",
+ "schwarzschild": "spherical",
+ "kerr-schild": "spherical"}
+
+_cis = np.fromiter(chain.from_iterable(product([0,1], [0,1], [0,1])),
+ dtype=np.int64, count = 8*3)
+_cis.shape = (8, 3)
+
+class AthenaPPLogarithmicMesh(SemiStructuredMesh):
+ _index_offset = 0
+
+ def __init__(self, mesh_id, filename, connectivity_indices,
+ connectivity_coords, index, blocks, dims):
+ super(AthenaPPLogarithmicMesh, self).__init__(mesh_id, filename,
+ connectivity_indices,
+ connectivity_coords, index)
+ self.mesh_blocks = blocks
+ self.mesh_dims = dims
+
+class AthenaPPLogarithmicIndex(UnstructuredIndex):
+ def __init__(self, ds, dataset_type = 'athena_pp'):
+ self._handle = ds._handle
+ super(AthenaPPLogarithmicIndex, self).__init__(ds, dataset_type)
+ self.index_filename = self.dataset.filename
+ self.directory = os.path.dirname(self.dataset.filename)
+ self.dataset_type = dataset_type
+
+ def _initialize_mesh(self):
+ mylog.debug("Setting up meshes.")
+ num_blocks = self._handle.attrs["NumMeshBlocks"]
+ log_loc = self._handle['LogicalLocations']
+ levels = self._handle["Levels"]
+ x1f = self._handle["x1f"]
+ x2f = self._handle["x2f"]
+ x3f = self._handle["x3f"]
+ nbx, nby, nbz = tuple(np.max(log_loc, axis=0)+1)
+ nlevel = self._handle.attrs["MaxLevel"]+1
+
+ nb = np.array([nbx, nby, nbz], dtype='int')
+ self.mesh_factors = np.ones(3, dtype='int')*((nb > 1).astype("int")+1)
+
+ block_grid = -np.ones((nbx,nby,nbz,nlevel), dtype=np.int)
+ block_grid[log_loc[:,0],log_loc[:,1],log_loc[:,2],levels[:]] = np.arange(num_blocks)
+
+ block_list = np.arange(num_blocks, dtype='int64')
+ bc = []
+ for i in range(num_blocks):
+ if block_list[i] >= 0:
+ ii, jj, kk = log_loc[i]
+ neigh = block_grid[ii:ii+2,jj:jj+2,kk:kk+2,levels[i]]
+ if np.all(neigh > -1):
+ loc_ids = neigh.transpose().flatten()
+ bc.append(loc_ids)
+ block_list[loc_ids] = -1
+ else:
+ bc.append(np.array(i))
+ block_list[i] = -1
+
+ num_meshes = len(bc)
+
+ self.meshes = []
+ pbar = get_pbar("Constructing meshes", num_meshes)
+ for i in range(num_meshes):
+ ob = bc[i][0]
+ x = x1f[ob,:]
+ y = x2f[ob,:]
+ z = x3f[ob,:]
+ if nbx > 1:
+ x = np.concatenate([x, x1f[bc[i][1],1:]])
+ if nby > 1:
+ y = np.concatenate([y, x2f[bc[i][2],1:]])
+ if nbz > 1:
+ z = np.concatenate([z, x3f[bc[i][4],1:]])
+ nxm = x.size
+ nym = y.size
+ nzm = z.size
+ coords = np.zeros((nxm, nym, nzm, 3), dtype="float64", order="C")
+ coords[:,:,:,0] = x[:,None,None]
+ coords[:,:,:,1] = y[None,:,None]
+ coords[:,:,:,2] = z[None,None,:]
+ coords.shape = (nxm * nym * nzm, 3)
+ cycle = np.rollaxis(np.indices((nxm-1,nym-1,nzm-1)), 0, 4)
+ cycle.shape = ((nxm-1)*(nym-1)*(nzm-1), 3)
+ off = _cis + cycle[:, np.newaxis]
+ connectivity = ((off[:,:,0] * nym) + off[:,:,1]) * nzm + off[:,:,2]
+ mesh = AthenaPPLogarithmicMesh(i, self.index_filename, connectivity,
+ coords, self, bc[i],
+ np.array([nxm-1, nym-1, nzm-1]))
+ self.meshes.append(mesh)
+ pbar.update(i)
+ pbar.finish()
+ mylog.debug("Done setting up meshes.")
+
+ def _detect_output_fields(self):
+ self.field_list = [("athena_pp", k) for k in self.ds._field_map]
+
+ def _chunk_io(self, dobj, cache = True, local_only = False):
+ gobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+ for subset in gobjs:
+ yield YTDataChunk(dobj, "io", [subset],
+ self._count_selection(dobj, [subset]),
+ cache = cache)
+
+class AthenaPPGrid(AMRGridPatch):
+ _id_offset = 0
+
+ def __init__(self, id, index, level):
+ AMRGridPatch.__init__(self, id, filename = index.index_filename,
+ index = index)
+ self.Parent = None
+ self.Children = []
+ self.Level = level
+
+ def _setup_dx(self):
+ # So first we figure out what the index is. We don't assume
+ # that dx=dy=dz , at least here. We probably do elsewhere.
+ id = self.id - self._id_offset
+ LE, RE = self.index.grid_left_edge[id,:], \
+ self.index.grid_right_edge[id,:]
+ self.dds = self.ds.arr((RE-LE)/self.ActiveDimensions, "code_length")
+ if self.ds.dimensionality < 2: self.dds[1] = 1.0
+ if self.ds.dimensionality < 3: self.dds[2] = 1.0
+ self.field_data['dx'], self.field_data['dy'], self.field_data['dz'] = self.dds
+
+ def __repr__(self):
+ return "AthenaPPGrid_%04i (%s)" % (self.id, self.ActiveDimensions)
+
+class AthenaPPHierarchy(GridIndex):
+
+ grid = AthenaPPGrid
+ _dataset_type='athena_pp'
+ _data_file = None
+
+ def __init__(self, ds, dataset_type='athena_pp'):
+ self.dataset = weakref.proxy(ds)
+ self.directory = os.path.dirname(self.dataset.filename)
+ self.dataset_type = dataset_type
+ # for now, the index file is the dataset!
+ self.index_filename = self.dataset.filename
+ self._handle = ds._handle
+ GridIndex.__init__(self, ds, dataset_type)
+
+ def _detect_output_fields(self):
+ self.field_list = [("athena_pp", k) for k in self.dataset._field_map]
+
+ def _count_grids(self):
+ self.num_grids = self._handle.attrs["NumMeshBlocks"]
+
+ def _parse_index(self):
+ num_grids = self._handle.attrs["NumMeshBlocks"]
+
+ self.grid_left_edge = np.zeros((num_grids, 3), dtype='float64')
+ self.grid_right_edge = np.zeros((num_grids, 3), dtype='float64')
+ self.grid_dimensions = np.zeros((num_grids, 3), dtype='int32')
+
+ for i in range(num_grids):
+ x = self._handle["x1f"][i,:]
+ y = self._handle["x2f"][i,:]
+ z = self._handle["x3f"][i,:]
+ self.grid_left_edge[i] = np.array([x[0], y[0], z[0]], dtype='float64')
+ self.grid_right_edge[i] = np.array([x[-1], y[-1], z[-1]], dtype='float64')
+ self.grid_dimensions[i] = self._handle.attrs["MeshBlockSize"]
+ levels = self._handle["Levels"][:]
+
+ self.grid_left_edge = self.ds.arr(self.grid_left_edge, "code_length")
+ self.grid_right_edge = self.ds.arr(self.grid_right_edge, "code_length")
+
+ self.grids = np.empty(self.num_grids, dtype='object')
+ for i in range(num_grids):
+ self.grids[i] = self.grid(i, self, levels[i])
+
+ if self.dataset.dimensionality <= 2:
+ self.grid_right_edge[:,2] = self.dataset.domain_right_edge[2]
+ if self.dataset.dimensionality == 1:
+ self.grid_right_edge[:,1:] = self.dataset.domain_right_edge[1:]
+ self.grid_particle_count = np.zeros([self.num_grids, 1], dtype='int64')
+
+ def _populate_grid_objects(self):
+ for g in self.grids:
+ g._prepare_grid()
+ g._setup_dx()
+ self.max_level = self._handle.attrs["MaxLevel"]
+
+ def _chunk_io(self, dobj, cache = True, local_only = False):
+ gobjs = getattr(dobj._current_chunk, "objs", dobj._chunk_info)
+ for subset in gobjs:
+ yield YTDataChunk(dobj, "io", [subset],
+ self._count_selection(dobj, [subset]),
+ cache = cache)
+
+class AthenaPPDataset(Dataset):
+ _field_info_class = AthenaPPFieldInfo
+ _dataset_type = "athena_pp"
+
+ def __init__(self, filename, dataset_type='athena_pp',
+ storage_filename=None, parameters=None,
+ units_override=None, unit_system="code"):
+ self.fluid_types += ("athena_pp",)
+ if parameters is None:
+ parameters = {}
+ self.specified_parameters = parameters
+ if units_override is None:
+ units_override = {}
+ self._handle = HDF5FileHandler(filename)
+ xrat = self._handle.attrs["RootGridX1"][2]
+ yrat = self._handle.attrs["RootGridX2"][2]
+ zrat = self._handle.attrs["RootGridX3"][2]
+ if xrat != 1.0 or yrat != 1.0 or zrat != 1.0:
+ self._index_class = AthenaPPLogarithmicIndex
+ self.logarithmic = True
+ else:
+ self._index_class = AthenaPPHierarchy
+ self.logarithmic = False
+ Dataset.__init__(self, filename, dataset_type, units_override=units_override,
+ unit_system=unit_system)
+ self.filename = filename
+ if storage_filename is None:
+ storage_filename = '%s.yt' % filename.split('/')[-1]
+ self.storage_filename = storage_filename
+ self.backup_filename = self.filename[:-4] + "_backup.gdf"
+
+ def _set_code_unit_attributes(self):
+ """
+ Generates the conversion to various physical _units based on the
+ parameter file
+ """
+ if "length_unit" not in self.units_override:
+ self.no_cgs_equiv_length = True
+ for unit, cgs in [("length", "cm"), ("time", "s"), ("mass", "g"),
+ ("temperature", "K")]:
+ # We set these to cgs for now, but they may have been overriden
+ if getattr(self, unit+'_unit', None) is not None:
+ continue
+ mylog.warning("Assuming 1.0 = 1.0 %s", cgs)
+ setattr(self, "%s_unit" % unit, self.quan(1.0, cgs))
+
+ def set_code_units(self):
+ super(AthenaPPDataset, self).set_code_units()
+ mag_unit = getattr(self, "magnetic_unit", None)
+ if mag_unit is None:
+ self.magnetic_unit = np.sqrt(4*np.pi * self.mass_unit /
+ (self.time_unit**2 * self.length_unit))
+ self.magnetic_unit.convert_to_units("gauss")
+ temp_unit = getattr(self, "temperature_unit", None)
+ if temp_unit is None:
+ self.temperature_unit = self.quan(1.0, "K")
+ self.velocity_unit = self.length_unit / self.time_unit
+
+ def _parse_parameter_file(self):
+
+ xmin, xmax = self._handle.attrs["RootGridX1"][:2]
+ ymin, ymax = self._handle.attrs["RootGridX2"][:2]
+ zmin, zmax = self._handle.attrs["RootGridX3"][:2]
+
+ self.domain_left_edge = np.array([xmin, ymin, zmin], dtype='float64')
+ self.domain_right_edge = np.array([xmax, ymax, zmax], dtype='float64')
+
+ self.geometry = geom_map[self._handle.attrs["Coordinates"].decode('utf-8')]
+ self.domain_width = self.domain_right_edge-self.domain_left_edge
+ self.domain_dimensions = self._handle.attrs["RootGridSize"]
+
+ self._field_map = {}
+ k = 0
+ for (i, dname), num_var in zip(enumerate(self._handle.attrs["DatasetNames"]),
+ self._handle.attrs["NumVariables"]):
+ for j in range(num_var):
+ fname = self._handle.attrs["VariableNames"][k].decode("ascii","ignore")
+ self._field_map[fname] = (dname.decode("ascii","ignore"), j)
+ k += 1
+
+ self.refine_by = 2
+ dimensionality = 3
+ if self.domain_dimensions[2] == 1:
+ dimensionality = 2
+ if self.domain_dimensions[1] == 1:
+ dimensionality = 1
+ self.dimensionality = dimensionality
+ self.current_time = self._handle.attrs["Time"]
+ self.unique_identifier = self.parameter_filename.__hash__()
+ self.cosmological_simulation = False
+ self.num_ghost_zones = 0
+ self.field_ordering = 'fortran'
+ self.boundary_conditions = [1]*6
+ if 'periodicity' in self.specified_parameters:
+ self.periodicity = ensure_tuple(self.specified_parameters['periodicity'])
+ else:
+ self.periodicity = (True,True,True,)
+ if 'gamma' in self.specified_parameters:
+ self.gamma = float(self.specified_parameters['gamma'])
+ else:
+ self.gamma = 5./3.
+
+ self.current_redshift = self.omega_lambda = self.omega_matter = \
+ self.hubble_constant = self.cosmological_simulation = 0.0
+ self.parameters['Time'] = self.current_time # Hardcode time conversion for now.
+ self.parameters["HydroMethod"] = 0 # Hardcode for now until field staggering is supported.
+ if "gamma" in self.specified_parameters:
+ self.parameters["Gamma"] = self.specified_parameters["gamma"]
+ else:
+ self.parameters["Gamma"] = 5./3.
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ if args[0].endswith('athdf'):
+ return True
+ except:
+ pass
+ return False
+
+ @property
+ def _skip_cache(self):
+ return True
+
+ def __repr__(self):
+ return self.basename.rsplit(".", 1)[0]
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena_pp/definitions.py
--- /dev/null
+++ b/yt/frontends/athena_pp/definitions.py
@@ -0,0 +1,14 @@
+"""
+Various definitions for various other modules and routines
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena_pp/fields.py
--- /dev/null
+++ b/yt/frontends/athena_pp/fields.py
@@ -0,0 +1,91 @@
+"""
+Athena++-specific fields
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+from yt.fields.field_info_container import \
+ FieldInfoContainer
+from yt.utilities.physical_constants import \
+ kboltz, mh
+
+b_units = "code_magnetic"
+pres_units = "code_mass/(code_length*code_time**2)"
+rho_units = "code_mass / code_length**3"
+vel_units = "code_length / code_time"
+
+def velocity_field(j):
+ def _velocity(field, data):
+ return data["athena_pp", "mom%d" % j]/data["athena_pp","dens"]
+ return _velocity
+
+class AthenaPPFieldInfo(FieldInfoContainer):
+ known_other_fields = (
+ ("rho", (rho_units, ["density"], None)),
+ ("dens", (rho_units, ["density"], None)),
+ ("Bcc1", (b_units, [], None)),
+ ("Bcc2", (b_units, [], None)),
+ ("Bcc3", (b_units, [], None)),
+ )
+
+ def setup_fluid_fields(self):
+ from yt.fields.magnetic_field import \
+ setup_magnetic_field_aliases
+ unit_system = self.ds.unit_system
+ # Add velocity fields
+ vel_prefix = "velocity"
+ for i, comp in enumerate(self.ds.coordinates.axis_order):
+ vel_field = ("athena_pp", "vel%d" % (i+1))
+ mom_field = ("athena_pp", "mom%d" % (i+1))
+ if vel_field in self.field_list:
+ self.add_output_field(vel_field, sampling_type="cell", units="code_length/code_time")
+ self.alias(("gas","%s_%s" % (vel_prefix, comp)), vel_field,
+ units=unit_system["velocity"])
+ elif mom_field in self.field_list:
+ self.add_output_field(mom_field, sampling_type="cell",
+ units="code_mass/code_time/code_length**2")
+ self.add_field(("gas","%s_%s" % (vel_prefix, comp)), sampling_type="cell",
+ function=velocity_field(i+1), units=unit_system["velocity"])
+ # Figure out thermal energy field
+ if ("athena_pp","press") in self.field_list:
+ self.add_output_field(("athena_pp","press"), sampling_type="cell",
+ units=pres_units)
+ self.alias(("gas","pressure"),("athena_pp","press"),
+ units=unit_system["pressure"])
+ def _thermal_energy(field, data):
+ return data["athena_pp","press"] / \
+ (data.ds.gamma-1.)/data["athena_pp","rho"]
+ self.add_field(("gas","thermal_energy"), sampling_type="cell",
+ function=_thermal_energy,
+ units=unit_system["specific_energy"])
+ elif ("athena_pp","Etot") in self.field_list:
+ self.add_output_field(("athena_pp","Etot"), sampling_type="cell",
+ units=pres_units)
+ def _thermal_energy(field, data):
+ eint = data["athena_pp", "Etot"] - data["gas","kinetic_energy"]
+ if ("athena_pp", "B1") in self.field_list:
+ eint -= data["gas","magnetic_energy"]
+ return eint/data["athena_pp","dens"]
+ self.add_field(("gas","thermal_energy"), sampling_type="cell",
+ function=_thermal_energy,
+ units=unit_system["specific_energy"])
+ # Add temperature field
+ def _temperature(field, data):
+ if data.has_field_parameter("mu"):
+ mu = data.get_field_parameter("mu")
+ else:
+ mu = 0.6
+ return (data["gas","pressure"]/data["gas","density"])*mu*mh/kboltz
+ self.add_field(("gas", "temperature"), sampling_type="cell", function=_temperature,
+ units=unit_system["temperature"])
+
+ setup_magnetic_field_aliases(self, "athena_pp", ["Bcc%d" % ax for ax in (1,2,3)])
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena_pp/io.py
--- /dev/null
+++ b/yt/frontends/athena_pp/io.py
@@ -0,0 +1,100 @@
+"""
+Athena++-specific IO functions
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+from itertools import groupby
+
+from yt.utilities.io_handler import \
+ BaseIOHandler
+from yt.utilities.logger import ytLogger as mylog
+
+# http://stackoverflow.com/questions/2361945/detecting-consecutive-integers-in-a-list
+def grid_sequences(grids):
+ g_iter = sorted(grids, key = lambda g: g.id)
+ for k, g in groupby(enumerate(g_iter), lambda i_x1:i_x1[0]-i_x1[1].id):
+ seq = list(v[1] for v in g)
+ yield seq
+
+ii = [0, 1, 0, 1, 0, 1, 0, 1]
+jj = [0, 0, 1, 1, 0, 0, 1, 1]
+kk = [0, 0, 0, 0, 1, 1, 1, 1]
+
+class IOHandlerAthenaPP(BaseIOHandler):
+ _particle_reader = False
+ _dataset_type = "athena_pp"
+
+ def __init__(self, ds):
+ super(IOHandlerAthenaPP, self).__init__(ds)
+ self._handle = ds._handle
+
+ def _read_particles(self, fields_to_read, type, args, grid_list,
+ count_list, conv_factors):
+ pass
+
+ def _read_fluid_selection(self, chunks, selector, fields, size):
+ chunks = list(chunks)
+ if any((ftype != "athena_pp" for ftype, fname in fields)):
+ raise NotImplementedError
+ f = self._handle
+ rv = {}
+ for field in fields:
+ # Always use *native* 64-bit float.
+ rv[field] = np.empty(size, dtype="=f8")
+ ng = sum(len(c.objs) for c in chunks)
+ mylog.debug("Reading %s cells of %s fields in %s blocks",
+ size, [f2 for f1, f2 in fields], ng)
+ for field in fields:
+ ftype, fname = field
+ dname, fdi = self.ds._field_map[fname]
+ ds = f["/%s" % dname]
+ ind = 0
+ for chunk in chunks:
+ if self.ds.logarithmic:
+ for mesh in chunk.objs:
+ nx, ny, nz = mesh.mesh_dims // self.ds.index.mesh_factors
+ data = np.empty(mesh.mesh_dims, dtype="=f8")
+ for n, id in enumerate(mesh.mesh_blocks):
+ data[ii[n]*nx:(ii[n]+1)*nx,jj[n]*ny:(jj[n]+1)*ny,kk[n]*nz:(kk[n]+1)*nz] = \
+ ds[fdi,id,:,:,:].transpose()
+ ind += mesh.select(selector, data, rv[field], ind) # caches
+ else:
+ for gs in grid_sequences(chunk.objs):
+ start = gs[0].id - gs[0]._id_offset
+ end = gs[-1].id - gs[-1]._id_offset + 1
+ data = ds[fdi,start:end,:,:,:].transpose()
+ for i, g in enumerate(gs):
+ ind += g.select(selector, data[...,i], rv[field], ind)
+ return rv
+
+ def _read_chunk_data(self, chunk, fields):
+ if self.ds.logarithmic:
+ pass
+ f = self._handle
+ rv = {}
+ for g in chunk.objs:
+ rv[g.id] = {}
+ if len(fields) == 0:
+ return rv
+ for field in fields:
+ ftype, fname = field
+ dname, fdi = self.ds._field_map[fname]
+ ds = f["/%s" % dname]
+ for gs in grid_sequences(chunk.objs):
+ start = gs[0].id - gs[0]._id_offset
+ end = gs[-1].id - gs[-1]._id_offset + 1
+ data = ds[fdi,start:end,:,:,:].transpose()
+ for i, g in enumerate(gs):
+ rv[g.id][field] = np.asarray(data[...,i], "=f8")
+ return rv
\ No newline at end of file
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/frontends/athena_pp/tests/test_outputs.py
--- /dev/null
+++ b/yt/frontends/athena_pp/tests/test_outputs.py
@@ -0,0 +1,79 @@
+"""
+Athena++ frontend tests
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2013, yt Development Team.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+import numpy as np
+from yt.testing import \
+ assert_equal, \
+ requires_file, \
+ units_override_check, \
+ assert_allclose
+from yt.utilities.answer_testing.framework import \
+ requires_ds, \
+ small_patch_amr, \
+ data_dir_load, \
+ GenericArrayTest
+from yt.frontends.athena_pp.api import AthenaPPDataset
+from yt.convenience import load
+
+_fields_disk = ("density", "velocity_r")
+
+disk = "KeplerianDisk/disk.out1.00000.athdf"
+ at requires_ds(disk)
+def test_disk():
+ ds = data_dir_load(disk)
+ assert_equal(str(ds), "disk.out1.00000")
+ dd = ds.all_data()
+ vol = (ds.domain_right_edge[0]**3-ds.domain_left_edge[0]**3)/3.0
+ vol *= np.cos(ds.domain_left_edge[1])-np.cos(ds.domain_right_edge[1])
+ vol *= ds.domain_right_edge[2].v-ds.domain_left_edge[2].v
+ assert_allclose(dd.quantities.total_quantity("cell_volume"), vol)
+ for field in _fields_disk:
+ def field_func(name):
+ return dd[field]
+ yield GenericArrayTest(ds, field_func, args=[field])
+
+_fields_AM06 = ("temperature", "density", "velocity_magnitude", "magnetic_field_x")
+
+AM06 = "AM06/AM06.out1.00400.athdf"
+ at requires_ds(AM06)
+def test_AM06():
+ ds = data_dir_load(AM06)
+ assert_equal(str(ds), "AM06.out1.00400")
+ for test in small_patch_amr(ds, _fields_AM06):
+ test_AM06.__name__ = test.description
+ yield test
+
+uo_AM06 = {
+ 'length_unit': (1.0, 'kpc'),
+ 'mass_unit': (1.0, 'Msun'),
+ 'time_unit': (1.0, 'Myr'),
+}
+
+ at requires_file(AM06)
+def test_AM06_override():
+ # verify that overriding units causes derived unit values to be updated.
+ # see issue #1259
+ ds = load(AM06, units_override=uo_AM06)
+ assert_equal(float(ds.magnetic_unit.in_units('gauss')),
+ 9.01735778342523e-08)
+
+ at requires_file(AM06)
+def test_units_override():
+ for test in units_override_check(AM06):
+ yield test
+
+ at requires_file(AM06)
+def test_AthenaPPDataset():
+ assert isinstance(data_dir_load(AM06), AthenaPPDataset)
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/units/tests/test_unit_systems.py
--- a/yt/units/tests/test_unit_systems.py
+++ b/yt/units/tests/test_unit_systems.py
@@ -108,7 +108,7 @@
ad = ds.all_data()
dens = ad['density']
magx = ad['magx']
- magnetic_field_x = ad['magnetic_field_x']
+ magnetic_field_x = ad['magnetic_field_r']
if us == 'cgs':
assert str(dens.units) == 'g/cm**3'
diff -r b7b71a7c2e00c7a43b7e736a9111fe7d4051b7e8 -r e03273b808d64c8e09a1914149bcea38e16e4831 yt/visualization/plot_window.py
--- a/yt/visualization/plot_window.py
+++ b/yt/visualization/plot_window.py
@@ -1314,6 +1314,10 @@
if field_parameters is None:
field_parameters = {}
+ if ds.geometry == "spherical" or ds.geometry == "cylindrical":
+ mylog.info("Setting origin='native' for %s geometry." % ds.geometry)
+ origin = 'native'
+
if isinstance(ds, YTSpatialPlotDataset):
slc = ds.all_data()
slc.axis = axis
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