Merge branch 'indep-octree-raytracing' into support-vertex-centred-data

doc/source/cookbook/render_two_fields.py

@@ -15,13 +15,13 @@
 cam.switch_orientation()
 
 # add rendering of density field
-dens = VolumeSource(ds, field='dens')
+dens = create_volume_source(ds, field='dens')
 dens.use_ghost_zones = True
 sc.add_source(dens)
 sc.save('density.png', sigma_clip=6)
 
 # add rendering of x-velocity field
-vel = VolumeSource(ds, field='velx')
+vel = create_volume_source(ds, field='velx')
 vel.use_ghost_zones = True
 sc.add_source(vel)
 sc.save('density_any_velocity.png', sigma_clip=6)

doc/source/cookbook/various_lens.py

@@ -14,7 +14,7 @@
 # Follow the simple_volume_rendering cookbook for the first part of this.
 ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
 sc = Scene()
-vol = VolumeSource(ds, field=field)
+vol = create_volume_source(ds, field=field)
 tf = vol.transfer_function
 tf.grey_opacity = True
 

setup.py

@@ -153,7 +153,7 @@ def _compile(
               ],
               libraries=std_libs,
               language="c++",
-              extra_compile_arg=["-std=c++03"]),
+              extra_compile_args=["-std=c++03"]),
     Extension("yt.utilities.lib.cykdtree.utils",
               [
                   "yt/utilities/lib/cykdtree/utils.pyx",
@@ -162,7 +162,7 @@ def _compile(
               depends=["yt/utilities/lib/cykdtree/c_utils.hpp"],
               libraries=std_libs,
               language="c++",
-              extra_compile_arg=["-std=c++03"]),    
+              extra_compile_args=["-std=c++03"]),
     Extension("yt.utilities.lib.fnv_hash",
               ["yt/utilities/lib/fnv_hash.pyx"],
               include_dirs=["yt/utilities/lib/"],
@@ -174,10 +174,11 @@ def _compile(
               libraries=std_libs),
     Extension("yt.utilities.lib.marching_cubes",
               ["yt/utilities/lib/marching_cubes.pyx",
-               "yt/utilities/lib/fixed_interpolator.c"],
+               "yt/utilities/lib/fixed_interpolator.cpp"],
               include_dirs=["yt/utilities/lib/"],
               libraries=std_libs,
-              depends=["yt/utilities/lib/fixed_interpolator.h"]),
+              language="c++",
+              depends=["yt/utilities/lib/fixed_interpolator.hpp"]),
     Extension("yt.utilities.lib.mesh_triangulation",
               ["yt/utilities/lib/mesh_triangulation.pyx"],
               depends=["yt/utilities/lib/mesh_triangulation.h"]),
@@ -201,6 +202,14 @@ def _compile(
               extra_link_args=omp_args,
               libraries=std_libs,
               language='c++'),
+    Extension("yt.utilities.lib.cyoctree_raytracing",
+              ["yt/utilities/lib/cyoctree_raytracing.pyx"],
+              include_dirs=["yt/utilities/lib"],
+              extra_compile_args=omp_args + ["-std=c++11"],
+              extra_link_args=omp_args,
+              libraries=std_libs,
+              language='c++',
+              depends=["yt/utilities/lib/octree_raytracing.hpp"]),
     Extension("yt.utilities.lib.primitives",
               ["yt/utilities/lib/primitives.pyx"],
               libraries=std_libs),
@@ -213,29 +222,32 @@ def _compile(
               depends=["yt/utilities/lib/origami_tags.h"]),
     Extension("yt.utilities.lib.grid_traversal",
               ["yt/utilities/lib/grid_traversal.pyx",
-               "yt/utilities/lib/fixed_interpolator.c"],
+               "yt/utilities/lib/fixed_interpolator.cpp"],
               include_dirs=["yt/utilities/lib/"],
               libraries=std_libs,
-              depends=["yt/utilities/lib/fixed_interpolator.h"]),
+              language="c++",
+              depends=["yt/utilities/lib/fixed_interpolator.hpp"]),
     Extension("yt.utilities.lib.ewah_bool_wrap",
               ["yt/utilities/lib/ewah_bool_wrap.pyx"],
               include_dirs=["yt/utilities/lib/",
                             "yt/utilities/lib/ewahboolarray"],
               language="c++"),
     Extension("yt.utilities.lib.image_samplers",
               ["yt/utilities/lib/image_samplers.pyx",
-               "yt/utilities/lib/fixed_interpolator.c"],
+               "yt/utilities/lib/fixed_interpolator.cpp"],
               include_dirs=["yt/utilities/lib/"],
               libraries=std_libs,
               extra_compile_args=omp_args,
               extra_link_args=omp_args,
