Ray object dl for SPH data: interpolate kernel table for (b/hsml)**2 instead of b/hsml

yt/data_objects/selection_objects/ray.py

@@ -222,6 +222,7 @@ def _generate_container_field_sph(self, field):
 
  # Use an interpolation table to evaluate the integrated 2D
  # kernel from the dimensionless impact parameter b/hsml.
+ # The table tabulates integrals for values of (b/hsml)**2
  itab = SPHKernelInterpolationTable(self.ds.kernel_name)
- dl = itab.interpolate_array(b / hsml) * mass / dens / hsml**2
+ dl = itab.interpolate_array((b / hsml) ** 2) * mass / dens / hsml**2
  return dl / length

yt/data_objects/tests/test_rays.py

@@ -4,7 +4,10 @@
 from yt import load
 from yt.testing import (
  assert_rel_equal,
+ cubicspline_python,
  fake_random_ds,
+ fake_sph_grid_ds,
+ integrate_kernel,
  requires_file,
  requires_module,
 )
@@ -71,3 +74,80 @@ def test_ray_particle():
  ray = ds.ray(ds.domain_left_edge, ds.domain_right_edge)
  assert_equal(ray["t"].shape, (1451,))
  assert ray["dts"].sum(dtype="f8") > 0
+
+
+## test that kernels integrate correctly
+# (1) including the right particles
+# (2) correct t and dts values for those particles
+## fake_sph_grid_ds:
+# This dataset should have 27 particles with the particles arranged
+# uniformly on a 3D grid. The bottom left corner is (0.5,0.5,0.5) and
+# the top right corner is (2.5,2.5,2.5). All particles will have
+# non-overlapping smoothing regions with a radius of 0.05, masses of
+# 1, and densities of 1, and zero velocity.
+
+
+def test_ray_particle2():
+ kernelfunc = cubicspline_python
+ ds = fake_sph_grid_ds(hsml_factor=1.0)
+ ds.kernel_name = "cubic"
+
+ ## Ray through the one particle at (0.5, 0.5, 0.5):
+ ## test basic kernel integration
+ eps = 0.0 # 1e-7
+ start0 = np.array((1.0 + eps, 0.0, 0.5))
+ end0 = np.array((0.0, 1.0 + eps, 0.5))
+ ray0 = ds.ray(start0, end0)
+ b0 = np.array([np.sqrt(2.0) * eps])
+ hsml0 = np.array([0.05])
+ len0 = np.sqrt(np.sum((end0 - start0) ** 2))
+ # for a ParticleDataset like this one, the Ray object attempts
+ # to generate the 't' and 'dts' fields using the grid method
+ ray0.field_data["t"] = ray0.ds.arr(ray0._generate_container_field_sph("t"))
+ ray0.field_data["dts"] = ray0.ds.arr(ray0._generate_container_field_sph("dts"))
+ # not demanding too much precision;
+ # from kernel volume integrals, the linear interpolation
+ # restricts you to 4 -- 5 digits precision
+ assert_equal(ray0["t"].shape, (1,))
+ assert_rel_equal(ray0["t"], np.array([0.5]), 5)
+ assert_rel_equal(ray0[("gas", "position")].v, np.array([[0.5, 0.5, 0.5]]), 5)
+ dl0 = integrate_kernel(kernelfunc, b0, hsml0)
+ dl0 *= ray0[("gas", "mass")].v / ray0[("gas", "density")].v
+ assert_rel_equal(ray0[("dts")].v, dl0 / len0, 4)
+
+ ## Ray in the middle of the box:
+ ## test end points, >1 particle
+ start1 = np.array((1.53, 0.53, 1.0))
+ end1 = np.array((1.53, 0.53, 3.0))
+ ray1 = ds.ray(start1, end1)
+ b1 = np.array([np.sqrt(2.0) * 0.03] * 2)
+ hsml1 = np.array([0.05] * 2)
+ len1 = np.sqrt(np.sum((end1 - start1) ** 2))
+ # for a ParticleDataset like this one, the Ray object attempts
+ # to generate the 't' and 'dts' fields using the grid method
+ ray1.field_data["t"] = ray1.ds.arr(ray1._generate_container_field_sph("t"))
+ ray1.field_data["dts"] = ray1.ds.arr(ray1._generate_container_field_sph("dts"))
+ # not demanding too much precision;
+ # from kernel volume integrals, the linear interpolation
+ # restricts you to 4 -- 5 digits precision
+ assert_equal(ray1["t"].shape, (2,))
+ assert_rel_equal(ray1["t"], np.array([0.25, 0.75]), 5)
+ assert_rel_equal(
+ ray1[("gas", "position")].v, np.array([[1.5, 0.5, 1.5], [1.5, 0.5, 2.5]]), 5
+ )
+ dl1 = integrate_kernel(kernelfunc, b1, hsml1)
+ dl1 *= ray1[("gas", "mass")].v / ray1[("gas", "density")].v
+ assert_rel_equal(ray1[("dts")].v, dl1 / len1, 4)
+
+ ## Ray missing all particles:
+ ## test handling of size-0 selections
+ start2 = np.array((1.0, 2.0, 0.0))
+ end2 = np.array((1.0, 2.0, 3.0))
+ ray2 = ds.ray(start2, end2)
+ # for a ParticleDataset like this one, the Ray object attempts
+ # to generate the 't' and 'dts' fields using the grid method
+ ray2.field_data["t"] = ray2.ds.arr(ray2._generate_container_field_sph("t"))
+ ray2.field_data["dts"] = ray2.ds.arr(ray2._generate_container_field_sph("dts"))
+ assert_equal(ray2["t"].shape, (0,))
+ assert_equal(ray2["dts"].shape, (0,))
+ assert_equal(ray2[("gas", "position")].v.shape, (0, 3))