Merge pull request #2610 from cphyc/indep-octree-raytracing

Volume rendering for octrees

appveyor.yml

@@ -11,7 +11,7 @@ environment:
       - PYTHON_VERSION: "3.8"
 
 platform:
-    -x64
+    - x64
 
 install:
     - "if not exist \"%userprofile%\\.config\\yt\" mkdir %userprofile%\\.config\\yt"

doc/source/cookbook/render_two_fields.py

@@ -1,5 +1,5 @@
 import yt
-from yt.visualization.volume_rendering.api import Scene, VolumeSource
+from yt.visualization.volume_rendering.api import Scene, create_volume_source
 
 filePath = "Sedov_3d/sedov_hdf5_chk_0003"
 ds = yt.load(filePath)
@@ -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

@@ -1,5 +1,5 @@
 import yt
-from yt.visualization.volume_rendering.api import Scene, VolumeSource
+from yt.visualization.volume_rendering.api import Scene, create_volume_source
 import numpy as np
 
 field = ("gas", "density")
@@ -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

@@ -1,6 +1,7 @@
 import glob
 import os
 import sys
+from distutils.ccompiler import get_default_compiler
 from distutils.version import LooseVersion
 
 import pkg_resources
@@ -45,6 +46,11 @@
 else:
     std_libs = ["m"]
 
+if get_default_compiler() == "msvc":
+    CPP14_FLAG = ["/std:c++14"]
+else:
+    CPP14_FLAG = ["--std=c++14"]
+
 cythonize_aliases = {
     "LIB_DIR": "yt/utilities/lib/",
     "LIB_DIR_EWAH": ["yt/utilities/lib/", "yt/utilities/lib/ewahboolarray/"],
@@ -56,8 +62,9 @@
     ],
     "STD_LIBS": std_libs,
     "OMP_ARGS": omp_args,
-    "FIXED_INTERP": "yt/utilities/lib/fixed_interpolator.c",
+    "FIXED_INTERP": "yt/utilities/lib/fixed_interpolator.cpp",
     "ARTIO_SOURCE": glob.glob("yt/frontends/artio/artio_headers/*.c"),
+    "CPP14_FLAG": CPP14_FLAG,
 }
 
 lib_exts = [

yt/data_objects/index_subobjects/octree_subset.py

@@ -6,7 +6,7 @@
 import yt.geometry.particle_deposit as particle_deposit
 import yt.geometry.particle_smooth as particle_smooth
 from yt.data_objects.data_containers import YTSelectionContainer
-from yt.funcs import mylog
+from yt.funcs import issue_deprecation_warning, mylog
 from yt.geometry.particle_oct_container import ParticleOctreeContainer
 from yt.units.dimensions import length
 from yt.units.yt_array import YTArray
@@ -522,6 +522,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."
+            )
+            issue_deprecation_warning(message)
+            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:
     def __init__(self, ind, block_slice):

yt/frontends/ramses/data_structures.py

@@ -276,14 +276,19 @@ def _fill_with_ghostzones(
             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(
-            selector, self.domain_id, cell_count
-        )
+        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:

yt/frontends/stream/data_structures.py

@@ -646,7 +646,10 @@ class StreamOctreeSubset(OctreeSubset):
     domain_id = 1
     _domain_offset = 1
 
-    def __init__(self, base_region, ds, oct_handler, over_refine_factor=1):
+    def __init__(
+        self, base_region, ds, oct_handler, over_refine_factor=1, num_ghost_zones=0
+    ):
+        self._over_refine_factor = over_refine_factor
         self._num_zones = 1 << (over_refine_factor)
         self.field_data = YTFieldData()
         self.field_parameters = {}
@@ -659,7 +662,33 @@ def __init__(self, base_region, ds, oct_handler, over_refine_factor=1):
         self.base_region = base_region
         self.base_selector = base_region.selector
 
-    def fill(self, content, dest, selector, offset):
+        self._num_ghost_zones = num_ghost_zones
+
+        if num_ghost_zones > 0:
+            if not all(ds.periodicity):
+                mylog.warn("Ghost zones will wrongly assume the domain to be periodic.")
+            base_grid = StreamOctreeSubset(
+                base_region, ds, oct_handler, over_refine_factor
+            )
+            self._base_grid = base_grid
+
+    def retrieve_ghost_zones(self, ngz, fields, smoothed=False):
+        try:
+            new_subset = self._subset_with_gz
+            mylog.debug("Reusing previous subset with ghost zone.")
+        except AttributeError:
+            new_subset = StreamOctreeSubset(
+                self.base_region,
+                self.ds,
+                self.oct_handler,
+                self._over_refine_factor,
+                num_ghost_zones=ngz,
+            )
+            self._subset_with_gz = new_subset
+
+        return new_subset
+
+    def _fill_no_ghostzones(self, content, dest, selector, offset):
         # Here we get a copy of the file, which we skip through and read the
         # bits we want.
         oct_handler = self.oct_handler
@@ -675,6 +704,37 @@ def fill(self, content, dest, selector, offset):
         )
         return count
 
