[yt-svn] commit/yt-3.0: 12 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Jun 11 06:32:57 PDT 2013
12 new commits in yt-3.0:
https://bitbucket.org/yt_analysis/yt-3.0/commits/7b55808f11cf/
Changeset: 7b55808f11cf
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-04 23:28:54
Summary: Minor fix to get __repr__ for data_collection to work properly.
Affected #: 1 file
diff -r 5e745eda811bc188d30587bafa6ca9fa01878835 -r 7b55808f11cfd8cdee936fdd95bc547b70be4482 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -894,7 +894,7 @@
Child cells are not returned.
"""
_type_name = "data_collection"
- _con_args = ("obj_list",)
+ _con_args = ("_obj_list",)
def __init__(self, center, obj_list, pf = None, field_parameters = None):
YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
self._obj_ids = np.array([o.id - o._id_offset for o in obj_list],
https://bitbucket.org/yt_analysis/yt-3.0/commits/cf67430914e7/
Changeset: cf67430914e7
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-05 00:16:27
Summary: This change enables deferral of fields that cannot be generated.
The next step will be refactoring the "can I generate" code into its own
routine and attempting to speed it up. After that we can address the idea of
chunks that are IO oriented preloading spatial chunks.
Affected #: 1 file
diff -r 7b55808f11cfd8cdee936fdd95bc547b70be4482 -r cf67430914e72b7fb9892980f856f029122983da yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -477,8 +477,23 @@
if fields is None: return
fields = self._determine_fields(fields)
# Now we collect all our fields
- fields_to_get = [f for f in fields if f not in self.field_data]
- if len(fields_to_get) == 0:
+ # Here is where we need to perform a validation step, so that if we
+ # have a field requested that we actually *can't* yet get, we put it
+ # off until the end. This prevents double-reading fields that will
+ # need to be used in spatial fields later on.
+ fields_to_get = []
+ # This will be pre-populated with spatial fields
+ fields_to_generate = []
+ for field in self._determine_fields(fields):
+ if field in self.field_data: continue
+ finfo = self.pf._get_field_info(*field)
+ try:
+ finfo.check_available(self)
+ except NeedsGridType:
+ fields_to_generate.append(field)
+ continue
+ fields_to_get.append(field)
+ if len(fields_to_get) == 0 and fields_to_generate == 0:
return
elif self._locked == True:
raise GenerationInProgress(fields)
@@ -502,12 +517,18 @@
read_particles, gen_particles = self.hierarchy._read_particle_fields(
particles, self, self._current_chunk)
self.field_data.update(read_particles)
- fields_to_generate = gen_fluids + gen_particles
+ fields_to_generate += gen_fluids + gen_particles
self._generate_fields(fields_to_generate)
def _generate_fields(self, fields_to_generate):
index = 0
with self._field_lock():
+ # At this point, we assume that any fields that are necessary to
+ # *generate* a field are in fact already available to us. Note
+ # that we do not make any assumption about whether or not the
+ # fields have a spatial requirement. This will be checked inside
+ # _generate_field, at which point additional dependencies may
+ # actually be noted.
while any(f not in self.field_data for f in fields_to_generate):
field = fields_to_generate[index % len(fields_to_generate)]
index += 1
https://bitbucket.org/yt_analysis/yt-3.0/commits/10c24e37517f/
Changeset: 10c24e37517f
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-05 20:55:17
Summary: Moving the OWLSStaticOutput below Gadget, subclassing, and setting up fields.
This brings OWLS mostly in feature parity with Gadget, although I noticed that
I had to make some odd changes to the cm units to make it work the same way.
Affected #: 2 files
diff -r 443f359e25980fe2574652a4c873b5ca9a702169 -r 10c24e37517f0d437f700786f50e5068bc1ad78d yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -163,71 +163,6 @@
for subset in oobjs:
yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
-class OWLSStaticOutput(StaticOutput):
- _hierarchy_class = ParticleGeometryHandler
- _domain_class = ParticleDomainFile
- _fieldinfo_fallback = OWLSFieldInfo
- _fieldinfo_known = KnownOWLSFields
-
- def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
- self._root_dimensions = root_dimensions
- # Set up the template for domain files
- self.storage_filename = None
- super(OWLSStaticOutput, self).__init__(filename, data_style)
-
- def __repr__(self):
- return os.path.basename(self.parameter_filename).split(".")[0]
-
- def _set_units(self):
- self.units = {}
- self.time_units = {}
- self.conversion_factors = {}
- DW = self.domain_right_edge - self.domain_left_edge
- self.units["unitary"] = 1.0 / DW.max()
-
- def _parse_parameter_file(self):
- handle = h5py.File(self.parameter_filename)
- hvals = {}
- hvals.update(handle["/Header"].attrs)
-
- self.dimensionality = 3
- self.refine_by = 2
- self.parameters["HydroMethod"] = "sph"
- self.unique_identifier = \
- int(os.stat(self.parameter_filename)[stat.ST_CTIME])
- # Set standard values
- self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
- self.domain_left_edge = np.zeros(3, "float64")
- self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
- self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
- self.cosmological_simulation = 1
- self.periodicity = (True, True, True)
- self.current_redshift = hvals["Redshift"]
- self.omega_lambda = hvals["OmegaLambda"]
- self.omega_matter = hvals["Omega0"]
- self.hubble_constant = hvals["HubbleParam"]
- self.parameters = hvals
-
- prefix = self.parameter_filename.split(".", 1)[0]
- suffix = self.parameter_filename.rsplit(".", 1)[-1]
- self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
- self.domain_count = hvals["NumFilesPerSnapshot"]
-
- handle.close()
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- try:
- fileh = h5py.File(args[0],'r')
- if "Constants" in fileh["/"].keys() and \
- "Header" in fileh["/"].keys():
- fileh.close()
- return True
- fileh.close()
- except:
- pass
- return False
-
class GadgetBinaryDomainFile(ParticleDomainFile):
def __init__(self, pf, io, domain_filename, domain_id):
with open(domain_filename, "rb") as f:
@@ -404,6 +339,74 @@
# We do not allow load() of these files.
return False
+class OWLSStaticOutput(GadgetStaticOutput):
+ _hierarchy_class = ParticleGeometryHandler
+ _domain_class = ParticleDomainFile
+ _fieldinfo_fallback = OWLSFieldInfo # For now we have separate from Gadget
+ _fieldinfo_known = KnownOWLSFields
+ _header_spec = None # Override so that there's no confusion
+
+ def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
+ self._root_dimensions = root_dimensions
+ # Set up the template for domain files
+ self.storage_filename = None
+ super(OWLSStaticOutput, self).__init__(filename, data_style,
+ root_dimensions,
+ unit_base = None)
+
+ def __repr__(self):
+ return os.path.basename(self.parameter_filename).split(".")[0]
+
+ def _parse_parameter_file(self):
+ handle = h5py.File(self.parameter_filename)
+ hvals = {}
+ hvals.update((str(k), v) for k, v in handle["/Header"].attrs.items())
+
+ self.dimensionality = 3
+ self.refine_by = 2
+ self.parameters["HydroMethod"] = "sph"
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[stat.ST_CTIME])
+
+ # Set standard values
+ self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
+ self.domain_left_edge = np.zeros(3, "float64")
+ self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
+ self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
+ self.cosmological_simulation = 1
+ self.periodicity = (True, True, True)
+ self.current_redshift = hvals["Redshift"]
+ self.omega_lambda = hvals["OmegaLambda"]
+ self.omega_matter = hvals["Omega0"]
+ self.hubble_constant = hvals["HubbleParam"]
+ self.parameters = hvals
+
+ prefix = self.parameter_filename.split(".", 1)[0]
+ suffix = self.parameter_filename.rsplit(".", 1)[-1]
+ self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
+ self.domain_count = hvals["NumFilesPerSnapshot"]
+
+ # To avoid having to open files twice
+ self._unit_base = {}
+ self._unit_base.update((str(k), v) for k, v in handle["/Units"].attrs.items())
+ # Comoving cm is given in the Units
+ self._unit_base['cmcm'] = 1.0 / self._unit_base["UnitLength_in_cm"]
+
+ handle.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = h5py.File(args[0],'r')
+ if "Constants" in fileh["/"].keys() and \
+ "Header" in fileh["/"].keys():
+ fileh.close()
+ return True
+ fileh.close()
+ except:
+ pass
+ return False
+
class TipsyDomainFile(ParticleDomainFile):
def _calculate_offsets(self, field_list):
diff -r 443f359e25980fe2574652a4c873b5ca9a702169 -r 10c24e37517f0d437f700786f50e5068bc1ad78d yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -247,3 +247,60 @@
GadgetFieldInfo.add_field(("all", oname + ax), function=func,
particle_type = True)
+# OWLS
+# ====
+
+# I am optimistic that some day we will be able to get rid of much of this, and
+# make OWLS a subclass of Gadget fields.
+
+_owls_ptypes = ("PartType0", "PartType1", "PartType2", "PartType3",
+ "PartType4")
+
+for fname in ["Coordinates", "Velocities", "ParticleIDs",
+ # Note: Mass, not Masses
+ "Mass"]:
+ func = _field_concat(fname)
+ OWLSFieldInfo.add_field(("all", fname), function=func,
+ particle_type = True)
+
+def _owls_particle_fields(ptype):
+ def _Mass(field, data):
+ pind = _owls_ptypes.index(ptype)
+ if data.pf["MassTable"][pind] == 0.0:
+ raise RuntimeError
+ mass = np.ones(data[ptype, "ParticleIDs"].shape[0], dtype="float64")
+ # Note that this is an alias, which is why we need to apply conversion
+ # here. Otherwise we'd have an asymmetry.
+ mass *= data.pf["MassTable"][pind]
+ return mass
+ OWLSFieldInfo.add_field((ptype, "Mass"), function=_Mass,
+ convert_function = _get_conv("mass"),
+ particle_type = True)
+
+for ptype in _owls_ptypes:
+ # Note that this adds a "Known" Mass field and a "Derived" Mass field.
+ # This way the "Known" will get used, and if it's not there, it will use
+ # the derived.
+ KnownOWLSFields.add_field((ptype, "Mass"), function=NullFunc,
+ particle_type = True,
+ convert_function=_get_conv("mass"),
+ units = r"\mathrm{g}")
+ _owls_particle_fields(ptype)
+ KnownOWLSFields.add_field((ptype, "Velocities"), function=NullFunc,
+ particle_type = True,
+ convert_function=_get_conv("velocity"),
+ units = r"\mathrm{cm}/\mathrm{s}")
+ _particle_functions(ptype, "Coordinates", "Mass", OWLSFieldInfo)
+ KnownOWLSFields.add_field((ptype, "Coordinates"), function=NullFunc,
+ particle_type = True)
+_particle_functions("all", "Coordinates", "Mass", OWLSFieldInfo)
+
+# Now we have to manually apply the splits for "all", since we don't want to
+# use the splits defined above.
+
+for iname, oname in [("Coordinates", "particle_position_"),
+ ("Velocities", "particle_velocity_")]:
+ for axi, ax in enumerate("xyz"):
+ func = _field_concat_slice(iname, axi)
+ OWLSFieldInfo.add_field(("all", oname + ax), function=func,
+ particle_type = True)
https://bitbucket.org/yt_analysis/yt-3.0/commits/2dd46f58a1f4/
Changeset: 2dd46f58a1f4
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-05 21:42:54
Summary: Fixing "MaxValue" and "MaxPosition" naming for fields that are tuples.
Affected #: 1 file
diff -r 10c24e37517f0d437f700786f50e5068bc1ad78d -r 2dd46f58a1f4e08f360680cb83954ce207873328 yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -78,7 +78,7 @@
source.quantities["MaxLocation"](field)
mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f",
max_val, mx, my, mz)
- self.pf.parameters["Max%sValue" % (field)] = max_val
- self.pf.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
+ self.pf.parameters["Max%sValue" % (field,)] = max_val
+ self.pf.parameters["Max%sPos" % (field,)] = "%s" % ((mx,my,mz),)
return max_val, np.array((mx,my,mz), dtype='float64')
https://bitbucket.org/yt_analysis/yt-3.0/commits/8ee5a448456b/
Changeset: 8ee5a448456b
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-05 22:26:38
Summary: Enable deposition into covering grids.
This removes the deposit() method on DataContainers and moves it to the
CoveringGrid. I've added a LeftEdge property as well for inside the
particle_deposit routine. Particle fields can now be obtained for
CoveringGrids, and deposition can be conducted as well.
Affected #: 2 files
diff -r 2dd46f58a1f4e08f360680cb83954ce207873328 -r 8ee5a448456bea82a56021dc11992d3c3d6dd6e0 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -53,6 +53,7 @@
MinimalProjectionData
from yt.utilities.parallel_tools.parallel_analysis_interface import \
parallel_objects, parallel_root_only, ParallelAnalysisInterface
+import yt.geometry.particle_deposit as particle_deposit
from .field_info_container import\
NeedsGridType,\
@@ -433,20 +434,31 @@
fields_to_get = [f for f in fields if f not in self.field_data]
fields_to_get = self._identify_dependencies(fields_to_get)
if len(fields_to_get) == 0: return
- fill, gen = self._split_fields(fields_to_get)
+ fill, gen, part = self._split_fields(fields_to_get)
+ if len(part) > 0: self._fill_particles(part)
if len(fill) > 0: self._fill_fields(fill)
if len(gen) > 0: self._generate_fields(gen)
def _split_fields(self, fields_to_get):
fill, gen = self.pf.h._split_fields(fields_to_get)
+ particles = []
for field in gen:
finfo = self.pf._get_field_info(*field)
try:
finfo.check_available(self)
except NeedsOriginalGrid:
fill.append(field)
+ for field in fill:
+ finfo = self.pf._get_field_info(*field)
+ if finfo.particle_type:
+ particles.append(field)
gen = [f for f in gen if f not in fill]
- return fill, gen
+ fill = [f for f in fill if f not in particles]
+ return fill, gen, particles
+
+ def _fill_particles(self, part):
+ for p in part:
+ self[p] = self._data_source[p]
def _fill_fields(self, fields):
output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
@@ -485,6 +497,20 @@
raise KeyError(field)
return rv
+ @property
+ def LeftEdge(self):
+ return self.left_edge
+
+ def deposit(self, positions, fields = None, method = None):
+ cls = getattr(particle_deposit, "deposit_%s" % method, None)
+ if cls is None:
+ raise YTParticleDepositionNotImplemented(method)
+ op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+ op.initialize()
+ op.process_grid(self, positions, fields)
+ vals = op.finalize()
+ return vals.reshape(self.ActiveDimensions, order="F")
+
class LevelState(object):
current_dx = None
current_dims = None
diff -r 2dd46f58a1f4e08f360680cb83954ce207873328 -r 8ee5a448456bea82a56021dc11992d3c3d6dd6e0 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -519,11 +519,6 @@
if f not in fields_to_generate:
fields_to_generate.append(f)
- def deposit(self, positions, fields, op):
- assert(self._current_chunk.chunk_type == "spatial")
- fields = ensure_list(fields)
- self.hierarchy._deposit_particle_fields(self, positions, fields, op)
-
@contextmanager
def _field_lock(self):
self._locked = True
https://bitbucket.org/yt_analysis/yt-3.0/commits/eb2f5ab4af3f/
Changeset: eb2f5ab4af3f
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-05 22:39:39
Summary: First attempt at porting Bobby's smoothing kernel to particle_deposit.
Note that combined with the covering_grid fixes, one can now apply smoothing
kernels to particles that are collected from diverse regions. I've tested this
for SPH kernels applied to full-domain SPH simulations, and even at relatively
small counts it's unmanageably slow.
Affected #: 2 files
diff -r 8ee5a448456bea82a56021dc11992d3c3d6dd6e0 -r eb2f5ab4af3fe683977390a2fc10943883b158f5 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -29,6 +29,7 @@
import numpy as np
from libc.stdlib cimport malloc, free
cimport cython
+from libc.math cimport sqrt
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -39,6 +40,21 @@
cdef inline int gind(int i, int j, int k, int dims[3]):
return ((k*dims[1])+j)*dims[0]+i
+
+####################################################
+# Standard SPH kernel for use with the Grid method #
+####################################################
+
+cdef inline np.float64_t sph_kernel(np.float64_t x) nogil:
+ cdef np.float64_t kernel
+ if x <= 0.5:
+ kernel = 1.-6.*x*x*(1.-x)
+ elif x>0.5 and x<=1.0:
+ kernel = 2.*(1.-x)*(1.-x)*(1.-x)
+ else:
+ kernel = 0.
+ return kernel
+
cdef class ParticleDepositOperation:
# We assume each will allocate and define their own temporary storage
cdef np.int64_t nvals
diff -r 8ee5a448456bea82a56021dc11992d3c3d6dd6e0 -r eb2f5ab4af3fe683977390a2fc10943883b158f5 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -29,6 +29,7 @@
import numpy as np
from libc.stdlib cimport malloc, free
cimport cython
+from libc.math cimport sqrt
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, \
@@ -104,6 +105,7 @@
left_edge[i] = gobj.LeftEdge[i]
dims[i] = gobj.ActiveDimensions[i]
for i in range(positions.shape[0]):
+ if i % 10000 == 0: print i, positions.shape[0]
# Now we process
for j in range(nf):
field_vals[j] = field_pointers[j][i]
@@ -145,6 +147,68 @@
deposit_count = CountParticles
+cdef class SimpleSmooth(ParticleDepositOperation):
+ # Note that this does nothing at the edges. So it will give a poor
+ # estimate there, and since Octrees are mostly edges, this will be a very
+ # poor SPH kernel.
+ cdef np.float64_t *data
+ cdef public object odata
+ cdef np.float64_t *temp
+ cdef public object otemp
+
+ def initialize(self):
+ self.odata = np.zeros(self.nvals, dtype="float64")
+ cdef np.ndarray arr = self.odata
+ self.data = <np.float64_t*> arr.data
+ self.otemp = np.zeros(self.nvals, dtype="float64")
+ arr = self.otemp
+ self.temp = <np.float64_t*> arr.data
+
+ @cython.cdivision(True)
+ cdef void process(self, int dim[3],
+ np.float64_t left_edge[3],
+ np.float64_t dds[3],
+ np.int64_t offset,
+ np.float64_t ppos[3],
+ np.float64_t *fields
+ ):
+ cdef int ii[3], half_len, ib0[3], ib1[3]
+ cdef int i, j, k
+ cdef np.float64_t idist[3], kernel_sum, dist
+ # Smoothing length is fields[0]
+ kernel_sum = 0.0
+ for i in range(3):
+ ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
+ half_len = <int>(fields[0]/dds[i]) + 1
+ ib0[i] = ii[i] - half_len
+ ib1[i] = ii[i] + half_len
+ if ib0[i] >= dim[i] or ib1[i] <0:
+ return
+ ib0[i] = iclip(ib0[i], 0, dim[i] - 1)
+ ib1[i] = iclip(ib1[i], 0, dim[i] - 1)
+ for i from ib0[0] <= i <= ib1[0]:
+ idist[0] = (ii[0] - i) * (ii[0] - i) * dds[0]
+ for j from ib0[1] <= j <= ib1[1]:
+ idist[1] = (ii[1] - j) * (ii[1] - j) * dds[1]
+ for k from ib0[2] <= k <= ib1[2]:
+ idist[2] = (ii[2] - k) * (ii[2] - k) * dds[2]
+ dist = idist[0] + idist[1] + idist[2]
+ # Calculate distance in multiples of the smoothing length
+ dist = sqrt(dist) / fields[0]
+ self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist)
+ kernel_sum += self.temp[gind(i,j,k,dim) + offset]
+ # Having found the kernel, deposit accordingly into gdata
+ for i from ib0[0] <= i <= ib1[0]:
+ for j from ib0[1] <= j <= ib1[1]:
+ for k from ib0[2] <= k <= ib1[2]:
+ dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
+ self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
+
+ def finalize(self):
+ return self.odata
+
+deposit_simple_smooth = SimpleSmooth
+
cdef class SumParticleField(ParticleDepositOperation):
cdef np.float64_t *sum
cdef public object osum
https://bitbucket.org/yt_analysis/yt-3.0/commits/8c8ef2d6c645/
Changeset: 8c8ef2d6c645
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-05 22:48:42
Summary: Removed print statement.
Affected #: 1 file
diff -r eb2f5ab4af3fe683977390a2fc10943883b158f5 -r 8c8ef2d6c64544d46670ce506b8e8a5d07a004f3 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -105,7 +105,6 @@
left_edge[i] = gobj.LeftEdge[i]
dims[i] = gobj.ActiveDimensions[i]
for i in range(positions.shape[0]):
- if i % 10000 == 0: print i, positions.shape[0]
# Now we process
for j in range(nf):
field_vals[j] = field_pointers[j][i]
https://bitbucket.org/yt_analysis/yt-3.0/commits/5260b868fa0c/
Changeset: 5260b868fa0c
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-07 16:09:18
Summary: Merging in some misc changes.
Affected #: 2 files
diff -r a392d0b00fa31a49aeb9e81c8e629ca64cdf8ca1 -r 5260b868fa0c432905808d6180aee828fd4bdbe8 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -477,8 +477,23 @@
if fields is None: return
fields = self._determine_fields(fields)
# Now we collect all our fields
- fields_to_get = [f for f in fields if f not in self.field_data]
- if len(fields_to_get) == 0:
+ # Here is where we need to perform a validation step, so that if we
+ # have a field requested that we actually *can't* yet get, we put it
+ # off until the end. This prevents double-reading fields that will
+ # need to be used in spatial fields later on.
+ fields_to_get = []
+ # This will be pre-populated with spatial fields
+ fields_to_generate = []
+ for field in self._determine_fields(fields):
+ if field in self.field_data: continue
+ finfo = self.pf._get_field_info(*field)
+ try:
+ finfo.check_available(self)
+ except NeedsGridType:
+ fields_to_generate.append(field)
+ continue
+ fields_to_get.append(field)
+ if len(fields_to_get) == 0 and fields_to_generate == 0:
return
elif self._locked == True:
raise GenerationInProgress(fields)
@@ -502,12 +517,18 @@
read_particles, gen_particles = self.hierarchy._read_particle_fields(
particles, self, self._current_chunk)
self.field_data.update(read_particles)
- fields_to_generate = gen_fluids + gen_particles
+ fields_to_generate += gen_fluids + gen_particles
self._generate_fields(fields_to_generate)
def _generate_fields(self, fields_to_generate):
index = 0
with self._field_lock():
+ # At this point, we assume that any fields that are necessary to
+ # *generate* a field are in fact already available to us. Note
+ # that we do not make any assumption about whether or not the
+ # fields have a spatial requirement. This will be checked inside
+ # _generate_field, at which point additional dependencies may
+ # actually be noted.
while any(f not in self.field_data for f in fields_to_generate):
field = fields_to_generate[index % len(fields_to_generate)]
index += 1
diff -r a392d0b00fa31a49aeb9e81c8e629ca64cdf8ca1 -r 5260b868fa0c432905808d6180aee828fd4bdbe8 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -894,7 +894,7 @@
Child cells are not returned.
"""
_type_name = "data_collection"
- _con_args = ("obj_list",)
+ _con_args = ("_obj_list",)
def __init__(self, center, obj_list, pf = None, field_parameters = None):
YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
self._obj_ids = np.array([o.id - o._id_offset for o in obj_list],
https://bitbucket.org/yt_analysis/yt-3.0/commits/ef8a2cfed946/
Changeset: ef8a2cfed946
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-07 16:41:18
Summary: Merging in changes to the OWLS frontend and simple_smooth.
Affected #: 7 files
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -53,6 +53,7 @@
MinimalProjectionData
from yt.utilities.parallel_tools.parallel_analysis_interface import \
parallel_objects, parallel_root_only, ParallelAnalysisInterface
+import yt.geometry.particle_deposit as particle_deposit
from .field_info_container import\
NeedsGridType,\
@@ -433,20 +434,31 @@
fields_to_get = [f for f in fields if f not in self.field_data]
fields_to_get = self._identify_dependencies(fields_to_get)
if len(fields_to_get) == 0: return
- fill, gen = self._split_fields(fields_to_get)
+ fill, gen, part = self._split_fields(fields_to_get)
+ if len(part) > 0: self._fill_particles(part)
if len(fill) > 0: self._fill_fields(fill)
if len(gen) > 0: self._generate_fields(gen)
def _split_fields(self, fields_to_get):
fill, gen = self.pf.h._split_fields(fields_to_get)
+ particles = []
for field in gen:
finfo = self.pf._get_field_info(*field)
try:
finfo.check_available(self)
except NeedsOriginalGrid:
fill.append(field)
+ for field in fill:
+ finfo = self.pf._get_field_info(*field)
+ if finfo.particle_type:
+ particles.append(field)
gen = [f for f in gen if f not in fill]
- return fill, gen
+ fill = [f for f in fill if f not in particles]
+ return fill, gen, particles
+
+ def _fill_particles(self, part):
+ for p in part:
+ self[p] = self._data_source[p]
def _fill_fields(self, fields):
output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
@@ -485,6 +497,20 @@
raise KeyError(field)
return rv
+ @property
+ def LeftEdge(self):
+ return self.left_edge
+
+ def deposit(self, positions, fields = None, method = None):
+ cls = getattr(particle_deposit, "deposit_%s" % method, None)
+ if cls is None:
+ raise YTParticleDepositionNotImplemented(method)
+ op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+ op.initialize()
+ op.process_grid(self, positions, fields)
+ vals = op.finalize()
+ return vals.reshape(self.ActiveDimensions, order="F")
+
class LevelState(object):
current_dx = None
current_dims = None
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -540,11 +540,6 @@
if f not in fields_to_generate:
fields_to_generate.append(f)
- def deposit(self, positions, fields, op):
- assert(self._current_chunk.chunk_type == "spatial")
- fields = ensure_list(fields)
- self.hierarchy._deposit_particle_fields(self, positions, fields, op)
-
@contextmanager
def _field_lock(self):
self._locked = True
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -163,71 +163,6 @@
for subset in oobjs:
yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
-class OWLSStaticOutput(StaticOutput):
- _hierarchy_class = ParticleGeometryHandler
- _domain_class = ParticleDomainFile
- _fieldinfo_fallback = OWLSFieldInfo
- _fieldinfo_known = KnownOWLSFields
-
- def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
- self._root_dimensions = root_dimensions
- # Set up the template for domain files
- self.storage_filename = None
- super(OWLSStaticOutput, self).__init__(filename, data_style)
-
- def __repr__(self):
- return os.path.basename(self.parameter_filename).split(".")[0]
-
- def _set_units(self):
- self.units = {}
- self.time_units = {}
- self.conversion_factors = {}
- DW = self.domain_right_edge - self.domain_left_edge
- self.units["unitary"] = 1.0 / DW.max()
-
- def _parse_parameter_file(self):
- handle = h5py.File(self.parameter_filename)
- hvals = {}
- hvals.update(handle["/Header"].attrs)
-
- self.dimensionality = 3
- self.refine_by = 2
- self.parameters["HydroMethod"] = "sph"
- self.unique_identifier = \
- int(os.stat(self.parameter_filename)[stat.ST_CTIME])
- # Set standard values
- self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
- self.domain_left_edge = np.zeros(3, "float64")
- self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
- self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
- self.cosmological_simulation = 1
- self.periodicity = (True, True, True)
- self.current_redshift = hvals["Redshift"]
- self.omega_lambda = hvals["OmegaLambda"]
- self.omega_matter = hvals["Omega0"]
- self.hubble_constant = hvals["HubbleParam"]
- self.parameters = hvals
-
- prefix = self.parameter_filename.split(".", 1)[0]
- suffix = self.parameter_filename.rsplit(".", 1)[-1]
- self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
- self.domain_count = hvals["NumFilesPerSnapshot"]
-
- handle.close()
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- try:
- fileh = h5py.File(args[0],'r')
- if "Constants" in fileh["/"].keys() and \
- "Header" in fileh["/"].keys():
- fileh.close()
- return True
- fileh.close()
- except:
- pass
- return False
-
class GadgetBinaryDomainFile(ParticleDomainFile):
def __init__(self, pf, io, domain_filename, domain_id):
with open(domain_filename, "rb") as f:
@@ -404,6 +339,74 @@
# We do not allow load() of these files.
return False
+class OWLSStaticOutput(GadgetStaticOutput):
+ _hierarchy_class = ParticleGeometryHandler
+ _domain_class = ParticleDomainFile
+ _fieldinfo_fallback = OWLSFieldInfo # For now we have separate from Gadget
+ _fieldinfo_known = KnownOWLSFields
+ _header_spec = None # Override so that there's no confusion
+
+ def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
+ self._root_dimensions = root_dimensions
+ # Set up the template for domain files
+ self.storage_filename = None
+ super(OWLSStaticOutput, self).__init__(filename, data_style,
+ root_dimensions,
+ unit_base = None)
+
+ def __repr__(self):
+ return os.path.basename(self.parameter_filename).split(".")[0]
+
+ def _parse_parameter_file(self):
+ handle = h5py.File(self.parameter_filename)
+ hvals = {}
+ hvals.update((str(k), v) for k, v in handle["/Header"].attrs.items())
+
+ self.dimensionality = 3
+ self.refine_by = 2
+ self.parameters["HydroMethod"] = "sph"
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[stat.ST_CTIME])
+
+ # Set standard values
+ self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
+ self.domain_left_edge = np.zeros(3, "float64")
+ self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
+ self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
+ self.cosmological_simulation = 1
+ self.periodicity = (True, True, True)
+ self.current_redshift = hvals["Redshift"]
+ self.omega_lambda = hvals["OmegaLambda"]
+ self.omega_matter = hvals["Omega0"]
+ self.hubble_constant = hvals["HubbleParam"]
+ self.parameters = hvals
+
+ prefix = self.parameter_filename.split(".", 1)[0]
+ suffix = self.parameter_filename.rsplit(".", 1)[-1]
+ self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
+ self.domain_count = hvals["NumFilesPerSnapshot"]
+
+ # To avoid having to open files twice
+ self._unit_base = {}
+ self._unit_base.update((str(k), v) for k, v in handle["/Units"].attrs.items())
+ # Comoving cm is given in the Units
+ self._unit_base['cmcm'] = 1.0 / self._unit_base["UnitLength_in_cm"]
+
+ handle.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = h5py.File(args[0],'r')
+ if "Constants" in fileh["/"].keys() and \
+ "Header" in fileh["/"].keys():
+ fileh.close()
+ return True
+ fileh.close()
+ except:
+ pass
+ return False
+
class TipsyDomainFile(ParticleDomainFile):
def _calculate_offsets(self, field_list):
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -247,3 +247,60 @@
GadgetFieldInfo.add_field(("all", oname + ax), function=func,
particle_type = True)
+# OWLS
+# ====
+
+# I am optimistic that some day we will be able to get rid of much of this, and
+# make OWLS a subclass of Gadget fields.
+
+_owls_ptypes = ("PartType0", "PartType1", "PartType2", "PartType3",
+ "PartType4")
+
+for fname in ["Coordinates", "Velocities", "ParticleIDs",
+ # Note: Mass, not Masses
+ "Mass"]:
+ func = _field_concat(fname)
+ OWLSFieldInfo.add_field(("all", fname), function=func,
+ particle_type = True)
+
+def _owls_particle_fields(ptype):
+ def _Mass(field, data):
+ pind = _owls_ptypes.index(ptype)
+ if data.pf["MassTable"][pind] == 0.0:
+ raise RuntimeError
+ mass = np.ones(data[ptype, "ParticleIDs"].shape[0], dtype="float64")
+ # Note that this is an alias, which is why we need to apply conversion
+ # here. Otherwise we'd have an asymmetry.
+ mass *= data.pf["MassTable"][pind]
+ return mass
+ OWLSFieldInfo.add_field((ptype, "Mass"), function=_Mass,
+ convert_function = _get_conv("mass"),
+ particle_type = True)
+
+for ptype in _owls_ptypes:
+ # Note that this adds a "Known" Mass field and a "Derived" Mass field.
+ # This way the "Known" will get used, and if it's not there, it will use
+ # the derived.
+ KnownOWLSFields.add_field((ptype, "Mass"), function=NullFunc,
+ particle_type = True,
+ convert_function=_get_conv("mass"),
+ units = r"\mathrm{g}")
+ _owls_particle_fields(ptype)
+ KnownOWLSFields.add_field((ptype, "Velocities"), function=NullFunc,
+ particle_type = True,
+ convert_function=_get_conv("velocity"),
+ units = r"\mathrm{cm}/\mathrm{s}")
+ _particle_functions(ptype, "Coordinates", "Mass", OWLSFieldInfo)
+ KnownOWLSFields.add_field((ptype, "Coordinates"), function=NullFunc,
+ particle_type = True)
+_particle_functions("all", "Coordinates", "Mass", OWLSFieldInfo)
+
+# Now we have to manually apply the splits for "all", since we don't want to
+# use the splits defined above.
+
+for iname, oname in [("Coordinates", "particle_position_"),
+ ("Velocities", "particle_velocity_")]:
+ for axi, ax in enumerate("xyz"):
+ func = _field_concat_slice(iname, axi)
+ OWLSFieldInfo.add_field(("all", oname + ax), function=func,
+ particle_type = True)
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -78,7 +78,7 @@
source.quantities["MaxLocation"](field)
mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f",
max_val, mx, my, mz)
- self.pf.parameters["Max%sValue" % (field)] = max_val
- self.pf.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
+ self.pf.parameters["Max%sValue" % (field,)] = max_val
+ self.pf.parameters["Max%sPos" % (field,)] = "%s" % ((mx,my,mz),)
return max_val, np.array((mx,my,mz), dtype='float64')
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -29,6 +29,7 @@
import numpy as np
from libc.stdlib cimport malloc, free
cimport cython
+from libc.math cimport sqrt
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -39,6 +40,21 @@
cdef inline int gind(int i, int j, int k, int dims[3]):
return ((k*dims[1])+j)*dims[0]+i
+
+####################################################
+# Standard SPH kernel for use with the Grid method #
+####################################################
+
+cdef inline np.float64_t sph_kernel(np.float64_t x) nogil:
+ cdef np.float64_t kernel
+ if x <= 0.5:
+ kernel = 1.-6.*x*x*(1.-x)
+ elif x>0.5 and x<=1.0:
+ kernel = 2.*(1.-x)*(1.-x)*(1.-x)
+ else:
+ kernel = 0.
+ return kernel
+
cdef class ParticleDepositOperation:
# We assume each will allocate and define their own temporary storage
cdef np.int64_t nvals
diff -r 5260b868fa0c432905808d6180aee828fd4bdbe8 -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -29,6 +29,7 @@
import numpy as np
from libc.stdlib cimport malloc, free
cimport cython
+from libc.math cimport sqrt
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, \
@@ -145,6 +146,68 @@
deposit_count = CountParticles
+cdef class SimpleSmooth(ParticleDepositOperation):
+ # Note that this does nothing at the edges. So it will give a poor
+ # estimate there, and since Octrees are mostly edges, this will be a very
+ # poor SPH kernel.
+ cdef np.float64_t *data
+ cdef public object odata
+ cdef np.float64_t *temp
+ cdef public object otemp
+
+ def initialize(self):
+ self.odata = np.zeros(self.nvals, dtype="float64")
+ cdef np.ndarray arr = self.odata
+ self.data = <np.float64_t*> arr.data
+ self.otemp = np.zeros(self.nvals, dtype="float64")
+ arr = self.otemp
+ self.temp = <np.float64_t*> arr.data
+
+ @cython.cdivision(True)
+ cdef void process(self, int dim[3],
+ np.float64_t left_edge[3],
+ np.float64_t dds[3],
+ np.int64_t offset,
+ np.float64_t ppos[3],
+ np.float64_t *fields
+ ):
+ cdef int ii[3], half_len, ib0[3], ib1[3]
+ cdef int i, j, k
+ cdef np.float64_t idist[3], kernel_sum, dist
+ # Smoothing length is fields[0]
+ kernel_sum = 0.0
+ for i in range(3):
+ ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
+ half_len = <int>(fields[0]/dds[i]) + 1
+ ib0[i] = ii[i] - half_len
+ ib1[i] = ii[i] + half_len
+ if ib0[i] >= dim[i] or ib1[i] <0:
+ return
+ ib0[i] = iclip(ib0[i], 0, dim[i] - 1)
+ ib1[i] = iclip(ib1[i], 0, dim[i] - 1)
+ for i from ib0[0] <= i <= ib1[0]:
+ idist[0] = (ii[0] - i) * (ii[0] - i) * dds[0]
+ for j from ib0[1] <= j <= ib1[1]:
+ idist[1] = (ii[1] - j) * (ii[1] - j) * dds[1]
+ for k from ib0[2] <= k <= ib1[2]:
+ idist[2] = (ii[2] - k) * (ii[2] - k) * dds[2]
+ dist = idist[0] + idist[1] + idist[2]
+ # Calculate distance in multiples of the smoothing length
+ dist = sqrt(dist) / fields[0]
+ self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist)
+ kernel_sum += self.temp[gind(i,j,k,dim) + offset]
+ # Having found the kernel, deposit accordingly into gdata
+ for i from ib0[0] <= i <= ib1[0]:
+ for j from ib0[1] <= j <= ib1[1]:
+ for k from ib0[2] <= k <= ib1[2]:
+ dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
+ self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
+
+ def finalize(self):
+ return self.odata
+
+deposit_simple_smooth = SimpleSmooth
+
cdef class SumParticleField(ParticleDepositOperation):
cdef np.float64_t *sum
cdef public object osum
https://bitbucket.org/yt_analysis/yt-3.0/commits/570136cac0ee/
Changeset: 570136cac0ee
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-07 18:45:55
Summary: Adding arbitrary_grid object.
This object will return particles as would a .region, but will also be able to
deposit fields from those particles as well. It will not return fluid fields
that would require any fields on disk.
Affected #: 1 file
diff -r ef8a2cfed946999667f5e45eb7cea5d4634d36ab -r 570136cac0ee3abd07d3310f09b41819e38e844a yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -511,6 +511,53 @@
vals = op.finalize()
return vals.reshape(self.ActiveDimensions, order="F")
+class YTArbitraryGridBase(YTCoveringGridBase):
+ """A 3D region with arbitrary bounds and dimensions.
+
+ In contrast to the Covering Grid, this object accepts a left edge, a right
+ edge, and dimensions. This allows it to be used for creating 3D particle
+ deposition fields that are independent of the underlying mesh, whether that
+ is yt-generated or from the simulation data. For example, arbitrary boxes
+ around particles can be drawn and particle deposition fields can be
+ created. This object will refuse to generate any fluid fields.
+
+ Parameters
+ ----------
+ left_edge : array_like
+ The left edge of the region to be extracted
+ rigth_edge : array_like
+ The left edge of the region to be extracted
+ dims : array_like
+ Number of cells along each axis of resulting grid.
+
+ Examples
+ --------
+ >>> obj = pf.h.arbitrary_grid([0.0, 0.0, 0.0], [0.99, 0.99, 0.99],
+ ... dims=[128, 128, 128])
+ """
+ _spatial = True
+ _type_name = "arbitrary_grid"
+ _con_args = ('left_edge', 'right_edge', 'ActiveDimensions')
+ _container_fields = ("dx", "dy", "dz", "x", "y", "z")
+ def __init__(self, left_edge, right_edge, dims,
+ pf = None, field_parameters = None):
+ if field_parameters is None:
+ center = None
+ else:
+ center = field_parameters.get("center", None)
+ YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
+ self.left_edge = np.array(left_edge)
+ self.right_edge = np.array(right_edge)
+ self.ActiveDimensions = np.array(dims, dtype='int32')
+ if self.ActiveDimensions.size == 1:
+ self.ActiveDimensions = np.array([dims, dims, dims], dtype="int32")
+ self.dds = (self.right_edge - self.left_edge)/self.ActiveDimensions
+ self.level = 99
+ self._setup_data_source()
+
+ def _fill_fields(self, fields):
+ raise NotImplementedError
+
class LevelState(object):
current_dx = None
current_dims = None
https://bitbucket.org/yt_analysis/yt-3.0/commits/5e472283912f/
Changeset: 5e472283912f
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-07 18:46:27
Summary: Adding CIC deposit method to particle_deposit.
This gives bitwise identical results to particle_density in those frontends
that support it.
Affected #: 3 files
diff -r 570136cac0ee3abd07d3310f09b41819e38e844a -r 5e472283912f5343de5a875e2a3f68d7f907ef6c yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -696,3 +696,15 @@
units = r"\mathrm{g}/\mathrm{cm}^{3}",
projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
projection_conversion = 'cm')
+
+def particle_cic(field, data):
+ pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
+ d = data.deposit(pos, [data["ParticleMass"]], method="cic")
+ d /= data["CellVolume"]
+ return d
+add_field(("deposit", "all_cic"), function=particle_cic,
+ validators=[ValidateSpatial()],
+ display_name = "\\mathrm{All CIC Density}",
+ units = r"\mathrm{g}/\mathrm{cm}^{3}",
+ projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+ projection_conversion = 'cm')
diff -r 570136cac0ee3abd07d3310f09b41819e38e844a -r 5e472283912f5343de5a875e2a3f68d7f907ef6c yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -96,6 +96,20 @@
projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
projection_conversion = 'cm')
+ def particle_cic(field, data):
+ pos = data[ptype, coord_name]
+ d = data.deposit(pos, [data[ptype, mass_name]], method = "cic")
+ d /= data["CellVolume"]
+ return d
+
+ registry.add_field(("deposit", "%s_cic" % ptype),
+ function = particle_cic,
+ validators = [ValidateSpatial()],
+ display_name = "\\mathrm{%s CIC Density}" % ptype,
+ units = r"\mathrm{g}/\mathrm{cm}^{3}",
+ projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+ projection_conversion = 'cm')
+
# Now some translation functions.
registry.add_field((ptype, "ParticleMass"),
diff -r 570136cac0ee3abd07d3310f09b41819e38e844a -r 5e472283912f5343de5a875e2a3f68d7f907ef6c yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -295,3 +295,44 @@
deposit_std = StdParticleField
+cdef class CICDeposit(ParticleDepositOperation):
+ cdef np.float64_t *field
+ cdef public object ofield
+ def initialize(self):
+ self.ofield = np.zeros(self.nvals, dtype="float64")
+ cdef np.ndarray arr = self.ofield
+ self.field = <np.float64_t *> arr.data
+
+ cdef void process(self, int dim[3],
+ np.float64_t left_edge[3],
+ np.float64_t dds[3],
+ np.int64_t offset, # offset into IO field
+ np.float64_t ppos[3], # this particle's position
+ np.float64_t *fields # any other fields we need
+ ):
+
+ cdef int i, j, k, ind[3], ii
+ cdef np.float64_t rpos[3], rdds[3][2]
+ cdef np.float64_t fact, edge0, edge1, edge2
+ cdef np.float64_t le0, le1, le2
+ cdef np.float64_t dx, dy, dz, dx2, dy2, dz2
+
+ # Compute the position of the central cell
+ for i in range(3):
+ rpos[i] = (ppos[i]-left_edge[i])/dds[i]
+ rpos[i] = fclip(rpos[i], 0.5001, dim[i]-0.5001)
+ ind[i] = <int> (rpos[i] + 0.5)
+ # Note these are 1, then 0
+ rdds[i][1] = (<np.float64_t> ind[i]) + 0.5 - rpos[i]
+ rdds[i][0] = 1.0 - rdds[i][1]
+
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ ii = gind(ind[0] - i, ind[1] - j, ind[2] - k, dim) + offset
+ self.field[ii] += fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
+
+ def finalize(self):
+ return self.ofield
+
+deposit_cic = CICDeposit
https://bitbucket.org/yt_analysis/yt-3.0/commits/ea5291e8fa1b/
Changeset: ea5291e8fa1b
Branch: yt-3.0
User: MatthewTurk
Date: 2013-06-11 15:32:40
Summary: Merged in MatthewTurk/yt-3.0 (pull request #46)
OWLS, kernel smoothing and arbitrary grids
Affected #: 9 files
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -53,6 +53,7 @@
MinimalProjectionData
from yt.utilities.parallel_tools.parallel_analysis_interface import \
parallel_objects, parallel_root_only, ParallelAnalysisInterface
+import yt.geometry.particle_deposit as particle_deposit
from .field_info_container import\
NeedsGridType,\
@@ -433,20 +434,31 @@
fields_to_get = [f for f in fields if f not in self.field_data]
fields_to_get = self._identify_dependencies(fields_to_get)
if len(fields_to_get) == 0: return
- fill, gen = self._split_fields(fields_to_get)
+ fill, gen, part = self._split_fields(fields_to_get)
+ if len(part) > 0: self._fill_particles(part)
if len(fill) > 0: self._fill_fields(fill)
if len(gen) > 0: self._generate_fields(gen)
def _split_fields(self, fields_to_get):
fill, gen = self.pf.h._split_fields(fields_to_get)
+ particles = []
for field in gen:
finfo = self.pf._get_field_info(*field)
try:
finfo.check_available(self)
except NeedsOriginalGrid:
fill.append(field)
+ for field in fill:
+ finfo = self.pf._get_field_info(*field)
+ if finfo.particle_type:
+ particles.append(field)
gen = [f for f in gen if f not in fill]
- return fill, gen
+ fill = [f for f in fill if f not in particles]
+ return fill, gen, particles
+
+ def _fill_particles(self, part):
+ for p in part:
+ self[p] = self._data_source[p]
def _fill_fields(self, fields):
output_fields = [np.zeros(self.ActiveDimensions, dtype="float64")
@@ -485,6 +497,67 @@
raise KeyError(field)
return rv
+ @property
+ def LeftEdge(self):
+ return self.left_edge
+
+ def deposit(self, positions, fields = None, method = None):
+ cls = getattr(particle_deposit, "deposit_%s" % method, None)
+ if cls is None:
+ raise YTParticleDepositionNotImplemented(method)
+ op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+ op.initialize()
+ op.process_grid(self, positions, fields)
+ vals = op.finalize()
+ return vals.reshape(self.ActiveDimensions, order="F")
+
+class YTArbitraryGridBase(YTCoveringGridBase):
+ """A 3D region with arbitrary bounds and dimensions.
+
+ In contrast to the Covering Grid, this object accepts a left edge, a right
+ edge, and dimensions. This allows it to be used for creating 3D particle
+ deposition fields that are independent of the underlying mesh, whether that
+ is yt-generated or from the simulation data. For example, arbitrary boxes
+ around particles can be drawn and particle deposition fields can be
+ created. This object will refuse to generate any fluid fields.
+
+ Parameters
+ ----------
+ left_edge : array_like
+ The left edge of the region to be extracted
+ rigth_edge : array_like
+ The left edge of the region to be extracted
+ dims : array_like
+ Number of cells along each axis of resulting grid.
+
+ Examples
+ --------
+ >>> obj = pf.h.arbitrary_grid([0.0, 0.0, 0.0], [0.99, 0.99, 0.99],
+ ... dims=[128, 128, 128])
+ """
+ _spatial = True
+ _type_name = "arbitrary_grid"
+ _con_args = ('left_edge', 'right_edge', 'ActiveDimensions')
+ _container_fields = ("dx", "dy", "dz", "x", "y", "z")
+ def __init__(self, left_edge, right_edge, dims,
+ pf = None, field_parameters = None):
+ if field_parameters is None:
+ center = None
+ else:
+ center = field_parameters.get("center", None)
+ YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
+ self.left_edge = np.array(left_edge)
+ self.right_edge = np.array(right_edge)
+ self.ActiveDimensions = np.array(dims, dtype='int32')
+ if self.ActiveDimensions.size == 1:
+ self.ActiveDimensions = np.array([dims, dims, dims], dtype="int32")
+ self.dds = (self.right_edge - self.left_edge)/self.ActiveDimensions
+ self.level = 99
+ self._setup_data_source()
+
+ def _fill_fields(self, fields):
+ raise NotImplementedError
+
class LevelState(object):
current_dx = None
current_dims = None
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/data_objects/data_containers.py
--- a/yt/data_objects/data_containers.py
+++ b/yt/data_objects/data_containers.py
@@ -477,8 +477,23 @@
if fields is None: return
fields = self._determine_fields(fields)
# Now we collect all our fields
- fields_to_get = [f for f in fields if f not in self.field_data]
- if len(fields_to_get) == 0:
+ # Here is where we need to perform a validation step, so that if we
+ # have a field requested that we actually *can't* yet get, we put it
+ # off until the end. This prevents double-reading fields that will
+ # need to be used in spatial fields later on.
+ fields_to_get = []
+ # This will be pre-populated with spatial fields
+ fields_to_generate = []
+ for field in self._determine_fields(fields):
+ if field in self.field_data: continue
+ finfo = self.pf._get_field_info(*field)
+ try:
+ finfo.check_available(self)
+ except NeedsGridType:
+ fields_to_generate.append(field)
+ continue
+ fields_to_get.append(field)
+ if len(fields_to_get) == 0 and fields_to_generate == 0:
return
elif self._locked == True:
raise GenerationInProgress(fields)
@@ -502,12 +517,18 @@
read_particles, gen_particles = self.hierarchy._read_particle_fields(
particles, self, self._current_chunk)
self.field_data.update(read_particles)
- fields_to_generate = gen_fluids + gen_particles
+ fields_to_generate += gen_fluids + gen_particles
self._generate_fields(fields_to_generate)
def _generate_fields(self, fields_to_generate):
index = 0
with self._field_lock():
+ # At this point, we assume that any fields that are necessary to
+ # *generate* a field are in fact already available to us. Note
+ # that we do not make any assumption about whether or not the
+ # fields have a spatial requirement. This will be checked inside
+ # _generate_field, at which point additional dependencies may
+ # actually be noted.
while any(f not in self.field_data for f in fields_to_generate):
field = fields_to_generate[index % len(fields_to_generate)]
index += 1
@@ -519,11 +540,6 @@
if f not in fields_to_generate:
fields_to_generate.append(f)
- def deposit(self, positions, fields, op):
- assert(self._current_chunk.chunk_type == "spatial")
- fields = ensure_list(fields)
- self.hierarchy._deposit_particle_fields(self, positions, fields, op)
-
@contextmanager
def _field_lock(self):
self._locked = True
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/data_objects/selection_data_containers.py
--- a/yt/data_objects/selection_data_containers.py
+++ b/yt/data_objects/selection_data_containers.py
@@ -894,7 +894,7 @@
Child cells are not returned.
"""
_type_name = "data_collection"
- _con_args = ("obj_list",)
+ _con_args = ("_obj_list",)
def __init__(self, center, obj_list, pf = None, field_parameters = None):
YTSelectionContainer3D.__init__(self, center, pf, field_parameters)
self._obj_ids = np.array([o.id - o._id_offset for o in obj_list],
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/frontends/enzo/fields.py
--- a/yt/frontends/enzo/fields.py
+++ b/yt/frontends/enzo/fields.py
@@ -696,3 +696,15 @@
units = r"\mathrm{g}/\mathrm{cm}^{3}",
projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
projection_conversion = 'cm')
+
+def particle_cic(field, data):
+ pos = np.column_stack([data["particle_position_%s" % ax] for ax in 'xyz'])
+ d = data.deposit(pos, [data["ParticleMass"]], method="cic")
+ d /= data["CellVolume"]
+ return d
+add_field(("deposit", "all_cic"), function=particle_cic,
+ validators=[ValidateSpatial()],
+ display_name = "\\mathrm{All CIC Density}",
+ units = r"\mathrm{g}/\mathrm{cm}^{3}",
+ projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+ projection_conversion = 'cm')
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/frontends/sph/data_structures.py
--- a/yt/frontends/sph/data_structures.py
+++ b/yt/frontends/sph/data_structures.py
@@ -163,71 +163,6 @@
for subset in oobjs:
yield YTDataChunk(dobj, "io", [subset], subset.cell_count)
-class OWLSStaticOutput(StaticOutput):
- _hierarchy_class = ParticleGeometryHandler
- _domain_class = ParticleDomainFile
- _fieldinfo_fallback = OWLSFieldInfo
- _fieldinfo_known = KnownOWLSFields
-
- def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
- self._root_dimensions = root_dimensions
- # Set up the template for domain files
- self.storage_filename = None
- super(OWLSStaticOutput, self).__init__(filename, data_style)
-
- def __repr__(self):
- return os.path.basename(self.parameter_filename).split(".")[0]
-
- def _set_units(self):
- self.units = {}
- self.time_units = {}
- self.conversion_factors = {}
- DW = self.domain_right_edge - self.domain_left_edge
- self.units["unitary"] = 1.0 / DW.max()
-
- def _parse_parameter_file(self):
- handle = h5py.File(self.parameter_filename)
- hvals = {}
- hvals.update(handle["/Header"].attrs)
-
- self.dimensionality = 3
- self.refine_by = 2
- self.parameters["HydroMethod"] = "sph"
- self.unique_identifier = \
- int(os.stat(self.parameter_filename)[stat.ST_CTIME])
- # Set standard values
- self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
- self.domain_left_edge = np.zeros(3, "float64")
- self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
- self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
- self.cosmological_simulation = 1
- self.periodicity = (True, True, True)
- self.current_redshift = hvals["Redshift"]
- self.omega_lambda = hvals["OmegaLambda"]
- self.omega_matter = hvals["Omega0"]
- self.hubble_constant = hvals["HubbleParam"]
- self.parameters = hvals
-
- prefix = self.parameter_filename.split(".", 1)[0]
- suffix = self.parameter_filename.rsplit(".", 1)[-1]
- self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
- self.domain_count = hvals["NumFilesPerSnapshot"]
-
- handle.close()
-
- @classmethod
- def _is_valid(self, *args, **kwargs):
- try:
- fileh = h5py.File(args[0],'r')
- if "Constants" in fileh["/"].keys() and \
- "Header" in fileh["/"].keys():
- fileh.close()
- return True
- fileh.close()
- except:
- pass
- return False
-
class GadgetBinaryDomainFile(ParticleDomainFile):
def __init__(self, pf, io, domain_filename, domain_id):
with open(domain_filename, "rb") as f:
@@ -404,6 +339,74 @@
# We do not allow load() of these files.
return False
+class OWLSStaticOutput(GadgetStaticOutput):
+ _hierarchy_class = ParticleGeometryHandler
+ _domain_class = ParticleDomainFile
+ _fieldinfo_fallback = OWLSFieldInfo # For now we have separate from Gadget
+ _fieldinfo_known = KnownOWLSFields
+ _header_spec = None # Override so that there's no confusion
+
+ def __init__(self, filename, data_style="OWLS", root_dimensions = 64):
+ self._root_dimensions = root_dimensions
+ # Set up the template for domain files
+ self.storage_filename = None
+ super(OWLSStaticOutput, self).__init__(filename, data_style,
+ root_dimensions,
+ unit_base = None)
+
+ def __repr__(self):
+ return os.path.basename(self.parameter_filename).split(".")[0]
+
+ def _parse_parameter_file(self):
+ handle = h5py.File(self.parameter_filename)
+ hvals = {}
+ hvals.update((str(k), v) for k, v in handle["/Header"].attrs.items())
+
+ self.dimensionality = 3
+ self.refine_by = 2
+ self.parameters["HydroMethod"] = "sph"
+ self.unique_identifier = \
+ int(os.stat(self.parameter_filename)[stat.ST_CTIME])
+
+ # Set standard values
+ self.current_time = hvals["Time_GYR"] * sec_conversion["Gyr"]
+ self.domain_left_edge = np.zeros(3, "float64")
+ self.domain_right_edge = np.ones(3, "float64") * hvals["BoxSize"]
+ self.domain_dimensions = np.ones(3, "int32") * self._root_dimensions
+ self.cosmological_simulation = 1
+ self.periodicity = (True, True, True)
+ self.current_redshift = hvals["Redshift"]
+ self.omega_lambda = hvals["OmegaLambda"]
+ self.omega_matter = hvals["Omega0"]
+ self.hubble_constant = hvals["HubbleParam"]
+ self.parameters = hvals
+
+ prefix = self.parameter_filename.split(".", 1)[0]
+ suffix = self.parameter_filename.rsplit(".", 1)[-1]
+ self.domain_template = "%s.%%(num)i.%s" % (prefix, suffix)
+ self.domain_count = hvals["NumFilesPerSnapshot"]
+
+ # To avoid having to open files twice
+ self._unit_base = {}
+ self._unit_base.update((str(k), v) for k, v in handle["/Units"].attrs.items())
+ # Comoving cm is given in the Units
+ self._unit_base['cmcm'] = 1.0 / self._unit_base["UnitLength_in_cm"]
+
+ handle.close()
+
+ @classmethod
+ def _is_valid(self, *args, **kwargs):
+ try:
+ fileh = h5py.File(args[0],'r')
+ if "Constants" in fileh["/"].keys() and \
+ "Header" in fileh["/"].keys():
+ fileh.close()
+ return True
+ fileh.close()
+ except:
+ pass
+ return False
+
class TipsyDomainFile(ParticleDomainFile):
def _calculate_offsets(self, field_list):
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/frontends/sph/fields.py
--- a/yt/frontends/sph/fields.py
+++ b/yt/frontends/sph/fields.py
@@ -96,6 +96,20 @@
projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
projection_conversion = 'cm')
+ def particle_cic(field, data):
+ pos = data[ptype, coord_name]
+ d = data.deposit(pos, [data[ptype, mass_name]], method = "cic")
+ d /= data["CellVolume"]
+ return d
+
+ registry.add_field(("deposit", "%s_cic" % ptype),
+ function = particle_cic,
+ validators = [ValidateSpatial()],
+ display_name = "\\mathrm{%s CIC Density}" % ptype,
+ units = r"\mathrm{g}/\mathrm{cm}^{3}",
+ projected_units = r"\mathrm{g}/\mathrm{cm}^{2}",
+ projection_conversion = 'cm')
+
# Now some translation functions.
registry.add_field((ptype, "ParticleMass"),
@@ -247,3 +261,60 @@
GadgetFieldInfo.add_field(("all", oname + ax), function=func,
particle_type = True)
+# OWLS
+# ====
+
+# I am optimistic that some day we will be able to get rid of much of this, and
+# make OWLS a subclass of Gadget fields.
+
+_owls_ptypes = ("PartType0", "PartType1", "PartType2", "PartType3",
+ "PartType4")
+
+for fname in ["Coordinates", "Velocities", "ParticleIDs",
+ # Note: Mass, not Masses
+ "Mass"]:
+ func = _field_concat(fname)
+ OWLSFieldInfo.add_field(("all", fname), function=func,
+ particle_type = True)
+
+def _owls_particle_fields(ptype):
+ def _Mass(field, data):
+ pind = _owls_ptypes.index(ptype)
+ if data.pf["MassTable"][pind] == 0.0:
+ raise RuntimeError
+ mass = np.ones(data[ptype, "ParticleIDs"].shape[0], dtype="float64")
+ # Note that this is an alias, which is why we need to apply conversion
+ # here. Otherwise we'd have an asymmetry.
+ mass *= data.pf["MassTable"][pind]
+ return mass
+ OWLSFieldInfo.add_field((ptype, "Mass"), function=_Mass,
+ convert_function = _get_conv("mass"),
+ particle_type = True)
+
+for ptype in _owls_ptypes:
+ # Note that this adds a "Known" Mass field and a "Derived" Mass field.
+ # This way the "Known" will get used, and if it's not there, it will use
+ # the derived.
+ KnownOWLSFields.add_field((ptype, "Mass"), function=NullFunc,
+ particle_type = True,
+ convert_function=_get_conv("mass"),
+ units = r"\mathrm{g}")
+ _owls_particle_fields(ptype)
+ KnownOWLSFields.add_field((ptype, "Velocities"), function=NullFunc,
+ particle_type = True,
+ convert_function=_get_conv("velocity"),
+ units = r"\mathrm{cm}/\mathrm{s}")
+ _particle_functions(ptype, "Coordinates", "Mass", OWLSFieldInfo)
+ KnownOWLSFields.add_field((ptype, "Coordinates"), function=NullFunc,
+ particle_type = True)
+_particle_functions("all", "Coordinates", "Mass", OWLSFieldInfo)
+
+# Now we have to manually apply the splits for "all", since we don't want to
+# use the splits defined above.
+
+for iname, oname in [("Coordinates", "particle_position_"),
+ ("Velocities", "particle_velocity_")]:
+ for axi, ax in enumerate("xyz"):
+ func = _field_concat_slice(iname, axi)
+ OWLSFieldInfo.add_field(("all", oname + ax), function=func,
+ particle_type = True)
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/geometry/oct_geometry_handler.py
--- a/yt/geometry/oct_geometry_handler.py
+++ b/yt/geometry/oct_geometry_handler.py
@@ -78,7 +78,7 @@
source.quantities["MaxLocation"](field)
mylog.info("Max Value is %0.5e at %0.16f %0.16f %0.16f",
max_val, mx, my, mz)
- self.pf.parameters["Max%sValue" % (field)] = max_val
- self.pf.parameters["Max%sPos" % (field)] = "%s" % ((mx,my,mz),)
+ self.pf.parameters["Max%sValue" % (field,)] = max_val
+ self.pf.parameters["Max%sPos" % (field,)] = "%s" % ((mx,my,mz),)
return max_val, np.array((mx,my,mz), dtype='float64')
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -29,6 +29,7 @@
import numpy as np
from libc.stdlib cimport malloc, free
cimport cython
+from libc.math cimport sqrt
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
@@ -39,6 +40,21 @@
cdef inline int gind(int i, int j, int k, int dims[3]):
return ((k*dims[1])+j)*dims[0]+i
+
+####################################################
+# Standard SPH kernel for use with the Grid method #
+####################################################
+
+cdef inline np.float64_t sph_kernel(np.float64_t x) nogil:
+ cdef np.float64_t kernel
+ if x <= 0.5:
+ kernel = 1.-6.*x*x*(1.-x)
+ elif x>0.5 and x<=1.0:
+ kernel = 2.*(1.-x)*(1.-x)*(1.-x)
+ else:
+ kernel = 0.
+ return kernel
+
cdef class ParticleDepositOperation:
# We assume each will allocate and define their own temporary storage
cdef np.int64_t nvals
diff -r 7a79dd9c7ebefddba48dfe0440a47ca3fddd6b98 -r ea5291e8fa1bf7647c3dbdd0052dcb68444fb260 yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -29,6 +29,7 @@
import numpy as np
from libc.stdlib cimport malloc, free
cimport cython
+from libc.math cimport sqrt
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, \
@@ -145,6 +146,68 @@
deposit_count = CountParticles
+cdef class SimpleSmooth(ParticleDepositOperation):
+ # Note that this does nothing at the edges. So it will give a poor
+ # estimate there, and since Octrees are mostly edges, this will be a very
+ # poor SPH kernel.
+ cdef np.float64_t *data
+ cdef public object odata
+ cdef np.float64_t *temp
+ cdef public object otemp
+
+ def initialize(self):
+ self.odata = np.zeros(self.nvals, dtype="float64")
+ cdef np.ndarray arr = self.odata
+ self.data = <np.float64_t*> arr.data
+ self.otemp = np.zeros(self.nvals, dtype="float64")
+ arr = self.otemp
+ self.temp = <np.float64_t*> arr.data
+
+ @cython.cdivision(True)
+ cdef void process(self, int dim[3],
+ np.float64_t left_edge[3],
+ np.float64_t dds[3],
+ np.int64_t offset,
+ np.float64_t ppos[3],
+ np.float64_t *fields
+ ):
+ cdef int ii[3], half_len, ib0[3], ib1[3]
+ cdef int i, j, k
+ cdef np.float64_t idist[3], kernel_sum, dist
+ # Smoothing length is fields[0]
+ kernel_sum = 0.0
+ for i in range(3):
+ ii[i] = <int>((ppos[i] - left_edge[i])/dds[i])
+ half_len = <int>(fields[0]/dds[i]) + 1
+ ib0[i] = ii[i] - half_len
+ ib1[i] = ii[i] + half_len
+ if ib0[i] >= dim[i] or ib1[i] <0:
+ return
+ ib0[i] = iclip(ib0[i], 0, dim[i] - 1)
+ ib1[i] = iclip(ib1[i], 0, dim[i] - 1)
+ for i from ib0[0] <= i <= ib1[0]:
+ idist[0] = (ii[0] - i) * (ii[0] - i) * dds[0]
+ for j from ib0[1] <= j <= ib1[1]:
+ idist[1] = (ii[1] - j) * (ii[1] - j) * dds[1]
+ for k from ib0[2] <= k <= ib1[2]:
+ idist[2] = (ii[2] - k) * (ii[2] - k) * dds[2]
+ dist = idist[0] + idist[1] + idist[2]
+ # Calculate distance in multiples of the smoothing length
+ dist = sqrt(dist) / fields[0]
+ self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist)
+ kernel_sum += self.temp[gind(i,j,k,dim) + offset]
+ # Having found the kernel, deposit accordingly into gdata
+ for i from ib0[0] <= i <= ib1[0]:
+ for j from ib0[1] <= j <= ib1[1]:
+ for k from ib0[2] <= k <= ib1[2]:
+ dist = self.temp[gind(i,j,k,dim) + offset] / kernel_sum
+ self.data[gind(i,j,k,dim) + offset] += fields[1] * dist
+
+ def finalize(self):
+ return self.odata
+
+deposit_simple_smooth = SimpleSmooth
+
cdef class SumParticleField(ParticleDepositOperation):
cdef np.float64_t *sum
cdef public object osum
@@ -232,3 +295,44 @@
deposit_std = StdParticleField
+cdef class CICDeposit(ParticleDepositOperation):
+ cdef np.float64_t *field
+ cdef public object ofield
+ def initialize(self):
+ self.ofield = np.zeros(self.nvals, dtype="float64")
+ cdef np.ndarray arr = self.ofield
+ self.field = <np.float64_t *> arr.data
+
+ cdef void process(self, int dim[3],
+ np.float64_t left_edge[3],
+ np.float64_t dds[3],
+ np.int64_t offset, # offset into IO field
+ np.float64_t ppos[3], # this particle's position
+ np.float64_t *fields # any other fields we need
+ ):
+
+ cdef int i, j, k, ind[3], ii
+ cdef np.float64_t rpos[3], rdds[3][2]
+ cdef np.float64_t fact, edge0, edge1, edge2
+ cdef np.float64_t le0, le1, le2
+ cdef np.float64_t dx, dy, dz, dx2, dy2, dz2
+
+ # Compute the position of the central cell
+ for i in range(3):
+ rpos[i] = (ppos[i]-left_edge[i])/dds[i]
+ rpos[i] = fclip(rpos[i], 0.5001, dim[i]-0.5001)
+ ind[i] = <int> (rpos[i] + 0.5)
+ # Note these are 1, then 0
+ rdds[i][1] = (<np.float64_t> ind[i]) + 0.5 - rpos[i]
+ rdds[i][0] = 1.0 - rdds[i][1]
+
+ for i in range(2):
+ for j in range(2):
+ for k in range(2):
+ ii = gind(ind[0] - i, ind[1] - j, ind[2] - k, dim) + offset
+ self.field[ii] += fields[0]*rdds[0][i]*rdds[1][j]*rdds[2][k]
+
+ def finalize(self):
+ return self.ofield
+
+deposit_cic = CICDeposit
Repository URL: https://bitbucket.org/yt_analysis/yt-3.0/
--
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