-              depends=["yt/utilities/lib/fixed_interpolator.h"]),
+              language="c++",
+              depends=["yt/utilities/lib/fixed_interpolator.hpp"]),
     Extension("yt.utilities.lib.partitioned_grid",
               ["yt/utilities/lib/partitioned_grid.pyx",
-               "yt/utilities/lib/fixed_interpolator.c"],
+               "yt/utilities/lib/fixed_interpolator.cpp"],
               include_dirs=["yt/utilities/lib/"],
               libraries=std_libs,
-              depends=["yt/utilities/lib/fixed_interpolator.h"]),
+              language="c++",
+              depends=["yt/utilities/lib/fixed_interpolator.hpp"]),
     Extension("yt.utilities.lib.element_mappings",
               ["yt/utilities/lib/element_mappings.pyx"],
               libraries=std_libs),

yt/data_objects/octree_subset.py

@@ -1,3 +1,4 @@
+import warnings
 from contextlib import contextmanager
 import numpy as np
 
@@ -391,6 +392,35 @@ def select_particles(self, selector, x, y, z):
         mask = selector.select_points(x,y,z, 0.0)
         return mask
 
+    def get_vertex_centered_data(self, fields):
+        _old_api = isinstance(fields, (str, tuple))
+        if _old_api:
+            message = (
+                'get_vertex_centered_data() requires list of fields, rather than '
+                'a single field as an argument.'
+            )
+            warnings.warn(message, DeprecationWarning, stacklevel=2)
+            fields = [fields]
+
+        # Make sure the field list has only unique entries
+        fields = list(set(fields))
+        new_fields = {}
+        cg = self.retrieve_ghost_zones(1, fields)
+        for field in fields:
+            new_fields[field] = cg[field][1: ,1: ,1: ].copy()
+            np.add(new_fields[field], cg[field][:-1,1: ,1: ], new_fields[field])
+            np.add(new_fields[field], cg[field][1: ,:-1,1: ], new_fields[field])
+            np.add(new_fields[field], cg[field][1: ,1: ,:-1], new_fields[field])
+            np.add(new_fields[field], cg[field][:-1,1: ,:-1], new_fields[field])
+            np.add(new_fields[field], cg[field][1: ,:-1,:-1], new_fields[field])
+            np.add(new_fields[field], cg[field][:-1,:-1,1: ], new_fields[field])
+            np.add(new_fields[field], cg[field][:-1,:-1,:-1], new_fields[field])
+            np.multiply(new_fields[field], 0.125, new_fields[field])
+
+        if _old_api:
+            return new_fields[fields[0]]
+        return new_fields
+
 class OctreeSubsetBlockSlicePosition(object):
     def __init__(self, ind, block_slice):
         self.ind = ind
@@ -440,7 +470,6 @@ def get_vertex_centered_data(self, fields, smoothed=False, no_ghost=False):
         new_field = self.block_slice.get_vertex_centered_data(fields)[field]
         return {field: new_field[..., self.ind]}
 