+    def _fill_with_ghostzones(self, content, dest, selector, offset):
+        oct_handler = self.oct_handler
+        ndim = self.ds.dimensionality
+        cell_count = (
+            selector.count_octs(self.oct_handler, self.domain_id) * self.nz ** ndim
+        )
+
+        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,
+            ) = oct_handler.file_index_octs_with_ghost_zones(
+                selector, self.domain_id, cell_count
+            )
+            self._ghost_zone_cache = gz_cache
+        levels[:] = 0
+        dest.update((field, np.empty(cell_count, dtype="float64")) for field in content)
+        # Make references ...
+        oct_handler.fill_level(0, levels, cell_inds, file_inds, dest, content, offset)
+
+    def fill(self, content, dest, selector, offset):
+        if self._num_ghost_zones == 0:
+            return self._fill_no_ghostzones(content, dest, selector, offset)
+        else:
+            return self._fill_with_ghostzones(content, dest, selector, offset)
+
 
 class StreamOctreeHandler(OctreeIndex):
     def __init__(self, ds, dataset_type=None):
@@ -752,6 +812,22 @@ class StreamOctreeDataset(StreamDataset):
     _field_info_class = StreamFieldInfo
     _dataset_type = "stream_octree"
 
+    levelmax = None
+
+    def __init__(
+        self,
+        stream_handler,
+        storage_filename=None,
+        geometry="cartesian",
+        unit_system="cgs",
+    ):
+        super(StreamOctreeDataset, self).__init__(
+            stream_handler, storage_filename, geometry, unit_system
+        )
+        # Set up levelmax
+        self.max_level = stream_handler.levels.max()
+        self.min_level = stream_handler.levels.min()
+
 
 class StreamUnstructuredMesh(UnstructuredMesh):
     _index_offset = 0

yt/geometry/oct_container.pxd

@@ -74,7 +74,7 @@ cdef class OctreeContainer:
                         OctVisitor visitor,
                         int vc = ?, np.int64_t *indices = ?)
     cdef Oct *next_root(self, int domain_id, int ind[3])
-    cdef Oct *next_child(self, int domain_id, int ind[3], Oct *parent)
+    cdef Oct *next_child(self, int domain_id, int ind[3], Oct *parent) except? NULL
     cdef void append_domain(self, np.int64_t domain_count)
     # The fill_style is the ordering, C or F, of the octs in the file.  "o"
     # corresponds to C, and "r" is for Fortran.

yt/geometry/oct_container.pyx

@@ -643,7 +643,7 @@ cdef class OctreeContainer:
         self.nocts += 1
         return next
 
-    cdef Oct* next_child(self, int domain_id, int ind[3], Oct *parent):
+    cdef Oct* next_child(self, int domain_id, int ind[3], Oct *parent) except? NULL:
         cdef int i
         cdef Oct *next = NULL
         if parent.children != NULL:

yt/utilities/answer_testing/utils.py

@@ -19,7 +19,7 @@
 from yt.loaders import load, load_simulation
 from yt.units.yt_array import YTArray, YTQuantity
 from yt.visualization import particle_plots, plot_window as pw, profile_plotter
-from yt.visualization.volume_rendering.scene import Scene
+from yt.visualization.volume_rendering.api import Scene
 
 
 def _streamline_for_io(params):