-
     @contextmanager
     def _field_parameter_state(self, field_parameters):
         yield self.block_slice.octree_subset._field_parameter_state(

yt/fields/geometric_fields.py

@@ -93,11 +93,11 @@ def _zeros(field, data):
 
     def _ones(field, data):
         """Returns one for all cells"""
-        arr = np.ones(data.ires.shape, dtype="float64")
-        tmp = data.apply_units(arr, field.units)
+        tmp = np.ones(data.ires.shape, dtype="float64")
+        arr = data.apply_units(tmp, field.units)
         if data._spatial:
-            return data._reshape_vals(tmp)
-        return tmp
+            return data._reshape_vals(arr)
+        return arr
 
     registry.add_field(("index", "ones"),
                        sampling_type="cell",

yt/frontends/ramses/data_structures.py

@@ -234,8 +234,13 @@ def _fill_with_ghostzones(self, fd, fields, selector, file_handler, num_ghost_zo
 
         cell_count = selector.count_octs(self.oct_handler, self.domain_id) * self.nz**ndim
 
-        levels, cell_inds, file_inds, domains = self.oct_handler.file_index_octs_with_ghost_zones(
+        gz_cache = getattr(self, "_ghost_zone_cache", None)
+        if gz_cache:
+            levels, cell_inds, file_inds, domains = gz_cache
+        else:
+            gz_cache = levels, cell_inds, file_inds, domains = self.oct_handler.file_index_octs_with_ghost_zones(
             selector, self.domain_id, cell_count)
+            self._ghost_zone_cache = gz_cache
 
         # Initializing data container
         for field in fields:
@@ -291,8 +296,12 @@ def fill(self, fd, fields, selector, file_handler):
     def retrieve_ghost_zones(self, ngz, fields, smoothed=False):
         if smoothed is True:
             raise NotImplementedError
-
-        new_subset = RAMSESDomainSubset(self.base_region, self.domain, self.ds, num_ghost_zones=ngz, base_grid=self)
+        new_subset = getattr(self, '_subset_with_gz', None)
+        if not new_subset:
+            new_subset = RAMSESDomainSubset(self.base_region, self.domain, self.ds, num_ghost_zones=ngz, base_grid=self)
+        else:
+            mylog.debug('Reusing previous subset with ghost zone for domain %s' % self.domain_id)
+        self._subset_with_gz = new_subset
 
         # Cache the fields
         new_subset.get_data(fields)

yt/utilities/lib/cyoctree_raytracing.pxd

@@ -0,0 +1,30 @@
+"""This is a wrapper around the C++ class to efficiently cast rays into an octree.
+It relies on the seminal paper by  J. Revelles,, C.Ureña and M.Lastra.
+"""
+
+
+cimport numpy as np
+import numpy as np
+from libcpp.vector cimport vector
+cimport cython
+from cython.parallel import prange, parallel
+from libc.stdlib cimport free, malloc
+
+from .image_samplers cimport ImageAccumulator, ImageSampler
+from .grid_traversal cimport sampler_function
+from .volume_container cimport VolumeContainer
+from .partitioned_grid cimport PartitionedGrid
+
+cdef extern from "octree_raytracing.cpp":
+    cdef cppclass RayInfo[T]:
+        vector[T] keys
+        vector[double] t
+
+    cdef cppclass Octree3D[T] nogil:
+        Octree3D(int depth, double* LE, double* RE)
+        void insert_node_no_ret(const int* ipos, const int lvl, T key)
+        void cast_ray(double* origins, double* directions, vector[T] keyList, vector[double] tList)
+
+cdef class CythonOctreeRayTracing:
+    cdef Octree3D[int]* oct
+    cdef int depth

yt/utilities/lib/cyoctree_raytracing.pyx

@@ -0,0 +1,41 @@
+"""This is a wrapper around the C++ class to efficiently cast rays into an octree.
+It relies on the seminal paper by  J. Revelles,, C.Ureña and M.Lastra.
+"""
+
+
+cimport numpy as np
+import numpy as np
+from libcpp.vector cimport vector
+cimport cython
+from cython.parallel import prange, parallel
+from libc.stdlib cimport free, malloc
+
+from .image_samplers cimport ImageAccumulator, ImageSampler
+from .grid_traversal cimport sampler_function
+from .volume_container cimport VolumeContainer
+from .partitioned_grid cimport PartitionedGrid
+
+DEF Nch = 4
+
+
+cdef class CythonOctreeRayTracing:
+    def __init__(self, np.ndarray LE, np.ndarray RE, int depth):
+        cdef double* LE_ptr = <double *>LE.data
+        cdef double* RE_ptr = <double *>RE.data
+        self.oct = new Octree3D[int](depth, LE_ptr, RE_ptr)
+        self.depth = depth
+
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    def add_nodes(self, int[:, :] ipos_view, int[:] lvl_view, int[:] key):
+        cdef int i
+        cdef int ii[3]
+
+        for i in range(len(key)):
+            ii[0] = ipos_view[i, 0]
+            ii[1] = ipos_view[i, 1]
+            ii[2] = ipos_view[i, 2]
+            self.oct.insert_node_no_ret(ii, lvl_view[i], <int> key[i])
+
+    def __dealloc__(self):
+        del self.oct