diff --git a/doc/source/analyzing/generating_processed_data.rst b/doc/source/analyzing/generating_processed_data.rst index 84abf87eb7..5f748d9e25 100644 --- a/doc/source/analyzing/generating_processed_data.rst +++ b/doc/source/analyzing/generating_processed_data.rst @@ -54,7 +54,8 @@ the transformation of a variable mesh of points consisting of positions and sizes into a fixed-size array that appears like an image. This process is that of pixelization, which yt handles transparently internally. You can access this functionality by constructing a -:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` and supplying +:class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` +and supplying to it your :class:`~yt.data_objects.data_containers.YTSelectionContainer2D` object, as well as some information about how you want the final image to look. You can specify both the bounds of the image (in the appropriate x-y plane) and @@ -62,8 +63,8 @@ the resolution of the output image. You can then have yt pixelize any field you like. .. note:: In previous versions of yt, there was a special class of - FixedResolutionBuffer for off-axis slices. This is no longer - necessary. + FixedResolutionBuffer for off-axis slices. This is still used + for off-axis SPH data projections: OffAxisFixedResolutionBuffer. To create :class:`~yt.data_objects.data_containers.YTSelectionContainer2D` objects, you can access them as described in :ref:`data-objects`, specifically the section @@ -99,7 +100,10 @@ this, see :ref:`saving-grid-data-containers`. In the FITS case, there is an option for setting the ``units`` of the coordinate system in the file. If you want to overwrite a file with the same name, set ``clobber=True``. -The :class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` can even be exported +The :class:`~yt.visualization.fixed_resolution.FixedResolutionBuffer` +(and its +:class:`~yt.visualization.fixed_resolution.OffAxisProjectionFixedResolutionBuffer` +subclass) can even be exported as a 2D dataset itself, which may be operated on in the same way as any other dataset in yt: .. code-block:: python diff --git a/doc/source/quickstart/4)_Data_Objects_and_Time_Series.ipynb b/doc/source/quickstart/4)_Data_Objects_and_Time_Series.ipynb index 49592f19f1..351719dce7 100644 --- a/doc/source/quickstart/4)_Data_Objects_and_Time_Series.ipynb +++ b/doc/source/quickstart/4)_Data_Objects_and_Time_Series.ipynb @@ -337,6 +337,8 @@ "\n", "There are two different types of covering grids: unsmoothed and smoothed. Smoothed grids will be filled through a cascading interpolation process; they will be filled at level 0, interpolated to level 1, filled at level 1, interpolated to level 2, filled at level 2, etc. This will help to reduce edge effects. Unsmoothed covering grids will not be interpolated, but rather values will be duplicated multiple times.\n", "\n", + "For SPH datasets, the covering grid gives the SPH-interpolated value of a field at each grid cell center. This is done for unsmoothed grids; smoothed grids are not available for SPH data.\n", + "\n", "Here we create an unsmoothed covering grid at level 2, with the left edge at `[0.0, 0.0, 0.0]` and with dimensions equal to those that would cover the entire domain at level 2. We can then ask for the Density field, which will be a 3D array." ] }, @@ -385,13 +387,13 @@ ], "metadata": { "kernelspec": { - "name": "python3", "display_name": "Python 3.9.5 64-bit ('yt-dev': pyenv)", "metadata": { "interpreter": { "hash": "14363bd97bed451d1329fb3e06aa057a9e955a9421c5343dd7530f5497723a41" } - } + }, + "name": "python3" }, "language_info": { "codemirror_mode": { diff --git a/doc/source/visualizing/plots.rst b/doc/source/visualizing/plots.rst index ec7b6109a5..72e5aac758 100644 --- a/doc/source/visualizing/plots.rst +++ b/doc/source/visualizing/plots.rst @@ -293,8 +293,8 @@ argument. Optionally, a ``north_vector`` can be specified to fix the orientation of the image plane. .. note:: Not every data types have support for off-axis slices yet. - Currently, this operation is supported for grid based data with cartesian geometry. - In some cases (like SPH data) an off-axis projection over a thin region might be used instead. + Currently, this operation is supported for grid based and SPH data with cartesian geometry. + In some cases an off-axis projection over a thin region might be used instead. .. _projection-plots: @@ -433,6 +433,8 @@ by applying the In this use case, the volume renderer casts a set of plane parallel rays, one for each pixel in the image. The data values along each ray are summed, creating the final image buffer. +For SPH datsets, the coordinates are instead simply rotated before the axis-aligned +projection function is applied. .. _off-axis-projection-function: @@ -652,10 +654,6 @@ simply pass ``all`` as the first argument of the field tuple: Additional Notes for Plotting Particle Data ------------------------------------------- -An important caveat when visualizing particle data is that off-axis slice plotting is -not available for any particle data. However, axis-aligned slice plots (as described in -:ref:`slice-plots`) will work. - Since version 4.2.0, off-axis projections ares supported for non-SPH particle data. Previous to that, this operation was only supported for SPH particles. Two historical workaround methods were available for plotting non-SPH particles with off-axis diff --git a/nose_ignores.txt b/nose_ignores.txt index 115b917d59..b5529a6ecc 100644 --- a/nose_ignores.txt +++ b/nose_ignores.txt @@ -45,3 +45,5 @@ --ignore-file=test_set_log_level\.py --ignore-file=test_field_parsing\.py --ignore-file=test_disks\.py +--ignore-file=test_offaxisprojection_pytestonly\.py +--ignore-file=test_sph_pixelization_pytestonly\.py diff --git a/tests/tests.yaml b/tests/tests.yaml index d0f12ad1e5..40836b3381 100644 --- a/tests/tests.yaml +++ b/tests/tests.yaml @@ -37,7 +37,7 @@ answer_tests: - yt/frontends/amrvac/tests/test_outputs.py:test_riemann_cartesian_175D - yt/frontends/amrvac/tests/test_outputs.py:test_rmi_cartesian_dust_2D - local_arepo_011: # PR 4419 + local_arepo_013: # PR 4939 - yt/frontends/arepo/tests/test_outputs.py:test_arepo_bullet - yt/frontends/arepo/tests/test_outputs.py:test_arepo_tng59 - yt/frontends/arepo/tests/test_outputs.py:test_arepo_cr @@ -92,7 +92,7 @@ answer_tests: - yt/frontends/flash/tests/test_outputs.py:test_wind_tunnel - yt/frontends/flash/tests/test_outputs.py:test_fid_1to3_b1 - local_gadget_009: # PR 3258 + local_gadget_010: # PR 4939 - yt/frontends/gadget/tests/test_outputs.py:test_iso_collapse - yt/frontends/gadget/tests/test_outputs.py:test_pid_uniqueness - yt/frontends/gadget/tests/test_outputs.py:test_bigendian_field_access @@ -108,7 +108,7 @@ answer_tests: local_gdf_002: - yt/frontends/gdf/tests/test_outputs_nose.py:test_sedov_tunnel - local_gizmo_008: # PR 2909 + local_gizmo_009: # PR 4939 - yt/frontends/gizmo/tests/test_outputs.py:test_gizmo_64 local_halos_012: # PR 3325 @@ -118,7 +118,7 @@ answer_tests: - yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g5 - yt/frontends/gadget_fof/tests/test_outputs.py:test_fields_g42 - local_owls_008: # PR 2909 + local_owls_009: # PR 4939 - yt/frontends/owls/tests/test_outputs.py:test_snapshot_033 - yt/frontends/owls/tests/test_outputs.py:test_OWLS_particlefilter @@ -130,7 +130,7 @@ answer_tests: - yt/visualization/tests/test_particle_plot.py:test_particle_phase_answers - yt/visualization/tests/test_callbacks.py:test_axis_manipulations - local_tipsy_009: # PR 2909 + local_tipsy_010: # PR 4939 - yt/frontends/tipsy/tests/test_outputs.py:test_pkdgrav - yt/frontends/tipsy/tests/test_outputs.py:test_gasoline_dmonly - yt/frontends/tipsy/tests/test_outputs.py:test_tipsy_galaxy @@ -221,5 +221,7 @@ other_tests: - "--exclude-test=yt.frontends.gdf.tests.test_outputs.TestGDF" - "--exclude-test=yt.frontends.adaptahop.tests.test_outputs" - "--exclude-test=yt.frontends.stream.tests.test_stream_particles.test_stream_non_cartesian_particles" + - "--ignore-file=test_offaxisprojection_pytestonly\\.py" + - "--ignore-file=test_sph_pixelization_pytestonly\\.py" cookbook: - 'doc/source/cookbook/tests/test_cookbook.py' diff --git a/yt/data_objects/construction_data_containers.py b/yt/data_objects/construction_data_containers.py index 03c10f92c1..5f0c23ad58 100644 --- a/yt/data_objects/construction_data_containers.py +++ b/yt/data_objects/construction_data_containers.py @@ -649,6 +649,8 @@ class YTCoveringGrid(YTSelectionContainer3D): level : int The resolution level data to which data will be gridded. Level 0 is the root grid dx for that dataset. + (The grid resolution will be simulation size / 2**level along + each grid axis.) left_edge : array_like The left edge of the region to be extracted. Specify units by supplying a YTArray, otherwise code length units are assumed. @@ -995,14 +997,15 @@ def _fill_sph_particles(self, fields): smoothing_style = getattr(self.ds, "sph_smoothing_style", "scatter") normalize = getattr(self.ds, "use_sph_normalization", True) + kernel_name = getattr(self.ds, "kernel_name", "cubic") bounds, size = self._get_grid_bounds_size() period = self.ds.coordinates.period.copy() if hasattr(period, "in_units"): period = period.in_units("code_length").d - # TODO maybe there is a better way of handling this - is_periodic = int(any(self.ds.periodicity)) + # check periodicity per dimension + is_periodic = self.ds.periodicity if smoothing_style == "scatter": for field in fields: @@ -1036,6 +1039,7 @@ def _fill_sph_particles(self, fields): pbar=pbar, check_period=is_periodic, period=period, + kernel_name=kernel_name, ) if normalize: pixelize_sph_kernel_arbitrary_grid( @@ -1051,6 +1055,7 @@ def _fill_sph_particles(self, fields): pbar=pbar, check_period=is_periodic, period=period, + kernel_name=kernel_name, ) if normalize: diff --git a/yt/geometry/coordinates/cartesian_coordinates.py b/yt/geometry/coordinates/cartesian_coordinates.py index 775de8c6a5..4fed1977c4 100644 --- a/yt/geometry/coordinates/cartesian_coordinates.py +++ b/yt/geometry/coordinates/cartesian_coordinates.py @@ -12,6 +12,7 @@ pixelize_element_mesh, pixelize_element_mesh_line, pixelize_off_axis_cartesian, + pixelize_sph_kernel_cutting, pixelize_sph_kernel_projection, pixelize_sph_kernel_slice, ) @@ -323,11 +324,20 @@ def _ortho_pixelize( # We should be using fcoords field = data_source._determine_fields(field)[0] finfo = data_source.ds.field_info[field] - period = self.period[:2].copy() # dummy here - period[0] = self.period[self.x_axis[dim]] - period[1] = self.period[self.y_axis[dim]] - if hasattr(period, "in_units"): - period = period.in_units("code_length").d + # some coordinate handlers use only projection-plane periods, + # others need all box periods. + period2 = self.period[:2].copy() # dummy here + period2[0] = self.period[self.x_axis[dim]] + period2[1] = self.period[self.y_axis[dim]] + period3 = self.period[:].copy() # dummy here + period3[0] = self.period[self.x_axis[dim]] + period3[1] = self.period[self.y_axis[dim]] + zax = list({0, 1, 2} - {self.x_axis[dim], self.y_axis[dim]})[0] + period3[2] = self.period[zax] + if hasattr(period2, "in_units"): + period2 = period2.in_units("code_length").d + if hasattr(period3, "in_units"): + period3 = period3.in_units("code_length").d buff = np.full((size[1], size[0]), np.nan, dtype="float64") particle_datasets = (ParticleDataset, StreamParticlesDataset) @@ -349,26 +359,49 @@ def _ortho_pixelize( coord, bounds, int(antialias), - period, + period2, int(periodic), return_mask=True, ) elif isinstance(data_source.ds, particle_datasets) and is_sph_field: + # SPH handling ptype = field[0] if ptype == "gas": ptype = data_source.ds._sph_ptypes[0] px_name = self.axis_name[self.x_axis[dim]] py_name = self.axis_name[self.y_axis[dim]] + # need z coordinates for depth, + # but name isn't saved in the handler -> use the 'other one' + pz_name = list(set(self.axis_order) - {px_name, py_name})[0] + + # ignore default True periodic argument + # (not actually supplied by a call from + # FixedResolutionBuffer), and use the dataset periodicity + # instead + xa = self.x_axis[dim] + ya = self.y_axis[dim] + # axorder = data_source.ds.coordinates.axis_order + za = list({0, 1, 2} - {xa, ya})[0] + ds_periodic = data_source.ds.periodicity + _periodic = np.array(ds_periodic) + _periodic[0] = ds_periodic[xa] + _periodic[1] = ds_periodic[ya] + _periodic[2] = ds_periodic[za] ounits = data_source.ds.field_info[field].output_units bnds = data_source.ds.arr(bounds, "code_length").tolist() - if isinstance(data_source, YTParticleProj): + kernel_name = None + if hasattr(data_source.ds, "kernel_name"): + kernel_name = data_source.ds.kernel_name + if kernel_name is None: + kernel_name = "cubic" + + if isinstance(data_source, YTParticleProj): # projection weight = data_source.weight_field moment = data_source.moment le, re = data_source.data_source.get_bbox() - xa = self.x_axis[dim] - ya = self.y_axis[dim] # If we're not periodic, we need to clip to the boundary edges # or we get errors about extending off the edge of the region. + # (depth/z range is handled by region setting) if not self.ds.periodicity[xa]: le[xa] = max(bounds[0], self.ds.domain_left_edge[xa]) re[xa] = min(bounds[1], self.ds.domain_right_edge[xa]) @@ -389,6 +422,10 @@ def _ortho_pixelize( data_source=data_source.data_source, ) proj_reg.set_field_parameter("axis", data_source.axis) + # need some z bounds for SPH projection + # -> use source bounds + bnds3 = bnds + [le[za], re[za]] + buff = np.zeros(size, dtype="float64") mask_uint8 = np.zeros_like(buff, dtype="uint8") if weight is None: @@ -399,13 +436,15 @@ def _ortho_pixelize( mask_uint8, chunk[ptype, px_name].to("code_length"), chunk[ptype, py_name].to("code_length"), + chunk[ptype, pz_name].to("code_length"), chunk[ptype, "smoothing_length"].to("code_length"), chunk[ptype, "mass"].to("code_mass"), chunk[ptype, "density"].to("code_density"), chunk[field].in_units(ounits), - bnds, - check_period=int(periodic), - period=period, + bnds3, + _check_period=_periodic.astype("int"), + period=period3, + kernel_name=kernel_name, ) # We use code length here, but to get the path length right # we need to multiply by the conversion factor between @@ -430,14 +469,16 @@ def _ortho_pixelize( mask_uint8, chunk[ptype, px_name].to("code_length"), chunk[ptype, py_name].to("code_length"), + chunk[ptype, pz_name].to("code_length"), chunk[ptype, "smoothing_length"].to("code_length"), chunk[ptype, "mass"].to("code_mass"), chunk[ptype, "density"].to("code_density"), chunk[field].in_units(ounits), - bnds, - check_period=int(periodic), - period=period, + bnds3, + _check_period=_periodic.astype("int"), + period=period3, weight_field=chunk[weight].in_units(wounits), + kernel_name=kernel_name, ) mylog.info( "Making a fixed resolution buffer of (%s) %d by %d", @@ -452,13 +493,15 @@ def _ortho_pixelize( mask_uint8, chunk[ptype, px_name].to("code_length"), chunk[ptype, py_name].to("code_length"), + chunk[ptype, pz_name].to("code_length"), chunk[ptype, "smoothing_length"].to("code_length"), chunk[ptype, "mass"].to("code_mass"), chunk[ptype, "density"].to("code_density"), chunk[weight].in_units(wounits), - bnds, - check_period=int(periodic), - period=period, + bnds3, + _check_period=_periodic.astype("int"), + period=period3, + kernel_name=kernel_name, ) normalization_2d_utility(buff, weight_buff) if moment == 2: @@ -471,14 +514,16 @@ def _ortho_pixelize( mask_uint8, chunk[ptype, px_name].to("code_length"), chunk[ptype, py_name].to("code_length"), + chunk[ptype, pz_name].to("code_length"), chunk[ptype, "smoothing_length"].to("code_length"), chunk[ptype, "mass"].to("code_mass"), chunk[ptype, "density"].to("code_density"), chunk[field].in_units(ounits) ** 2, - bnds, - check_period=int(periodic), - period=period, + bnds3, + _check_period=_periodic.astype("int"), + period=period3, weight_field=chunk[weight].in_units(wounits), + kernel_name=kernel_name, ) normalization_2d_utility(buff2, weight_buff) buff = compute_stddev_image(buff2, buff) @@ -494,32 +539,39 @@ def _ortho_pixelize( buff_den = np.zeros(size, dtype="float64") for chunk in data_source.chunks([], "io"): + hsmlname = "smoothing_length" pixelize_sph_kernel_slice( buff, mask_uint8, - chunk[ptype, px_name].to("code_length"), - chunk[ptype, py_name].to("code_length"), - chunk[ptype, "smoothing_length"].to("code_length"), - chunk[ptype, "mass"].to("code_mass"), - chunk[ptype, "density"].to("code_density"), - chunk[field].in_units(ounits), + chunk[ptype, px_name].to("code_length").v, + chunk[ptype, py_name].to("code_length").v, + chunk[ptype, pz_name].to("code_length").v, + chunk[ptype, hsmlname].to("code_length").v, + chunk[ptype, "mass"].to("code_mass").v, + chunk[ptype, "density"].to("code_density").v, + chunk[field].in_units(ounits).v, bnds, - check_period=int(periodic), - period=period, + data_source.coord.to("code_length").v, + _check_period=_periodic.astype("int"), + period=period3, + kernel_name=kernel_name, ) if normalize: pixelize_sph_kernel_slice( buff_den, mask_uint8, - chunk[ptype, px_name].to("code_length"), - chunk[ptype, py_name].to("code_length"), - chunk[ptype, "smoothing_length"].to("code_length"), - chunk[ptype, "mass"].to("code_mass"), - chunk[ptype, "density"].to("code_density"), + chunk[ptype, px_name].to("code_length").v, + chunk[ptype, py_name].to("code_length").v, + chunk[ptype, pz_name].to("code_length").v, + chunk[ptype, hsmlname].to("code_length").v, + chunk[ptype, "mass"].to("code_mass").v, + chunk[ptype, "density"].to("code_density").v, np.ones(chunk[ptype, "density"].shape[0]), bnds, - check_period=int(periodic), - period=period, + data_source.coord.to("code_length").v, + _check_period=_periodic.astype("int"), + period=period3, + kernel_name=kernel_name, ) if normalize: @@ -608,7 +660,7 @@ def _ortho_pixelize( data_source[field], bounds, int(antialias), - period, + period2, int(periodic), return_mask=True, ) @@ -616,29 +668,130 @@ def _ortho_pixelize( return buff, mask def _oblique_pixelize(self, data_source, field, bounds, size, antialias): + from yt.data_objects.selection_objects.slices import YTCuttingPlane + from yt.frontends.sph.data_structures import ParticleDataset + from yt.frontends.stream.data_structures import StreamParticlesDataset from yt.frontends.ytdata.data_structures import YTSpatialPlotDataset - indices = np.argsort(data_source["pdx"])[::-1].astype("int64", copy=False) - buff = np.full((size[1], size[0]), np.nan, dtype="float64") - ftype = "index" - if isinstance(data_source.ds, YTSpatialPlotDataset): - ftype = "gas" - mask = pixelize_off_axis_cartesian( - buff, - data_source[ftype, "x"], - data_source[ftype, "y"], - data_source[ftype, "z"], - data_source["px"], - data_source["py"], - data_source["pdx"], - data_source["pdy"], - data_source["pdz"], - data_source.center, - data_source._inv_mat, - indices, - data_source[field], - bounds, - ) + # Determine what sort of data we're dealing with + # -> what backend to use + # copied from the _ortho_pixelize method + field = data_source._determine_fields(field)[0] + _finfo = data_source.ds.field_info[field] + is_sph_field = _finfo.is_sph_field + particle_datasets = (ParticleDataset, StreamParticlesDataset) + # finfo = self.ds._get_field_info(field) + + # SPH data + # only for slices: a function in off_axis_projection.py + # handles projections + if ( + isinstance(data_source.ds, particle_datasets) + and is_sph_field + and isinstance(data_source, YTCuttingPlane) + ): + normalize = getattr(self.ds, "use_sph_normalization", True) + le = data_source.ds.domain_left_edge.to("code_length") + re = data_source.ds.domain_right_edge.to("code_length") + boxbounds = np.array([le[0], re[0], le[1], re[1], le[2], re[2]]) + periodic = data_source.ds.periodicity + ptype = field[0] + if ptype == "gas": + ptype = data_source.ds._sph_ptypes[0] + axorder = data_source.ds.coordinates.axis_order + ounits = data_source.ds.field_info[field].output_units + # input bounds are in code length units already + widthxy = np.array((bounds[1] - bounds[0], bounds[3] - bounds[2])) + kernel_name = None + if hasattr(data_source.ds, "kernel_name"): + kernel_name = data_source.ds.kernel_name + if kernel_name is None: + kernel_name = "cubic" + # data_source should be a YTCuttingPlane object + # dimensionless unyt normal/north + # -> numpy array cython can deal with + normal_vector = data_source.normal.v + north_vector = data_source._y_vec.v + center = data_source.center.to("code_length") + + buff = np.zeros(size, dtype="float64") + mask_uint8 = np.zeros_like(buff, dtype="uint8") + if normalize: + buff_den = np.zeros(size, dtype="float64") + + for chunk in data_source.chunks([], "io"): + pixelize_sph_kernel_cutting( + buff, + mask_uint8, + chunk[ptype, axorder[0]].to("code_length").v, + chunk[ptype, axorder[1]].to("code_length").v, + chunk[ptype, axorder[2]].to("code_length").v, + chunk[ptype, "smoothing_length"].to("code_length").v, + chunk[ptype, "mass"].to("code_mass"), + chunk[ptype, "density"].to("code_density"), + chunk[field].in_units(ounits), + center, + widthxy, + normal_vector, + north_vector, + boxbounds, + periodic, + kernel_name=kernel_name, + check_period=1, + ) + if normalize: + pixelize_sph_kernel_cutting( + buff_den, + mask_uint8, + chunk[ptype, axorder[0]].to("code_length"), + chunk[ptype, axorder[1]].to("code_length"), + chunk[ptype, axorder[2]].to("code_length"), + chunk[ptype, "smoothing_length"].to("code_length"), + chunk[ptype, "mass"].to("code_mass"), + chunk[ptype, "density"].to("code_density"), + np.ones(chunk[ptype, "density"].shape[0]), + center, + widthxy, + normal_vector, + north_vector, + boxbounds, + periodic, + kernel_name=kernel_name, + check_period=1, + ) + + if normalize: + normalization_2d_utility(buff, buff_den) + + mask = mask_uint8.astype("bool", copy=False) + # swap axes for image plotting + mask = mask.swapaxes(0, 1) + buff = buff.swapaxes(0, 1) + + # whatever other data this code could handle before the + # SPH option was added + else: + indices = np.argsort(data_source["pdx"])[::-1].astype("int64", copy=False) + buff = np.full((size[1], size[0]), np.nan, dtype="float64") + ftype = "index" + if isinstance(data_source.ds, YTSpatialPlotDataset): + ftype = "gas" + mask = pixelize_off_axis_cartesian( + buff, + data_source[ftype, "x"], + data_source[ftype, "y"], + data_source[ftype, "z"], + data_source["px"], + data_source["py"], + data_source["pdx"], + data_source["pdy"], + data_source["pdz"], + data_source.center, + data_source._inv_mat, + indices, + data_source[field], + bounds, + ) return buff, mask def convert_from_cartesian(self, coord): diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization.py b/yt/geometry/coordinates/tests/test_sph_pixelization.py index 6de8b29f77..1d1c82be5c 100644 --- a/yt/geometry/coordinates/tests/test_sph_pixelization.py +++ b/yt/geometry/coordinates/tests/test_sph_pixelization.py @@ -1,7 +1,19 @@ +import numpy as np + import yt -from yt.testing import assert_rel_equal, requires_file +from yt.testing import ( + assert_rel_equal, + cubicspline_python, + fake_sph_flexible_grid_ds, + integrate_kernel, + requires_file, +) from yt.utilities.math_utils import compute_stddev_image +## off-axis projection tests for SPH data are in +## yt/visualization/tests/test_offaxisprojection.py + + magneticum = "MagneticumCluster/snap_132" mag_kwargs = { @@ -36,3 +48,172 @@ def _vysq(field, data): ) sigy = compute_stddev_image(prj1.frb["gas", "vysq"], prj1.frb["gas", "velocity_y"]) assert_rel_equal(sigy, prj2.frb["gas", "velocity_y"].d, 10) + + +def test_sph_projection_basic1(): + """ + small, uniform grid: expected values for given dl? + pixel centers at 0.5, 1., 1.5, 2., 2.5 + particles at 0.5, 1.5, 2.5 + """ + bbox = np.array([[0.0, 3.0]] * 3) + ds = fake_sph_flexible_grid_ds(hsml_factor=1.0, nperside=3, bbox=bbox) + # works, but no depth control (at least without specific filters) + proj = ds.proj(("gas", "density"), 2) + frb = proj.to_frb( + width=(2.5, "cm"), + resolution=(5, 5), + height=(2.5, "cm"), + center=np.array([1.5, 1.5, 1.5]), + periodic=False, + ) + out = frb.get_image(("gas", "density")) + + expected_out = np.zeros((5, 5), dtype=np.float64) + dl_1part = integrate_kernel(cubicspline_python, 0.0, 0.5) + linedens_1part = dl_1part * 1.0 # unit mass, density + linedens = 3.0 * linedens_1part + expected_out[::2, ::2] = linedens + + assert_rel_equal(expected_out, out.v, 5) + # return out + + +def test_sph_projection_basic2(): + """ + small, uniform grid: expected values for given dl? + pixel centers at 0.5, 1., 1.5, 2., 2.5 + particles at 0.5, 1.5, 2.5 + but hsml radii are 0.25 -> try non-zero impact parameters, + other pixels are still zero. + """ + bbox = np.array([[0.0, 3.0]] * 3) + ds = fake_sph_flexible_grid_ds(hsml_factor=0.5, nperside=3, bbox=bbox) + proj = ds.proj(("gas", "density"), 2) + frb = proj.to_frb( + width=(2.5, "cm"), + resolution=(5, 5), + height=(2.5, "cm"), + center=np.array([1.375, 1.375, 1.5]), + periodic=False, + ) + out = frb.get_image(("gas", "density")) + + expected_out = np.zeros((5, 5), dtype=np.float64) + dl_1part = integrate_kernel(cubicspline_python, np.sqrt(2) * 0.125, 0.25) + linedens_1part = dl_1part * 1.0 # unit mass, density + linedens = 3.0 * linedens_1part + expected_out[::2, ::2] = linedens + + # print(expected_out) + # print(out.v) + assert_rel_equal(expected_out, out.v, 4) + # return out + + +def get_dataset_sphrefine(reflevel: int = 1): + """ + constant density particle grid, + with increasing particle sampling + """ + lenfact = (1.0 / 3.0) ** (reflevel - 1) + massfact = lenfact**3 + nperside = 3**reflevel + + e1hat = np.array([lenfact, 0, 0]) + e2hat = np.array([0, lenfact, 0]) + e3hat = np.array([0, 0, lenfact]) + hsml_factor = lenfact + bbox = np.array([[0.0, 3.0]] * 3) + offsets = np.ones(3, dtype=np.float64) * 0.5 # in units of ehat + + def refmass(i: int, j: int, k: int) -> float: + return massfact + + unitrho = 1.0 / massfact # want density 1 for decreasing mass + + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=nperside, + periodic=True, + e1hat=e1hat, + e2hat=e2hat, + e3hat=e3hat, + offsets=offsets, + massgenerator=refmass, + unitrho=unitrho, + bbox=bbox, + ) + return ds + + +def getdata_test_gridproj2(): + # initial pixel centers at 0.5, 1., 1.5, 2., 2.5 + # particles at 0.5, 1.5, 2.5 + # refine particle grid, check if pixel values remain the + # same in the pixels passing through initial particle centers + outlist = [] + dss = [] + for rl in range(1, 4): + ds = get_dataset_sphrefine(reflevel=rl) + proj = ds.proj(("gas", "density"), 2) + frb = proj.to_frb( + width=(2.5, "cm"), + resolution=(5, 5), + height=(2.5, "cm"), + center=np.array([1.5, 1.5, 1.5]), + periodic=False, + ) + out = frb.get_image(("gas", "density")) + outlist.append(out) + dss.append(ds) + return outlist, dss + + +def test_sph_gridproj_reseffect1(): + """ + Comparing same pixel centers with higher particle resolution. + The pixel centers are at x/y coordinates [0.5, 1., 1.5, 2., 2.5] + at the first level, the spacing halves at each level. + Checking the pixels at [0.5, 1.5, 2.5], + which should have the same values at each resolution. + """ + imgs, _ = getdata_test_gridproj2() + ref = imgs[-1] + for img in imgs: + assert_rel_equal( + img[:: img.shape[0] // 2, :: img.shape[1] // 2], + ref[:: ref.shape[0] // 2, :: ref.shape[1] // 2], + 4, + ) + + +def test_sph_gridproj_reseffect2(): + """ + refine the pixel grid instead of the particle grid + """ + ds = get_dataset_sphrefine(reflevel=2) + proj = ds.proj(("gas", "density"), 2) + imgs = {} + maxrl = 5 + for rl in range(1, maxrl + 1): + npix = 1 + 2 ** (rl + 1) + margin = 0.5 - 0.5 ** (rl + 1) + frb = proj.to_frb( + width=(3.0 - 2.0 * margin, "cm"), + resolution=(npix, npix), + height=(3.0 - 2.0 * margin, "cm"), + center=np.array([1.5, 1.5, 1.5]), + periodic=False, + ) + out = frb.get_image(("gas", "density")) + imgs[rl] = out + ref = imgs[maxrl] + pixspace_ref = 2 ** (maxrl) + for rl in imgs: + img = imgs[rl] + pixspace = 2 ** (rl) + # print(f'Grid refinement level {rl}:') + assert_rel_equal( + img[::pixspace, ::pixspace], ref[::pixspace_ref, ::pixspace_ref], 4 + ) diff --git a/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py new file mode 100644 index 0000000000..75c7175e06 --- /dev/null +++ b/yt/geometry/coordinates/tests/test_sph_pixelization_pytestonly.py @@ -0,0 +1,547 @@ +from typing import Union + +import numpy as np +import pytest +import unyt + +import yt +from yt.data_objects.selection_objects.region import YTRegion +from yt.testing import ( + assert_rel_equal, + cubicspline_python, + distancematrix, + fake_random_sph_ds, + fake_sph_flexible_grid_ds, + integrate_kernel, +) + + +@pytest.mark.parametrize("weighted", [True, False]) +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm")]) +@pytest.mark.parametrize("shiftcenter", [False, True]) +@pytest.mark.parametrize("axis", [0, 1, 2]) +def test_sph_proj_general_alongaxes( + axis: int, + shiftcenter: bool, + depth: Union[float, None], + periodic: bool, + weighted: bool, +) -> None: + """ + The previous projection tests were for a specific issue. + Here, we test more functionality of the projections. + We just send lines of sight through pixel centers for convenience. + Particles at [0.5, 1.5, 2.5] (in each coordinate) + smoothing lengths 0.25 + all particles have mass 1., density 1.5, + except the single center particle, with mass 2., density 3. + + Parameters: + ----------- + axis: {0, 1, 2} + projection axis (aligned with sim. axis) + shiftcenter: bool + shift the coordinates to center the projection on. + (The grid is offset to this same center) + depth: float or None + depth of the projection slice + periodic: bool + assume periodic boundary conditions, or not + weighted: bool + make a weighted projection (density-weighted density), or not + + Returns: + -------- + None + """ + if shiftcenter: + center = unyt.unyt_array(np.array((0.625, 0.625, 0.625)), "cm") + else: + center = unyt.unyt_array(np.array((1.5, 1.5, 1.5)), "cm") + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") + hsml_factor = 0.5 + unitrho = 1.5 + + # test correct centering, particle selection + def makemasses(i, j, k): + if i == j == k == 1: + return 2.0 + else: + return 1.0 + + # m / rho, factor 1. / hsml**2 is included in the kernel integral + # (density is adjusted, so same for center particle) + prefactor = 1.0 / unitrho # / (0.5 * 0.5)**2 + dl_cen = integrate_kernel(cubicspline_python, 0.0, 0.25) + + # result shouldn't depend explicitly on the center if we re-center + # the data, unless we get cut-offs in the non-periodic case + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=3, + periodic=periodic, + offsets=np.full(3, 0.5), + massgenerator=makemasses, + unitrho=unitrho, + bbox=bbox.v, + recenter=center.v, + ) + if depth is None: + source = ds.all_data() + else: + depth = unyt.unyt_quantity(*depth) + le = np.array(ds.domain_left_edge) + re = np.array(ds.domain_right_edge) + le[axis] = center[axis] - 0.5 * depth + re[axis] = center[axis] + 0.5 * depth + cen = 0.5 * (le + re) + reg = YTRegion(center=cen, left_edge=le, right_edge=re, ds=ds) + source = reg + + # we don't actually want a plot, it's just a straightforward, + # common way to get an frb / image array + if weighted: + toweight_field = ("gas", "density") + else: + toweight_field = None + prj = yt.ProjectionPlot( + ds, + axis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=toweight_field, + buff_size=(5, 5), + center=center, + data_source=source, + ) + img = prj.frb.data[("gas", "density")] + if weighted: + expected_out = np.zeros( + ( + 5, + 5, + ), + dtype=img.v.dtype, + ) + expected_out[::2, ::2] = unitrho + if depth is None: + ## during shift, particle coords do wrap around edges + # if (not periodic) and shiftcenter: + # # weight 1. for unitrho, 2. for 2. * untrho + # expected_out[2, 2] *= 5. / 3. + # else: + # weight (2 * 1.) for unitrho, (1 * 2.) for 2. * unitrho + expected_out[2, 2] *= 1.5 + else: + # only 2 * unitrho element included + expected_out[2, 2] *= 2.0 + else: + expected_out = np.zeros( + ( + 5, + 5, + ), + dtype=img.v.dtype, + ) + expected_out[::2, ::2] = dl_cen * prefactor * unitrho + if depth is None: + # 3 particles per l.o.s., including the denser one + expected_out *= 3.0 + expected_out[2, 2] *= 4.0 / 3.0 + else: + # 1 particle per l.o.s., including the denser one + expected_out[2, 2] *= 2.0 + # grid is shifted to the left -> 'missing' stuff at the left + if (not periodic) and shiftcenter: + expected_out[:1, :] = 0.0 + expected_out[:, :1] = 0.0 + # print(axis, shiftcenter, depth, periodic, weighted) + # print(expected_out) + # print(img.v) + assert_rel_equal(expected_out, img.v, 5) + + +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("shiftcenter", [False, True]) +@pytest.mark.parametrize("zoff", [0.0, 0.1, 0.5, 1.0]) +@pytest.mark.parametrize("axis", [0, 1, 2]) +def test_sph_slice_general_alongaxes( + axis: int, + shiftcenter: bool, + periodic: bool, + zoff: float, +) -> None: + """ + Particles at [0.5, 1.5, 2.5] (in each coordinate) + smoothing lengths 0.25 + all particles have mass 1., density 1.5, + except the single center particle, with mass 2., density 3. + + Parameters: + ----------- + axis: {0, 1, 2} + projection axis (aligned with sim. axis) + northvector: tuple + y-axis direction in the final plot (direction vector) + shiftcenter: bool + shift the coordinates to center the projection on. + (The grid is offset to this same center) + zoff: float + offset of the slice plane from the SPH particle center plane + periodic: bool + assume periodic boundary conditions, or not + + Returns: + -------- + None + """ + if shiftcenter: + center = unyt.unyt_array(np.array((0.625, 0.625, 0.625)), "cm") + else: + center = unyt.unyt_array(np.array((1.5, 1.5, 1.5)), "cm") + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") + hsml_factor = 0.5 + unitrho = 1.5 + + # test correct centering, particle selection + def makemasses(i, j, k): + if i == j == k == 1: + return 2.0 + elif i == j == k == 2: + return 3.0 + else: + return 1.0 + + # result shouldn't depend explicitly on the center if we re-center + # the data, unless we get cut-offs in the non-periodic case + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=3, + periodic=periodic, + offsets=np.full(3, 0.5), + massgenerator=makemasses, + unitrho=unitrho, + bbox=bbox.v, + recenter=center.v, + ) + ad = ds.all_data() + # print(ad[('gas', 'position')]) + outgridsize = 10 + width = 2.5 + _center = center.to("cm").v.copy() + _center[axis] += zoff + + # we don't actually want a plot, it's just a straightforward, + # common way to get an frb / image array + slc = yt.SlicePlot( + ds, + axis, + ("gas", "density"), + width=(width, "cm"), + buff_size=(outgridsize,) * 2, + center=(_center, "cm"), + ) + img = slc.frb.data[("gas", "density")] + + # center is same in non-projection coords + if axis == 0: + ci = 1 + else: + ci = 0 + gridcens = ( + _center[ci] + - 0.5 * width + + 0.5 * width / outgridsize + + np.arange(outgridsize) * width / outgridsize + ) + xgrid = np.repeat(gridcens, outgridsize) + ygrid = np.tile(gridcens, outgridsize) + zgrid = np.full(outgridsize**2, _center[axis]) + gridcoords = np.empty((outgridsize**2, 3), dtype=xgrid.dtype) + if axis == 2: + gridcoords[:, 0] = xgrid + gridcoords[:, 1] = ygrid + gridcoords[:, 2] = zgrid + elif axis == 0: + gridcoords[:, 0] = zgrid + gridcoords[:, 1] = xgrid + gridcoords[:, 2] = ygrid + elif axis == 1: + gridcoords[:, 0] = ygrid + gridcoords[:, 1] = zgrid + gridcoords[:, 2] = xgrid + ad = ds.all_data() + sphcoords = np.array( + [ + (ad[("gas", "x")]).to("cm"), + (ad[("gas", "y")]).to("cm"), + (ad[("gas", "z")]).to("cm"), + ] + ).T + # print("sphcoords:") + # print(sphcoords) + # print("gridcoords:") + # print(gridcoords) + dists = distancematrix( + gridcoords, + sphcoords, + periodic=(periodic,) * 3, + periods=np.array([3.0, 3.0, 3.0]), + ) + # print("dists <= 1:") + # print(dists <= 1) + sml = (ad[("gas", "smoothing_length")]).to("cm") + normkern = cubicspline_python(dists / sml.v[np.newaxis, :]) + sphcontr = normkern / sml[np.newaxis, :] ** 3 * ad[("gas", "mass")] + contsum = np.sum(sphcontr, axis=1) + sphweights = ( + normkern + / sml[np.newaxis, :] ** 3 + * ad[("gas", "mass")] + / ad[("gas", "density")] + ) + weights = np.sum(sphweights, axis=1) + nzeromask = np.logical_not(weights == 0) + expected = np.zeros(weights.shape, weights.dtype) + expected[nzeromask] = contsum[nzeromask] / weights[nzeromask] + expected = expected.reshape((outgridsize, outgridsize)) + # expected[np.isnan(expected)] = 0.0 # convention in the slices + + # print("expected:\n", expected) + # print("recovered:\n", img.v) + assert_rel_equal(expected, img.v, 5) + + +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("shiftcenter", [False, True]) +@pytest.mark.parametrize("northvector", [None, (1.0e-4, 1.0, 0.0)]) +@pytest.mark.parametrize("zoff", [0.0, 0.1, 0.5, 1.0]) +def test_sph_slice_general_offaxis( + northvector: Union[tuple[float, float, float], None], + shiftcenter: bool, + zoff: float, + periodic: bool, +) -> None: + """ + Same as the on-axis slices, but we rotate the basis vectors + to test whether roations are handled ok. the rotation is chosen to + be small so that in/exclusion of particles within bboxes, etc. + works out the same way. + Particles at [0.5, 1.5, 2.5] (in each coordinate) + smoothing lengths 0.25 + all particles have mass 1., density 1.5, + except the single center particle, with mass 2., density 3. + + Parameters: + ----------- + northvector: tuple + y-axis direction in the final plot (direction vector) + shiftcenter: bool + shift the coordinates to center the projection on. + (The grid is offset to this same center) + zoff: float + offset of the slice plane from the SPH particle center plane + periodic: bool + assume periodic boundary conditions, or not + + Returns: + -------- + None + """ + if shiftcenter: + center = np.array((0.625, 0.625, 0.625)) # cm + else: + center = np.array((1.5, 1.5, 1.5)) # cm + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") + hsml_factor = 0.5 + unitrho = 1.5 + + # test correct centering, particle selection + def makemasses(i, j, k): + if i == j == k == 1: + return 2.0 + else: + return 1.0 + + # try to make sure dl differences from periodic wrapping are small + epsilon = 1e-4 + projaxis = np.array([epsilon, 0.00, np.sqrt(1.0 - epsilon**2)]) + e1dir = projaxis / np.sqrt(np.sum(projaxis**2)) + if northvector is None: + e2dir = np.array([0.0, 1.0, 0.0]) + else: + e2dir = np.asarray(northvector) + e2dir = e2dir - np.sum(e1dir * e2dir) * e2dir # orthonormalize + e2dir /= np.sqrt(np.sum(e2dir**2)) + e3dir = np.cross(e2dir, e1dir) + + outgridsize = 10 + width = 2.5 + _center = center.copy() + _center += zoff * e1dir + + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=3, + periodic=periodic, + offsets=np.full(3, 0.5), + massgenerator=makemasses, + unitrho=unitrho, + bbox=bbox.v, + recenter=center, + e1hat=e1dir, + e2hat=e2dir, + e3hat=e3dir, + ) + + # source = ds.all_data() + # couple to dataset -> right unit registry + center = ds.arr(center, "cm") + # print("position:\n", source["gas", "position"]) + slc = yt.SlicePlot( + ds, + e1dir, + ("gas", "density"), + width=(width, "cm"), + buff_size=(outgridsize,) * 2, + center=(_center, "cm"), + north_vector=e2dir, + ) + img = slc.frb.data[("gas", "density")] + + # center is same in x/y (e3dir/e2dir) + gridcenx = ( + np.dot(_center, e3dir) + - 0.5 * width + + 0.5 * width / outgridsize + + np.arange(outgridsize) * width / outgridsize + ) + gridceny = ( + np.dot(_center, e2dir) + - 0.5 * width + + 0.5 * width / outgridsize + + np.arange(outgridsize) * width / outgridsize + ) + xgrid = np.repeat(gridcenx, outgridsize) + ygrid = np.tile(gridceny, outgridsize) + zgrid = np.full(outgridsize**2, np.dot(_center, e1dir)) + gridcoords = ( + xgrid[:, np.newaxis] * e3dir[np.newaxis, :] + + ygrid[:, np.newaxis] * e2dir[np.newaxis, :] + + zgrid[:, np.newaxis] * e1dir[np.newaxis, :] + ) + # print("gridcoords:") + # print(gridcoords) + ad = ds.all_data() + sphcoords = np.array( + [ + (ad[("gas", "x")]).to("cm"), + (ad[("gas", "y")]).to("cm"), + (ad[("gas", "z")]).to("cm"), + ] + ).T + dists = distancematrix( + gridcoords, + sphcoords, + periodic=(periodic,) * 3, + periods=np.array([3.0, 3.0, 3.0]), + ) + sml = (ad[("gas", "smoothing_length")]).to("cm") + normkern = cubicspline_python(dists / sml.v[np.newaxis, :]) + sphcontr = normkern / sml[np.newaxis, :] ** 3 * ad[("gas", "mass")] + contsum = np.sum(sphcontr, axis=1) + sphweights = ( + normkern + / sml[np.newaxis, :] ** 3 + * ad[("gas", "mass")] + / ad[("gas", "density")] + ) + weights = np.sum(sphweights, axis=1) + nzeromask = np.logical_not(weights == 0) + expected = np.zeros(weights.shape, weights.dtype) + expected[nzeromask] = contsum[nzeromask] / weights[nzeromask] + expected = expected.reshape((outgridsize, outgridsize)) + expected = expected.T # transposed for image plotting + # expected[np.isnan(expected)] = 0.0 # convention in the slices + + # print(axis, shiftcenter, depth, periodic, weighted) + # print("expected:\n", expected) + # print("recovered:\n", img.v) + assert_rel_equal(expected, img.v, 4) + + +# only axis-aligned; testing YTArbitraryGrid, YTCoveringGrid +@pytest.mark.parametrize("periodic", [True, False, (True, True, False)]) +@pytest.mark.parametrize("wholebox", [True, False]) +def test_sph_grid( + periodic: Union[bool, tuple[bool, bool, bool]], + wholebox: bool, +) -> None: + bbox = np.array([[-1.0, 3.0], [1.0, 5.2], [-1.0, 3.0]]) + ds = fake_random_sph_ds(50, bbox, periodic=periodic) + + if not hasattr(periodic, "__len__"): + periodic = (periodic,) * 3 + + if wholebox: + left = bbox[:, 0].copy() + level = 2 + ncells = np.array([2**level] * 3) + # print("left: ", left) + # print("ncells: ", ncells) + resgrid = ds.covering_grid(level, tuple(left), ncells) + right = bbox[:, 1].copy() + xedges = np.linspace(left[0], right[0], ncells[0] + 1) + yedges = np.linspace(left[1], right[1], ncells[1] + 1) + zedges = np.linspace(left[2], right[2], ncells[2] + 1) + else: + left = np.array([-1.0, 1.8, -1.0]) + right = np.array([2.5, 5.2, 2.5]) + ncells = np.array([3, 4, 4]) + resgrid = ds.arbitrary_grid(left, right, dims=ncells) + xedges = np.linspace(left[0], right[0], ncells[0] + 1) + yedges = np.linspace(left[1], right[1], ncells[1] + 1) + zedges = np.linspace(left[2], right[2], ncells[2] + 1) + res = resgrid["gas", "density"] + xcens = 0.5 * (xedges[:-1] + xedges[1:]) + ycens = 0.5 * (yedges[:-1] + yedges[1:]) + zcens = 0.5 * (zedges[:-1] + zedges[1:]) + + ad = ds.all_data() + sphcoords = np.array( + [ + (ad[("gas", "x")]).to("cm"), + (ad[("gas", "y")]).to("cm"), + (ad[("gas", "z")]).to("cm"), + ] + ).T + gridx, gridy, gridz = np.meshgrid(xcens, ycens, zcens, indexing="ij") + outshape = gridx.shape + gridx = gridx.flatten() + gridy = gridy.flatten() + gridz = gridz.flatten() + gridcoords = np.array([gridx, gridy, gridz]).T + periods = bbox[:, 1] - bbox[:, 0] + dists = distancematrix(gridcoords, sphcoords, periodic=periodic, periods=periods) + sml = (ad[("gas", "smoothing_length")]).to("cm") + normkern = cubicspline_python(dists / sml.v[np.newaxis, :]) + sphcontr = normkern / sml[np.newaxis, :] ** 3 * ad[("gas", "mass")] + contsum = np.sum(sphcontr, axis=1) + sphweights = ( + normkern + / sml[np.newaxis, :] ** 3 + * ad[("gas", "mass")] + / ad[("gas", "density")] + ) + weights = np.sum(sphweights, axis=1) + nzeromask = np.logical_not(weights == 0) + expected = np.zeros(weights.shape, weights.dtype) + expected[nzeromask] = contsum[nzeromask] / weights[nzeromask] + expected = expected.reshape(outshape) + # expected[np.isnan(expected)] = 0.0 # convention in the slices + + # print(axis, shiftcenter, depth, periodic, weighted) + # print("expected:\n", expected) + # print("recovered:\n", res.v) + assert_rel_equal(expected, res.v, 4) diff --git a/yt/loaders.py b/yt/loaders.py index 54d52e5b83..8494649549 100644 --- a/yt/loaders.py +++ b/yt/loaders.py @@ -9,6 +9,7 @@ import time import types import warnings +from collections.abc import Mapping from pathlib import Path from typing import TYPE_CHECKING, Any, Optional, Union, cast from urllib.parse import urlsplit @@ -694,7 +695,7 @@ def load_amr_grids( def load_particles( - data: dict[AnyFieldKey, np.ndarray], + data: Mapping[AnyFieldKey, Union[np.ndarray, tuple[np.ndarray, str]]], length_unit=None, bbox=None, sim_time=None, @@ -833,7 +834,7 @@ def parse_unit(unit, dimension): field_units, data, _ = process_data(data) sfh = StreamDictFieldHandler() - pdata: dict[AnyFieldKey, np.ndarray] = {} + pdata: dict[AnyFieldKey, Union[np.ndarray, tuple[np.ndarray, str]]] = {} for key in data.keys(): field: FieldKey if not isinstance(key, tuple): diff --git a/yt/testing.py b/yt/testing.py index 2205439ea1..a5e6f10d83 100644 --- a/yt/testing.py +++ b/yt/testing.py @@ -6,9 +6,11 @@ import sys import tempfile import unittest +from collections.abc import Mapping from functools import wraps from importlib.util import find_spec from shutil import which +from typing import Callable, Union from unittest import SkipTest import matplotlib @@ -18,9 +20,11 @@ from unyt.exceptions import UnitOperationError from yt._maintenance.deprecation import issue_deprecation_warning +from yt._typing import AnyFieldKey from yt.config import ytcfg +from yt.frontends.stream.data_structures import StreamParticlesDataset from yt.funcs import is_sequence -from yt.loaders import load +from yt.loaders import load, load_particles from yt.units.yt_array import YTArray, YTQuantity ANSWER_TEST_TAG = "answer_test" @@ -80,6 +84,111 @@ def assert_rel_equal(a1, a2, decimals, err_msg="", verbose=True): ) +# tested: volume integral is 1. +def cubicspline_python( + x: Union[float, np.ndarray], +) -> np.ndarray: + """ + cubic spline SPH kernel function for testing against more + effiecient cython methods + + Parameters + ---------- + x: + impact parameter / smoothing length [dimenionless] + + Returns + ------- + value of the kernel function + """ + # C is 8/pi + _c = 8.0 / np.pi + x = np.asarray(x) + kernel = np.zeros(x.shape, dtype=x.dtype) + half1 = np.where(np.logical_and(x >= 0.0, x <= 0.5)) + kernel[half1] = 1.0 - 6.0 * x[half1] ** 2 * (1.0 - x[half1]) + half2 = np.where(np.logical_and(x > 0.5, x <= 1.0)) + kernel[half2] = 2.0 * (1.0 - x[half2]) ** 3 + return kernel * _c + + +def integrate_kernel( + kernelfunc: Callable[[float], float], b: float, hsml: float +) -> float: + """ + integrates a kernel function over a line passing entirely + through it + + Parameters: + ----------- + kernelfunc: + the kernel function to integrate + b: + impact parameter + hsml: + smoothing length [same units as impact parameter] + + Returns: + -------- + the integral of the SPH kernel function. + units: 1 / units of b and hsml + """ + pre = 1.0 / hsml**2 + x = b / hsml + xmax = np.sqrt(1.0 - x**2) + xmin = -1.0 * xmax + xe = np.linspace(xmin, xmax, 500) # shape: 500, x.shape + xc = 0.5 * (xe[:-1, ...] + xe[1:, ...]) + dx = np.diff(xe, axis=0) + spv = kernelfunc(np.sqrt(xc**2 + x**2)) + integral = np.sum(spv * dx, axis=0) + return pre * integral + + +_zeroperiods = np.array([0.0, 0.0, 0.0]) + + +def distancematrix( + pos3_i0: np.ndarray, + pos3_i1: np.ndarray, + periodic: tuple[bool, bool, bool] = (True,) * 3, + periods: np.ndarray = _zeroperiods, +) -> np.ndarray: + """ + Calculates the distances between two arrays of points. + + Parameters: + ---------- + pos3_i0: shape (first number of points, 3) + positions of the first set of points. The second index is + for positions along the different cartesian axes + pos3_i1: shape (second number of points, 3) + as pos3_i0, but for the second set of points + periodic: + are the positions along each axis periodic (True) or not + periods: + the periods along each axis. Ignored if positions in a given + direction are not periodic. + + Returns: + -------- + a 2D-array of distances between postions `pos3_i0` (changes along + index 0) and `pos3_i1` (changes along index 1) + + """ + d2 = np.zeros((len(pos3_i0), len(pos3_i1)), dtype=pos3_i0.dtype) + for ax in range(3): + # 'center on' pos3_i1 + _d = pos3_i0[:, ax, np.newaxis] - pos3_i1[np.newaxis, :, ax] + if periodic[ax]: + _period = periods[ax] + _d += 0.5 * _period # center on half box size + _d %= _period # wrap coordinate to 0 -- boxsize range + _d -= 0.5 * _period # center back to zero + d2 += _d**2 + return np.sqrt(d2) + + def amrspace(extent, levels=7, cells=8): """Creates two numpy arrays representing the left and right bounds of an AMR grid as well as an array for the AMR level of each cell. @@ -686,6 +795,220 @@ def fake_sph_grid_ds(hsml_factor=1.0): return load_particles(data=data, length_unit=1.0, bbox=bbox) +def constantmass(i: int, j: int, k: int) -> float: + return 1.0 + + +_xhat = np.array([1, 0, 0]) +_yhat = np.array([0, 1, 0]) +_zhat = np.array([0, 0, 1]) +_floathalves = 0.5 * np.ones((3,), dtype=np.float64) + + +def fake_sph_flexible_grid_ds( + hsml_factor: float = 1.0, + nperside: int = 3, + periodic: bool = True, + e1hat: np.ndarray = _xhat, + e2hat: np.ndarray = _yhat, + e3hat: np.ndarray = _zhat, + offsets: np.ndarray = _floathalves, + massgenerator: Callable[[int, int, int], float] = constantmass, + unitrho: float = 1.0, + bbox: Union[np.ndarray, None] = None, + recenter: Union[np.ndarray, None] = None, +) -> StreamParticlesDataset: + """Returns an in-memory SPH dataset useful for testing + + Parameters: + ----------- + hsml_factor: + all particles have smoothing lengths of `hsml_factor` * 0.5 + nperside: + the dataset will have `nperside`**3 particles, arranged + uniformly on a 3D grid + periodic: + are the positions taken to be periodic? (applies to all + coordinate axes) + e1hat: shape (3,) + the first basis vector defining the 3D grid. If the basis + vectors are not normalized to 1 or not orthogonal, the spacing + or overlap between SPH particles will be affected, but this is + allowed. + e2hat: shape (3,) + the second basis vector defining the 3D grid. (See `e1hat`.) + e3hat: shape (3,) + the third basis vector defining the 3D grid. (See `e1hat`.) + offsets: shape (3,) + the the zero point of the 3D grid along each coordinate axis + massgenerator: + a function assigning a mass to each particle, as a function of + the e[1-3]hat indices, in order + unitrho: + defines the density for a particle with mass 1 ('g'), and the + standard (uniform) grid `hsml_factor`. + bbox: if np.ndarray, shape is (2, 3) + the assumed enclosing volume of the particles. Should enclose + all the coordinate values. If not specified, a bbox is defined + which encloses all coordinates values with a margin. If + `periodic`, the size of the `bbox` along each coordinate is + also the period along that axis. + recenter: + if not `None`, after generating the grid, the positions are + periodically shifted to move the old center to this positions. + Useful for testing periodicity handling. + This shift is relative to the halfway positions of the bbox + edges. + + Returns: + -------- + A `StreamParticlesDataset` object with particle positions, masses, + velocities (zero), smoothing lengths, and densities specified. + Values are in cgs units. + """ + + npart = nperside**3 + + pos = np.empty((npart, 3), dtype=np.float64) + mass = np.empty((npart,), dtype=np.float64) + for i in range(0, nperside): + for j in range(0, nperside): + for k in range(0, nperside): + _pos = ( + (offsets[0] + i) * e1hat + + (offsets[1] + j) * e2hat + + (offsets[2] + k) * e3hat + ) + ind = nperside**2 * i + nperside * j + k + pos[ind, :] = _pos + mass[ind] = massgenerator(i, j, k) + rho = unitrho * mass + + if bbox is None: + eps = 1e-3 + margin = (1.0 + eps) * hsml_factor + bbox = np.array( + [ + [np.min(pos[:, 0]) - margin, np.max(pos[:, 0]) + margin], + [np.min(pos[:, 1]) - margin, np.max(pos[:, 1]) + margin], + [np.min(pos[:, 2]) - margin, np.max(pos[:, 2]) + margin], + ] + ) + + if recenter is not None: + periods = bbox[:, 1] - bbox[:, 0] + # old center -> new position + pos += -0.5 * periods[np.newaxis, :] + recenter[np.newaxis, :] + # wrap coordinates -> all in [0, boxsize) range + pos %= periods[np.newaxis, :] + # shift back to original bbox range + pos += (bbox[:, 0])[np.newaxis, :] + if not periodic: + # remove points outside bbox to avoid errors: + okinds = np.ones(len(mass), dtype=bool) + for ax in [0, 1, 2]: + okinds &= pos[:, ax] < bbox[ax, 1] + okinds &= pos[:, ax] >= bbox[ax, 0] + npart = sum(okinds) + else: + okinds = np.ones((npart,), dtype=bool) + + data: Mapping[AnyFieldKey, tuple[np.ndarray, str]] = { + "particle_position_x": (np.copy(pos[okinds, 0]), "cm"), + "particle_position_y": (np.copy(pos[okinds, 1]), "cm"), + "particle_position_z": (np.copy(pos[okinds, 2]), "cm"), + "particle_mass": (np.copy(mass[okinds]), "g"), + "particle_velocity_x": (np.zeros(npart), "cm/s"), + "particle_velocity_y": (np.zeros(npart), "cm/s"), + "particle_velocity_z": (np.zeros(npart), "cm/s"), + "smoothing_length": (np.ones(npart) * 0.5 * hsml_factor, "cm"), + "density": (np.copy(rho[okinds]), "g/cm**3"), + } + + ds = load_particles( + data=data, + bbox=bbox, + periodicity=(periodic,) * 3, + length_unit=1.0, + mass_unit=1.0, + time_unit=1.0, + velocity_unit=1.0, + ) + ds.kernel_name = "cubic" + return ds + + +def fake_random_sph_ds( + npart: int, + bbox: np.ndarray, + periodic: Union[bool, tuple[bool, bool, bool]] = True, + massrange: tuple[float, float] = (0.5, 2.0), + hsmlrange: tuple[float, float] = (0.5, 2.0), + unitrho: float = 1.0, +) -> StreamParticlesDataset: + """Returns an in-memory SPH dataset useful for testing + + Parameters: + ----------- + npart: + number of particles to generate + bbox: shape: (3, 2), units: "cm" + the assumed enclosing volume of the particles. Particle + positions are drawn uniformly from these ranges. + periodic: + are the positions taken to be periodic? If a single value, + that value is applied to all axes + massrange: + particle masses are drawn uniformly from this range (unit: "g") + hsmlrange: units: "cm" + particle smoothing lengths are drawn uniformly from this range + unitrho: + defines the density for a particle with mass 1 ("g"), and + smoothing length 1 ("cm"). + + Returns: + -------- + A `StreamParticlesDataset` object with particle positions, masses, + velocities (zero), smoothing lengths, and densities specified. + Values are in cgs units. + """ + + if not hasattr(periodic, "__len__"): + periodic = (periodic,) * 3 + gen = np.random.default_rng(seed=0) + + posx = gen.uniform(low=bbox[0][0], high=bbox[0][1], size=npart) + posy = gen.uniform(low=bbox[1][0], high=bbox[1][1], size=npart) + posz = gen.uniform(low=bbox[2][0], high=bbox[2][1], size=npart) + mass = gen.uniform(low=massrange[0], high=massrange[1], size=npart) + hsml = gen.uniform(low=hsmlrange[0], high=hsmlrange[1], size=npart) + dens = mass / hsml**3 * unitrho + + data: Mapping[AnyFieldKey, tuple[np.ndarray, str]] = { + "particle_position_x": (posx, "cm"), + "particle_position_y": (posy, "cm"), + "particle_position_z": (posz, "cm"), + "particle_mass": (mass, "g"), + "particle_velocity_x": (np.zeros(npart), "cm/s"), + "particle_velocity_y": (np.zeros(npart), "cm/s"), + "particle_velocity_z": (np.zeros(npart), "cm/s"), + "smoothing_length": (hsml, "cm"), + "density": (dens, "g/cm**3"), + } + + ds = load_particles( + data=data, + bbox=bbox, + periodicity=periodic, + length_unit=1.0, + mass_unit=1.0, + time_unit=1.0, + velocity_unit=1.0, + ) + ds.kernel_name = "cubic" + return ds + + def construct_octree_mask(prng=RandomState(0x1D3D3D3), refined=None): # noqa B008 # Implementation taken from url: # http://docs.hyperion-rt.org/en/stable/advanced/indepth_oct.html diff --git a/yt/utilities/answer_testing/framework.py b/yt/utilities/answer_testing/framework.py index 37358573dd..ad2db01c4f 100644 --- a/yt/utilities/answer_testing/framework.py +++ b/yt/utilities/answer_testing/framework.py @@ -776,9 +776,9 @@ def compare(self, new_result, old_result): def dump_images(new_result, old_result, decimals=10): - tmpfd, old_image = tempfile.mkstemp(suffix=".png") + tmpfd, old_image = tempfile.mkstemp(prefix="baseline_", suffix=".png") os.close(tmpfd) - tmpfd, new_image = tempfile.mkstemp(suffix=".png") + tmpfd, new_image = tempfile.mkstemp(prefix="thisPR_", suffix=".png") os.close(tmpfd) image_writer.write_projection(new_result, new_image) image_writer.write_projection(old_result, old_image) diff --git a/yt/utilities/lib/pixelization_routines.pyx b/yt/utilities/lib/pixelization_routines.pyx index 8cc37ab1d1..a207473e9e 100644 --- a/yt/utilities/lib/pixelization_routines.pyx +++ b/yt/utilities/lib/pixelization_routines.pyx @@ -1123,6 +1123,7 @@ def pixelize_sph_kernel_projection( np.uint8_t[:, :] mask, any_float[:] posx, any_float[:] posy, + any_float[:] posz, any_float[:] hsml, any_float[:] pmass, any_float[:] pdens, @@ -1130,21 +1131,24 @@ def pixelize_sph_kernel_projection( bounds, kernel_name="cubic", weight_field=None, - int check_period=1, + _check_period = (1, 1, 1), period=None): cdef np.intp_t xsize, ysize - cdef np.float64_t x_min, x_max, y_min, y_max, prefactor_j + cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max, prefactor_j cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi cdef np.float64_t q_ij2, posx_diff, posy_diff, ih_j2 - cdef np.float64_t x, y, dx, dy, idx, idy, h_j2, px, py - cdef np.float64_t period_x = 0, period_y = 0 - cdef int i, j, ii, jj + cdef np.float64_t x, y, dx, dy, idx, idy, h_j2, px, py, pz + cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 + cdef int i, j, ii, jj, kk cdef np.float64_t[:] _weight_field cdef int * xiter cdef int * yiter + cdef int * ziter cdef np.float64_t * xiterv cdef np.float64_t * yiterv + cdef np.float64_t * ziterv + cdef np.int8_t[3] check_period if weight_field is not None: _weight_field = weight_field @@ -1152,7 +1156,9 @@ def pixelize_sph_kernel_projection( if period is not None: period_x = period[0] period_y = period[1] - + period_z = period[2] + for i in range(3): + check_period[i] = np.int8(_check_period[i]) # we find the x and y range over which we have pixels and we find how many # pixels we have in each dimension xsize, ysize = buff.shape[0], buff.shape[1] @@ -1160,6 +1166,8 @@ def pixelize_sph_kernel_projection( x_max = bounds[1] y_min = bounds[2] y_max = bounds[3] + z_min = bounds[4] + z_max = bounds[5] dx = (x_max - x_min) / xsize dy = (y_max - y_min) / ysize @@ -1170,7 +1178,6 @@ def pixelize_sph_kernel_projection( if kernel_name not in kernel_tables: kernel_tables[kernel_name] = SPHKernelInterpolationTable(kernel_name) cdef SPHKernelInterpolationTable itab = kernel_tables[kernel_name] - with nogil, parallel(): # loop through every particle # NOTE: this loop can be quite time consuming. However it is easily @@ -1190,10 +1197,12 @@ def pixelize_sph_kernel_projection( local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) xiterv = malloc(sizeof(np.float64_t) * 2) yiterv = malloc(sizeof(np.float64_t) * 2) + ziterv = malloc(sizeof(np.float64_t) * 2) xiter = malloc(sizeof(int) * 2) yiter = malloc(sizeof(int) * 2) - xiter[0] = yiter[0] = 0 - xiterv[0] = yiterv[0] = 0.0 + ziter = malloc(sizeof(int) * 2) + xiter[0] = yiter[0] = ziter[0] = 0 + xiterv[0] = yiterv[0] = ziterv[0] = 0.0 for i in range(xsize * ysize): local_buff[i] = 0.0 @@ -1202,75 +1211,98 @@ def pixelize_sph_kernel_projection( with gil: PyErr_CheckSignals() - xiter[1] = yiter[1] = 999 + xiter[1] = yiter[1] = ziter[1] = 999 - if check_period == 1: + if check_period[0] == 1: if posx[j] - hsml[j] < x_min: xiter[1] = +1 xiterv[1] = period_x elif posx[j] + hsml[j] > x_max: xiter[1] = -1 xiterv[1] = -period_x + if check_period[1] == 1: if posy[j] - hsml[j] < y_min: yiter[1] = +1 yiterv[1] = period_y elif posy[j] + hsml[j] > y_max: yiter[1] = -1 yiterv[1] = -period_y + if check_period[2] == 1: + if posz[j] - hsml[j] < z_min: + ziter[1] = +1 + ziterv[1] = period_z + elif posz[j] + hsml[j] > z_max: + ziter[1] = -1 + ziterv[1] = -period_z # we set the smoothing length squared with lower limit of the pixel - h_j2 = fmax(hsml[j]*hsml[j], dx*dy) + # Nope! that causes weird grid resolution dependences and increases + # total values when resolution elements have hsml < grid spacing + h_j2 = hsml[j]*hsml[j] ih_j2 = 1.0/h_j2 prefactor_j = pmass[j] / pdens[j] / hsml[j]**2 * quantity_to_smooth[j] if weight_field is not None: prefactor_j *= _weight_field[j] - for ii in range(2): - if xiter[ii] == 999: continue - px = posx[j] + xiterv[ii] - if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): continue - for jj in range(2): - if yiter[jj] == 999: continue - py = posy[j] + yiterv[jj] - if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): continue - - # here we find the pixels which this particle contributes to - x0 = ((px - hsml[j] - x_min)*idx) - x1 = ((px + hsml[j] - x_min)*idx) - x0 = iclip(x0-1, 0, xsize) - x1 = iclip(x1+1, 0, xsize) + # Discussion point: do we want the hsml margin on the z direction? + # it's consistent with Ray and Region selections, I think, + # but does tend to 'tack on' stuff compared to the nominal depth + for kk in range(2): + # discard if z is outside bounds + if ziter[kk] == 999: continue + pz = posz[j] + ziterv[kk] + ## removed hsml 'margin' in the projection direction to avoid + ## double-counting particles near periodic edges + ## and adding extra 'depth' to projections + #if (pz + hsml[j] < z_min) or (pz - hsml[j] > z_max): continue + if (pz < z_min) or (pz > z_max): continue + + for ii in range(2): + if xiter[ii] == 999: continue + px = posx[j] + xiterv[ii] + if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): continue + for jj in range(2): + if yiter[jj] == 999: continue + py = posy[j] + yiterv[jj] + if (py + hsml[j] < y_min) or (py - hsml[j] > y_max): continue + + # here we find the pixels which this particle contributes to + x0 = ((px - hsml[j] - x_min)*idx) + x1 = ((px + hsml[j] - x_min)*idx) + x0 = iclip(x0-1, 0, xsize) + x1 = iclip(x1+1, 0, xsize) - y0 = ((py - hsml[j] - y_min)*idy) - y1 = ((py + hsml[j] - y_min)*idy) - y0 = iclip(y0-1, 0, ysize) - y1 = iclip(y1+1, 0, ysize) + y0 = ((py - hsml[j] - y_min)*idy) + y1 = ((py + hsml[j] - y_min)*idy) + y0 = iclip(y0-1, 0, ysize) + y1 = iclip(y1+1, 0, ysize) - # found pixels we deposit on, loop through those pixels - for xi in range(x0, x1): - # we use the centre of the pixel to calculate contribution - x = (xi + 0.5) * dx + x_min + # found pixels we deposit on, loop through those pixels + for xi in range(x0, x1): + # we use the centre of the pixel to calculate contribution + x = (xi + 0.5) * dx + x_min - posx_diff = px - x - posx_diff = posx_diff * posx_diff + posx_diff = px - x + posx_diff = posx_diff * posx_diff - if posx_diff > h_j2: continue + if posx_diff > h_j2: continue - for yi in range(y0, y1): - y = (yi + 0.5) * dy + y_min + for yi in range(y0, y1): + y = (yi + 0.5) * dy + y_min - posy_diff = py - y - posy_diff = posy_diff * posy_diff - if posy_diff > h_j2: continue + posy_diff = py - y + posy_diff = posy_diff * posy_diff + if posy_diff > h_j2: continue - q_ij2 = (posx_diff + posy_diff) * ih_j2 - if q_ij2 >= 1: - continue + q_ij2 = (posx_diff + posy_diff) * ih_j2 + if q_ij2 >= 1: + continue - # see equation 32 of the SPLASH paper - # now we just use the kernel projection - local_buff[xi + yi*xsize] += prefactor_j * itab.interpolate(q_ij2) - mask[xi, yi] = 1 + # see equation 32 of the SPLASH paper + # now we just use the kernel projection + local_buff[xi + yi*xsize] += prefactor_j * itab.interpolate(q_ij2) + mask[xi, yi] = 1 with gil: for xxi in range(xsize): @@ -1466,32 +1498,46 @@ def pixelize_sph_kernel_slice( np.float64_t[:, :] buff, np.uint8_t[:, :] mask, np.float64_t[:] posx, np.float64_t[:] posy, + np.float64_t[:] posz, np.float64_t[:] hsml, np.float64_t[:] pmass, np.float64_t[:] pdens, np.float64_t[:] quantity_to_smooth, - bounds, kernel_name="cubic", - int check_period=1, + bounds, + np.float64_t slicez, + kernel_name="cubic", + _check_period = (1, 1, 1), period=None): - + #print("bounds, slicez, kernel_name, check_period, period") + #print(bounds) + #print(slicez) + #print(kernel_name) + #print(check_period) + #print(period) + #print() + # bounds are [x0, x1, y0, y1], slicez is the single coordinate + # of the slice along the normal direction. # similar method to pixelize_sph_kernel_projection cdef np.intp_t xsize, ysize cdef np.float64_t x_min, x_max, y_min, y_max, prefactor_j cdef np.int64_t xi, yi, x0, x1, y0, y1, xxi, yyi - cdef np.float64_t q_ij, posx_diff, posy_diff, ih_j - cdef np.float64_t x, y, dx, dy, idx, idy, h_j2, h_j, px, py + cdef np.float64_t q_ij, posx_diff, posy_diff, posz_diff, ih_j + cdef np.float64_t x, y, dx, dy, idx, idy, h_j2, h_j, px, py, pz cdef int i, j, ii, jj - cdef np.float64_t period_x = 0, period_y = 0 + cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 cdef int * xiter cdef int * yiter cdef np.float64_t * xiterv cdef np.float64_t * yiterv + cdef np.int8_t[3] check_period if period is not None: period_x = period[0] period_y = period[1] + period_z = period[2] + for i in range(3): + check_period[i] = np.int8(_check_period[i]) xsize, ysize = buff.shape[0], buff.shape[1] - x_min = bounds[0] x_max = bounds[1] y_min = bounds[2] @@ -1503,7 +1549,7 @@ def pixelize_sph_kernel_slice( idy = 1.0/dy kernel = get_kernel_func(kernel_name) - + #print('particle index, ii, jj, px, py, pz') with nogil, parallel(): # NOTE see note in pixelize_sph_kernel_projection local_buff = malloc(sizeof(np.float64_t) * xsize * ysize) @@ -1520,27 +1566,41 @@ def pixelize_sph_kernel_slice( if j % 100000 == 0: with gil: PyErr_CheckSignals() - + #with gil: + # print(j) xiter[1] = yiter[1] = 999 - - if check_period == 1: + pz = posz[j] + if check_period[0] == 1: if posx[j] - hsml[j] < x_min: - xiter[1] = +1 + xiter[1] = 1 xiterv[1] = period_x elif posx[j] + hsml[j] > x_max: xiter[1] = -1 xiterv[1] = -period_x + if check_period[1] == 1: if posy[j] - hsml[j] < y_min: - yiter[1] = +1 + yiter[1] = 1 yiterv[1] = period_y elif posy[j] + hsml[j] > y_max: yiter[1] = -1 yiterv[1] = -period_y - - h_j2 = fmax(hsml[j]*hsml[j], dx*dy) - h_j = math.sqrt(h_j2) + if check_period[2] == 1: + # z of particle might be < hsml from the slice plane + # but across a periodic boundary + if posz[j] - hsml[j] > slicez: + pz = posz[j] - period_z + elif posz[j] + hsml[j] < slicez: + pz = posz[j] + period_z + + h_j2 = hsml[j] * hsml[j] #fmax(hsml[j]*hsml[j], dx*dy) + h_j = hsml[j] #math.sqrt(h_j2) ih_j = 1.0/h_j + posz_diff = pz - slicez + posz_diff = posz_diff * posz_diff + if posz_diff > h_j2: + continue + prefactor_j = pmass[j] / pdens[j] / hsml[j]**3 prefactor_j *= quantity_to_smooth[j] @@ -1562,7 +1622,8 @@ def pixelize_sph_kernel_slice( y1 = ( (py + hsml[j] - y_min) * idy) y0 = iclip(y0-1, 0, ysize) y1 = iclip(y1+1, 0, ysize) - + #with gil: + # print(ii, jj, px, py, pz) # Now we know which pixels to deposit onto for this particle, # so loop over them and add this particle's contribution for xi in range(x0, x1): @@ -1582,7 +1643,9 @@ def pixelize_sph_kernel_slice( continue # see equation 4 of the SPLASH paper - q_ij = math.sqrt(posx_diff + posy_diff) * ih_j + q_ij = math.sqrt(posx_diff + + posy_diff + + posz_diff) * ih_j if q_ij >= 1: continue @@ -1610,13 +1673,14 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff, np.float64_t[:] pdens, np.float64_t[:] quantity_to_smooth, bounds, pbar=None, kernel_name="cubic", - int check_period=1, period=None): + check_period=True, period=None): cdef np.intp_t xsize, ysize, zsize cdef np.float64_t x_min, x_max, y_min, y_max, z_min, z_max, prefactor_j cdef np.int64_t xi, yi, zi, x0, x1, y0, y1, z0, z1 cdef np.float64_t q_ij, posx_diff, posy_diff, posz_diff, px, py, pz - cdef np.float64_t x, y, z, dx, dy, dz, idx, idy, idz, h_j3, h_j2, h_j, ih_j + cdef np.float64_t x, y, z, dx, dy, dz, idx, idy, idz, h_j2, h_j, ih_j + # cdef np.float64_t h_j3 cdef int j, ii, jj, kk cdef np.float64_t period_x = 0, period_y = 0, period_z = 0 @@ -1626,10 +1690,21 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff, cdef np.float64_t xiterv[2] cdef np.float64_t yiterv[2] cdef np.float64_t ziterv[2] + cdef int[3] periodic + xiter[0] = yiter[0] = ziter[0] = 0 xiterv[0] = yiterv[0] = ziterv[0] = 0.0 + if hasattr(check_period, "__len__"): + periodic[0] = int(check_period[0]) + periodic[1] = int(check_period[1]) + periodic[2] = int(check_period[2]) + else: + _cp = int(check_period) + periodic[0] = _cp + periodic[1] = _cp + periodic[2] = _cp if period is not None: period_x = period[0] period_y = period[1] @@ -1652,6 +1727,16 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff, kernel = get_kernel_func(kernel_name) + # nogil seems dangerous here, but there are no actual parallel + # sections (e.g., prange instead of range) used here. + # However, for future writers: + # !! the final buff array mutation has no protections against + # !! race conditions (e.g., OpenMP's atomic read/write), and + # !! cython doesn't seem to provide such options. + # (other routines in this file use private variable buffer arrays + # and add everything together at the end, but grid arrays can get + # big fast, and having such a large array in each thread could + # cause memory use issues.) with nogil: # TODO make this parallel without using too much memory for j in range(0, posx.shape[0]): @@ -1660,23 +1745,26 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff, if(pbar is not None): pbar.update(50000) PyErr_CheckSignals() + # end with gil xiter[1] = yiter[1] = ziter[1] = 999 xiterv[1] = yiterv[1] = ziterv[1] = 0.0 - if check_period == 1: + if periodic[0] == 1: if posx[j] - hsml[j] < x_min: xiter[1] = +1 xiterv[1] = period_x elif posx[j] + hsml[j] > x_max: xiter[1] = -1 xiterv[1] = -period_x + if periodic[1] == 1: if posy[j] - hsml[j] < y_min: yiter[1] = +1 yiterv[1] = period_y elif posy[j] + hsml[j] > y_max: yiter[1] = -1 yiterv[1] = -period_y + if periodic[2] == 1: if posz[j] - hsml[j] < z_min: ziter[1] = +1 ziterv[1] = period_z @@ -1684,8 +1772,8 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff, ziter[1] = -1 ziterv[1] = -period_z - h_j3 = fmax(hsml[j]*hsml[j]*hsml[j], dx*dy*dz) - h_j = math.cbrt(h_j3) + #h_j3 = fmax(hsml[j]*hsml[j]*hsml[j], dx*dy*dz) + h_j = hsml[j] #math.cbrt(h_j3) h_j2 = h_j*h_j ih_j = 1/h_j @@ -1746,11 +1834,17 @@ def pixelize_sph_kernel_arbitrary_grid(np.float64_t[:, :, :] buff, continue # see equation 4 of the SPLASH paper - q_ij = math.sqrt(posx_diff + posy_diff + posz_diff) * ih_j + q_ij = math.sqrt(posx_diff + + posy_diff + + posz_diff) * ih_j if q_ij >= 1: continue - - buff[xi, yi, zi] += prefactor_j * kernel(q_ij) + # shared variable buff should not + # be mutatated in a nogil section + # where different threads may change + # the same array element + buff[xi, yi, zi] += prefactor_j \ + * kernel(q_ij) def pixelize_element_mesh_line(np.ndarray[np.float64_t, ndim=2] coords, @@ -1857,16 +1951,16 @@ def pixelize_element_mesh_line(np.ndarray[np.float64_t, ndim=2] coords, free(field_vals) return arc_length, plot_values - +# intended for use in ParticleImageBuffer @cython.boundscheck(False) @cython.wraparound(False) -def rotate_particle_coord(np.float64_t[:] px, - np.float64_t[:] py, - np.float64_t[:] pz, - center, - width, - normal_vector, - north_vector): +def rotate_particle_coord_pib(np.float64_t[:] px, + np.float64_t[:] py, + np.float64_t[:] pz, + center, + width, + normal_vector, + north_vector): # We want to do two rotations, one to first rotate our coordinates to have # the normal vector be the z-axis (i.e., the viewer's perspective), and then # another rotation to make the north-vector be the y-axis (i.e., north). @@ -1909,6 +2003,89 @@ def rotate_particle_coord(np.float64_t[:] px, return px_rotated, py_rotated, rot_bounds_x0, rot_bounds_x1, rot_bounds_y0, rot_bounds_y1 +# version intended for SPH off-axis slices/projections +# includes dealing with periodic boundaries, but also +# shifts particles so center -> origin. +# therefore, don't want to use this in the ParticleImageBuffer, +# which expects differently centered coordinates. +@cython.boundscheck(False) +@cython.wraparound(False) +def rotate_particle_coord(np.float64_t[:] px, + np.float64_t[:] py, + np.float64_t[:] pz, + center, + bounds, + periodic, + width, + depth, + normal_vector, + north_vector): + # We want to do two rotations, one to first rotate our coordinates to have + # the normal vector be the z-axis (i.e., the viewer's perspective), and then + # another rotation to make the north-vector be the y-axis (i.e., north). + # Fortunately, total_rotation_matrix = rotation_matrix_1 x rotation_matrix_2 + cdef np.int64_t num_particles = np.size(px) + cdef np.float64_t[:] z_axis = np.array([0., 0., 1.], dtype="float64") + cdef np.float64_t[:] y_axis = np.array([0., 1., 0.], dtype="float64") + cdef np.float64_t[:, :] normal_rotation_matrix + cdef np.float64_t[:] transformed_north_vector + cdef np.float64_t[:, :] north_rotation_matrix + cdef np.float64_t[:, :] rotation_matrix + + normal_rotation_matrix = get_rotation_matrix(normal_vector, z_axis) + transformed_north_vector = np.matmul(normal_rotation_matrix, north_vector) + north_rotation_matrix = get_rotation_matrix(transformed_north_vector, y_axis) + rotation_matrix = np.matmul(north_rotation_matrix, normal_rotation_matrix) + + cdef np.float64_t[:] px_rotated = np.empty(num_particles, dtype="float64") + cdef np.float64_t[:] py_rotated = np.empty(num_particles, dtype="float64") + cdef np.float64_t[:] pz_rotated = np.empty(num_particles, dtype="float64") + cdef np.float64_t[:] coordinate_matrix = np.empty(3, dtype="float64") + cdef np.float64_t[:] rotated_coordinates + cdef np.float64_t[:] rotated_center + cdef np.int64_t i + cdef int ax + #rotated_center = rotation_matmul( + # rotation_matrix, np.array([center[0], center[1], center[2]])) + rotated_center = np.zeros((3,), dtype=center.dtype) + # set up the rotated bounds + cdef np.float64_t rot_bounds_x0 = rotated_center[0] - 0.5 * width[0] + cdef np.float64_t rot_bounds_x1 = rotated_center[0] + 0.5 * width[0] + cdef np.float64_t rot_bounds_y0 = rotated_center[1] - 0.5 * width[1] + cdef np.float64_t rot_bounds_y1 = rotated_center[1] + 0.5 * width[1] + cdef np.float64_t rot_bounds_z0 = rotated_center[2] - 0.5 * depth + cdef np.float64_t rot_bounds_z1 = rotated_center[2] + 0.5 * depth + for i in range(num_particles): + coordinate_matrix[0] = px[i] + coordinate_matrix[1] = py[i] + coordinate_matrix[2] = pz[i] + + # centering: + # make sure this also works for centers close to periodic edges + # added consequence: the center is placed at the origin + # (might as well keep it there in these temporary coordinates) + for ax in range(3): + # assumed center is zero even if non-periodic + coordinate_matrix[ax] -= center[ax] + if not periodic[ax]: continue + period = bounds[2 * ax + 1] - bounds[2 * ax] + # abs. difference between points in the volume is <= period + if coordinate_matrix[ax] < -0.5 * period: + coordinate_matrix[ax] += period + if coordinate_matrix[ax] > 0.5 * period: + coordinate_matrix[ax] -= period + + rotated_coordinates = rotation_matmul( + rotation_matrix, coordinate_matrix) + px_rotated[i] = rotated_coordinates[0] + py_rotated[i] = rotated_coordinates[1] + pz_rotated[i] = rotated_coordinates[2] + + return (px_rotated, py_rotated, pz_rotated, + rot_bounds_x0, rot_bounds_x1, + rot_bounds_y0, rot_bounds_y1, + rot_bounds_z0, rot_bounds_z1) + @cython.boundscheck(False) @cython.wraparound(False) @@ -1921,33 +2098,95 @@ def off_axis_projection_SPH(np.float64_t[:] px, bounds, center, width, + periodic, np.float64_t[:] quantity_to_smooth, np.float64_t[:, :] projection_array, np.uint8_t[:, :] mask, normal_vector, north_vector, - weight_field=None): + weight_field=None, + depth=None, + kernel_name="cubic"): + # periodic: periodicity of the data set: # Do nothing in event of a 0 normal vector if np.allclose(normal_vector, 0.): return - - px_rotated, py_rotated, \ + if depth is None: + # set to volume diagonal + margin -> won't exclude anything + depth = 2. * np.sqrt((bounds[1] - bounds[0])**2 + + (bounds[3] - bounds[2])**2 + + (bounds[5] - bounds[4])**2) + px_rotated, py_rotated, pz_rotated, \ rot_bounds_x0, rot_bounds_x1, \ - rot_bounds_y0, rot_bounds_y1 = rotate_particle_coord(px, py, pz, - center, width, normal_vector, north_vector) - + rot_bounds_y0, rot_bounds_y1, \ + rot_bounds_z0, rot_bounds_z1 = rotate_particle_coord(px, py, pz, + center, bounds, + periodic, + width, depth, + normal_vector, + north_vector) + # check_period=0: assumed to be a small region compared to the box + # size. The rotation already ensures that a center close to a + # periodic edge works out fine. + # since the simple single-coordinate modulo math periodicity + # does not apply to the *rotated* coordinates, the periodicity + # approach implemented for this along-axis projection method + # would fail here + check_period = np.array([0, 0, 0], dtype="int") pixelize_sph_kernel_projection(projection_array, mask, px_rotated, py_rotated, + pz_rotated, smoothing_lengths, particle_masses, particle_densities, quantity_to_smooth, [rot_bounds_x0, rot_bounds_x1, - rot_bounds_y0, rot_bounds_y1], + rot_bounds_y0, rot_bounds_y1, + rot_bounds_z0, rot_bounds_z1], weight_field=weight_field, - check_period=0) + _check_period=check_period, + kernel_name=kernel_name) + +# like slice pixelization, but for off-axis planes +def pixelize_sph_kernel_cutting( + np.float64_t[:, :] buff, + np.uint8_t[:, :] mask, + np.float64_t[:] posx, np.float64_t[:] posy, + np.float64_t[:] posz, + np.float64_t[:] hsml, np.float64_t[:] pmass, + np.float64_t[:] pdens, + np.float64_t[:] quantity_to_smooth, + center, widthxy, + normal_vector, north_vector, + boxbounds, periodic, + kernel_name="cubic", + int check_period=1): + + if check_period == 0: + periodic = np.zeros(3, dtype=bool) + + posx_rot, posy_rot, posz_rot, \ + rot_bounds_x0, rot_bounds_x1, \ + rot_bounds_y0, rot_bounds_y1, \ + rot_bounds_z0, _ = rotate_particle_coord(posx, posy, posz, + center, boxbounds, + periodic, + widthxy, 0., + normal_vector, + north_vector) + bounds_rot = np.array([rot_bounds_x0, rot_bounds_x1, + rot_bounds_y0, rot_bounds_y1]) + slicez_rot = rot_bounds_z0 + pixelize_sph_kernel_slice(buff, mask, + posx_rot, posy_rot, posz_rot, + hsml, pmass, pdens, quantity_to_smooth, + bounds_rot, slicez_rot, + kernel_name=kernel_name, + _check_period=np.zeros(3, dtype="int"), + period=None) + @cython.boundscheck(False) diff --git a/yt/visualization/fits_image.py b/yt/visualization/fits_image.py index 85a098fa30..c740d91287 100644 --- a/yt/visualization/fits_image.py +++ b/yt/visualization/fits_image.py @@ -1648,6 +1648,7 @@ def __init__( north_vector=north_vector, method=method, weight=weight_field, + depth=depth, ).swapaxes(0, 1) if moment == 2: @@ -1675,6 +1676,7 @@ def _sq_field(field, data, item: FieldKey): north_vector=north_vector, method=method, weight=weight_field, + depth=depth, ).swapaxes(0, 1) buf[key] = compute_stddev_image(buff2, buf[key]) diff --git a/yt/visualization/fixed_resolution.py b/yt/visualization/fixed_resolution.py index d2da6a29e1..000c2dba8b 100644 --- a/yt/visualization/fixed_resolution.py +++ b/yt/visualization/fixed_resolution.py @@ -17,7 +17,7 @@ ) from yt.utilities.lib.pixelization_routines import ( pixelize_cylinder, - rotate_particle_coord, + rotate_particle_coord_pib, ) from yt.utilities.math_utils import compute_stddev_image from yt.utilities.on_demand_imports import _h5py as h5py @@ -630,6 +630,8 @@ def _generate_image_and_mask(self, item) -> None: self.buff_size[1], ) dd = self.data_source + # only need the first two for SPH, + # but need the third one for other data formats. width = self.ds.arr( ( self.bounds[1] - self.bounds[0], @@ -637,6 +639,7 @@ def _generate_image_and_mask(self, item) -> None: self.bounds[5] - self.bounds[4], ) ) + depth = dd.depth[0] if dd.depth is not None else None buff = off_axis_projection( dd.dd, dd.center, @@ -649,6 +652,7 @@ def _generate_image_and_mask(self, item) -> None: no_ghost=dd.no_ghost, interpolated=dd.interpolated, north_vector=dd.north_vector, + depth=depth, method=dd.method, ) if self.data_source.moment == 2: @@ -679,6 +683,7 @@ def _sq_field(field, data, item: FieldKey): no_ghost=dd.no_ghost, interpolated=dd.interpolated, north_vector=dd.north_vector, + depth=dd.depth, method=dd.method, ) buff = compute_stddev_image(buff2, buff) @@ -745,7 +750,7 @@ def _generate_image_and_mask(self, item) -> None: if hasattr(w, "to_value"): w = w.to_value("code_length") wd.append(w) - x_data, y_data, *bounds = rotate_particle_coord( + x_data, y_data, *bounds = rotate_particle_coord_pib( dd[ftype, "particle_position_x"].to_value("code_length"), dd[ftype, "particle_position_y"].to_value("code_length"), dd[ftype, "particle_position_z"].to_value("code_length"), diff --git a/yt/visualization/plot_window.py b/yt/visualization/plot_window.py index 4534d7325c..8176065f47 100644 --- a/yt/visualization/plot_window.py +++ b/yt/visualization/plot_window.py @@ -12,6 +12,8 @@ from yt._maintenance.deprecation import issue_deprecation_warning from yt._typing import AlphaT from yt.data_objects.image_array import ImageArray +from yt.frontends.sph.data_structures import ParticleDataset +from yt.frontends.stream.data_structures import StreamParticlesDataset from yt.frontends.ytdata.data_structures import YTSpatialPlotDataset from yt.funcs import ( fix_axis, @@ -20,6 +22,7 @@ iter_fields, mylog, obj_length, + parse_center_array, validate_moment, ) from yt.geometry.api import Geometry @@ -80,20 +83,33 @@ def get_window_parameters(axis, center, width, ds): return (bounds, center, display_center) -def get_oblique_window_parameters(normal, center, width, ds, depth=None): +def get_oblique_window_parameters( + normal, center, width, ds, depth=None, get3bounds=False +): center, display_center = ds.coordinates.sanitize_center(center, axis=None) width = ds.coordinates.sanitize_width(normal, width, depth) if len(width) == 2: # Transforming to the cutting plane coordinate system - center = (center - ds.domain_left_edge) / ds.domain_width - 0.5 + # the original dimensionless center messes up off-axis + # SPH projections though -> don't use this center there + center = ( + (center - ds.domain_left_edge) / ds.domain_width - 0.5 + ) * ds.domain_width (normal, perp1, perp2) = ortho_find(normal) mat = np.transpose(np.column_stack((perp1, perp2, normal))) center = np.dot(mat, center) w = tuple(el.in_units("code_length") for el in width) bounds = tuple(((2 * (i % 2)) - 1) * w[i // 2] / 2 for i in range(len(w) * 2)) - + if get3bounds and depth is None: + # off-axis projection, depth not specified + # -> set 'large enough' depth using half the box diagonal + margin + d2 = ds.domain_width[0].in_units("code_length") ** 2 + d2 += ds.domain_width[1].in_units("code_length") ** 2 + d2 += ds.domain_width[2].in_units("code_length") ** 2 + diag = np.sqrt(d2) + bounds = bounds + (-0.51 * diag, 0.51 * diag) return (bounds, center) @@ -1818,9 +1834,12 @@ def __init__( normal = self.sanitize_normal_vector(ds, normal) # this will handle time series data and controllers axis = fix_axis(normal, ds) + # print('center at SlicePlot init: ', center) + # print('current domain left edge: ', ds.domain_left_edge) (bounds, center, display_center) = get_window_parameters( axis, center, width, ds ) + # print('center after get_window_parameters: ', center) if field_parameters is None: field_parameters = {} @@ -2218,7 +2237,9 @@ def __init__( f"off-axis slices are not supported for {ds.geometry!r} geometry\n" f"currently supported geometries: {self._supported_geometries!r}" ) - + # bounds are in cutting plane coordinates, centered on 0: + # [xmin, xmax, ymin, ymax]. Can derive width/height back + # from these. unit is code_length (bounds, center_rot) = get_oblique_window_parameters(normal, center, width, ds) if field_parameters is None: field_parameters = {} @@ -2273,6 +2294,7 @@ def __init__( le=None, re=None, north_vector=None, + depth=None, method="integrate", data_source=None, *, @@ -2284,6 +2306,7 @@ def __init__( self.axis = None # always true for oblique data objects self.normal_vector = normal_vector self.width = width + self.depth = depth if data_source is None: self.dd = ds.all_data() else: @@ -2421,7 +2444,7 @@ def __init__( fields, center="center", width=None, - depth=(1, "1"), + depth=None, axes_unit=None, weight_field=None, max_level=None, @@ -2439,22 +2462,54 @@ def __init__( ): if ds.geometry not in self._supported_geometries: raise NotImplementedError( - f"off-axis slices are not supported for {ds.geometry!r} geometry\n" - f"currently supported geometries: {self._supported_geometries!r}" + "off-axis slices are not supported" + f" for {ds.geometry!r} geometry\n" + "currently supported geometries:" + f" {self._supported_geometries!r}" ) - + # center_rot normalizes the center to (0,0), + # units match bounds + # for SPH data, we want to input the original center + # the cython backend handles centering to this point and + # rotation. + # get3bounds gets a depth 0.5 * diagonal + margin in the + # depth=None case. (bounds, center_rot) = get_oblique_window_parameters( - normal, center, width, ds, depth=depth + normal, + center, + width, + ds, + depth=depth, + get3bounds=True, ) + # will probably fail if you try to project an SPH and non-SPH + # field in a single call + # checks for SPH fields copied from the + # _ortho_pixelize method in cartesian_coordinates.py + + ## data_source might be None here + ## (OffAxisProjectionDummyDataSource gets used later) + if data_source is None: + data_source = ds.all_data() + field = data_source._determine_fields(fields)[0] + finfo = data_source.ds.field_info[field] + is_sph_field = finfo.is_sph_field + particle_datasets = (ParticleDataset, StreamParticlesDataset) + + if isinstance(data_source.ds, particle_datasets) and is_sph_field: + center_use = parse_center_array(center, ds=data_source.ds, axis=None) + else: + center_use = center_rot fields = list(iter_fields(fields))[:] - oap_width = ds.arr( - (bounds[1] - bounds[0], bounds[3] - bounds[2], bounds[5] - bounds[4]) - ) + # oap_width = ds.arr( + # (bounds[1] - bounds[0], + # bounds[3] - bounds[2]) + # ) OffAxisProj = OffAxisProjectionDummyDataSource( - center_rot, + center_use, ds, normal, - oap_width, + width, fields, interpolated, weight=weight_field, @@ -2463,6 +2518,7 @@ def __init__( le=le, re=re, north_vector=north_vector, + depth=depth, method=method, data_source=data_source, moment=moment, diff --git a/yt/visualization/tests/test_offaxisprojection.py b/yt/visualization/tests/test_offaxisprojection.py index d317e31d8b..739e48d7ea 100644 --- a/yt/visualization/tests/test_offaxisprojection.py +++ b/yt/visualization/tests/test_offaxisprojection.py @@ -13,7 +13,10 @@ fake_octree_ds, fake_random_ds, ) -from yt.visualization.api import OffAxisProjectionPlot, OffAxisSlicePlot +from yt.visualization.api import ( + OffAxisProjectionPlot, + OffAxisSlicePlot, +) from yt.visualization.image_writer import write_projection from yt.visualization.volume_rendering.api import off_axis_projection @@ -235,10 +238,34 @@ def _vlos_sq(field, data): moment=2, buff_size=(400, 400), ) - assert_rel_equal( - np.sqrt( - p1.frb["gas", "velocity_los_squared"] - p1.frb["gas", "velocity_los"] ** 2 - ), - p2.frb["gas", "velocity_los"], - 10, + ## this failed because some - **2 values come out + ## marginally < 0, resulting in unmatched NaN values in the + ## first assert_rel_equal argument. The compute_stddev_image + ## function used in OffAxisProjectionPlot checks for and deals + ## with these cases. + # assert_rel_equal( + # np.sqrt( + # p1.frb["gas", "velocity_los_squared"] - p1.frb["gas", "velocity_los"] ** 2 + # ), + # p2.frb["gas", "velocity_los"], + # 10, + # ) + p1_expsq = p1.frb["gas", "velocity_los_squared"] + p1_sqexp = p1.frb["gas", "velocity_los"] ** 2 + # set values to zero that have **2 - **2 < 0, but + # the absolute values are much smaller than the smallest + # postive values of **2 and **2 + # (i.e., the difference is pretty much zero) + mindiff = 1e-10 * min( + np.min(p1_expsq[p1_expsq > 0]), np.min(p1_sqexp[p1_sqexp > 0]) + ) + # print(mindiff) + safeorbad = np.logical_not( + np.logical_and(p1_expsq - p1_sqexp < 0, p1_expsq - p1_sqexp > -1.0 * mindiff) ) + # avoid errors from sqrt(negative) + # sqrt in zeros_like insures correct units + p1res = np.zeros_like(np.sqrt(p1_expsq)) + p1res[safeorbad] = np.sqrt(p1_expsq[safeorbad] - p1_sqexp[safeorbad]) + p2res = p2.frb["gas", "velocity_los"] + assert_rel_equal(p1res, p2res, 10) diff --git a/yt/visualization/tests/test_offaxisprojection_pytestonly.py b/yt/visualization/tests/test_offaxisprojection_pytestonly.py new file mode 100644 index 0000000000..f63ac924c3 --- /dev/null +++ b/yt/visualization/tests/test_offaxisprojection_pytestonly.py @@ -0,0 +1,169 @@ +from typing import Union + +import numpy as np +import pytest +import unyt + +from yt.testing import ( + assert_rel_equal, + cubicspline_python, + fake_sph_flexible_grid_ds, + integrate_kernel, +) +from yt.visualization.api import ProjectionPlot + + +@pytest.mark.parametrize("weighted", [True, False]) +@pytest.mark.parametrize("periodic", [True, False]) +@pytest.mark.parametrize("depth", [None, (1.0, "cm"), (0.5, "cm")]) +@pytest.mark.parametrize("shiftcenter", [False, True]) +@pytest.mark.parametrize("northvector", [None, (1.0e-4, 1.0, 0.0)]) +def test_sph_proj_general_offaxis( + northvector: Union[tuple[float, float, float], None], + shiftcenter: bool, + depth: Union[tuple[float, str], None], + periodic: bool, + weighted: bool, +) -> None: + """ + Same as the on-axis projections, but we rotate the basis vectors + to test whether roations are handled ok. the rotation is chosen to + be small so that in/exclusion of particles within bboxes, etc. + works out the same way. + We just send lines of sight through pixel centers for convenience. + Particles at [0.5, 1.5, 2.5] (in each coordinate) + smoothing lengths 0.25 + all particles have mass 1., density 1.5, + except the single center particle, with mass 2., density 3. + + Parameters: + ----------- + northvector: tuple + y-axis direction in the final plot (direction vector) + shiftcenter: bool + shift the coordinates to center the projection on. + (The grid is offset to this same center) + depth: float or None + depth of the projection slice + periodic: bool + assume periodic boundary conditions, or not + weighted: bool + make a weighted projection (density-weighted density), or not + + Returns: + -------- + None + """ + if shiftcenter: + center = np.array((0.625, 0.625, 0.625)) # cm + else: + center = np.array((1.5, 1.5, 1.5)) # cm + bbox = unyt.unyt_array(np.array([[0.0, 3.0], [0.0, 3.0], [0.0, 3.0]]), "cm") + hsml_factor = 0.5 + unitrho = 1.5 + + # test correct centering, particle selection + def makemasses(i, j, k): + if i == j == k == 1: + return 2.0 + else: + return 1.0 + + # result shouldn't depend explicitly on the center if we re-center + # the data, unless we get cut-offs in the non-periodic case + # *almost* the z-axis + # try to make sure dl differences from periodic wrapping are small + epsilon = 1e-4 + projaxis = np.array([epsilon, 0.00, np.sqrt(1.0 - epsilon**2)]) + e1dir = projaxis / np.sqrt(np.sum(projaxis**2)) + # TODO: figure out other (default) axes for basis vectors here + if northvector is None: + e2dir = np.array([0.0, 1.0, 0.0]) + else: + e2dir = np.asarray(northvector) + e2dir = e2dir - np.sum(e1dir * e2dir) * e2dir # orthonormalize + e2dir /= np.sqrt(np.sum(e2dir**2)) + e3dir = np.cross(e1dir, e2dir) + + ds = fake_sph_flexible_grid_ds( + hsml_factor=hsml_factor, + nperside=3, + periodic=periodic, + offsets=np.full(3, 0.5), + massgenerator=makemasses, + unitrho=unitrho, + bbox=bbox.v, + recenter=center, + e1hat=e1dir, + e2hat=e2dir, + e3hat=e3dir, + ) + + source = ds.all_data() + # couple to dataset -> right unit registry + center = ds.arr(center, "cm") + # print('position:\n', source['gas','position']) + + # m / rho, factor 1. / hsml**2 is included in the kernel integral + # (density is adjusted, so same for center particle) + prefactor = 1.0 / unitrho # / (0.5 * 0.5)**2 + dl_cen = integrate_kernel(cubicspline_python, 0.0, 0.25) + + if weighted: + toweight_field = ("gas", "density") + else: + toweight_field = None + # we don't actually want a plot, it's just a straightforward, + # common way to get an frb / image array + prj = ProjectionPlot( + ds, + projaxis, + ("gas", "density"), + width=(2.5, "cm"), + weight_field=toweight_field, + buff_size=(5, 5), + center=center, + data_source=source, + north_vector=northvector, + depth=depth, + ) + img = prj.frb.data[("gas", "density")] + if weighted: + # periodic shifts will modify the (relative) dl values a bit + expected_out = np.zeros( + ( + 5, + 5, + ), + dtype=img.v.dtype, + ) + expected_out[::2, ::2] = unitrho + if depth is None: + expected_out[2, 2] *= 1.5 + else: + # only 2 * unitrho element included + expected_out[2, 2] *= 2.0 + else: + expected_out = np.zeros( + ( + 5, + 5, + ), + dtype=img.v.dtype, + ) + expected_out[::2, ::2] = dl_cen * prefactor * unitrho + if depth is None: + # 3 particles per l.o.s., including the denser one + expected_out *= 3.0 + expected_out[2, 2] *= 4.0 / 3.0 + else: + # 1 particle per l.o.s., including the denser one + expected_out[2, 2] *= 2.0 + # grid is shifted to the left -> 'missing' stuff at the left + if (not periodic) and shiftcenter: + expected_out[:1, :] = 0.0 + expected_out[:, :1] = 0.0 + # print(axis, shiftcenter, depth, periodic, weighted) + # print("expected:\n", expected_out) + # print("recovered:\n", img.v) + assert_rel_equal(expected_out, img.v, 4) diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 392749633e..3c2359a3f8 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -30,6 +30,7 @@ def off_axis_projection( no_ghost=False, interpolated=False, north_vector=None, + depth=None, num_threads=1, method="integrate", ): @@ -80,6 +81,10 @@ def off_axis_projection( north_vector : optional, array_like, default None A vector that, if specified, restricts the orientation such that the north vector dotted into the image plane points "up". Useful for rotations + depth: float, tuple[float, str], or unyt_array of size 1. + specify the depth of the projection region (size along the + line of sight). If no units are given (unyt_array or second + tuple element), code units are assumed. num_threads: integer, optional, default 1 Use this many OpenMP threads during projection. method : string @@ -145,6 +150,16 @@ def off_axis_projection( center = data_source.ds.arr(center, "code_length") if not hasattr(width, "units"): width = data_source.ds.arr(width, "code_length") + if depth is not None: + # handle units (intrinsic or as a tuple), + # then convert to code length + # float -> assumed to be in code units + if isinstance(depth, tuple): + depth = data_source.ds.arr(np.array([depth[0]]), depth[1]) + if hasattr(depth, "units"): + depth = depth.to("code_length").d + + # depth = data_source.ds.arr(depth, "code_length") if hasattr(data_source.ds, "_sph_ptypes"): if method != "integrate": @@ -203,16 +218,29 @@ def off_axis_projection( buf = np.zeros((resolution[0], resolution[1]), dtype="float64") mask = np.ones_like(buf, dtype="uint8") - x_min = center[0] - width[0] / 2 - x_max = center[0] + width[0] / 2 - y_min = center[1] - width[1] / 2 - y_max = center[1] + width[1] / 2 - z_min = center[2] - width[2] / 2 - z_max = center[2] + width[2] / 2 + ## width from fixed_resolution.py is just the size of the domain + # x_min = center[0] - width[0] / 2 + # x_max = center[0] + width[0] / 2 + # y_min = center[1] - width[1] / 2 + # y_max = center[1] + width[1] / 2 + # z_min = center[2] - width[2] / 2 + # z_max = center[2] + width[2] / 2 + + periodic = data_source.ds.periodicity + le = data_source.ds.domain_left_edge.to("code_length").d + re = data_source.ds.domain_right_edge.to("code_length").d + x_min, y_min, z_min = le + x_max, y_max, z_max = re + bounds = [x_min, x_max, y_min, y_max, z_min, z_max] + # only need (rotated) x/y widths + _width = (width.to("code_length").d)[:2] finfo = data_source.ds.field_info[item] ounits = finfo.output_units - bounds = [x_min, x_max, y_min, y_max, z_min, z_max] - + kernel_name = None + if hasattr(data_source.ds, "kernel_name"): + kernel_name = data_source.ds.kernel_name + if kernel_name is None: + kernel_name = "cubic" if weight is None: for chunk in data_source.chunks([], "io"): off_axis_projection_SPH( @@ -224,12 +252,15 @@ def off_axis_projection( chunk[ptype, "smoothing_length"].to("code_length").d, bounds, center.to("code_length").d, - width.to("code_length").d, + _width, + periodic, chunk[item].in_units(ounits), buf, mask, normal_vector, north, + depth=depth, + kernel_name=kernel_name, ) # Assure that the path length unit is in the default length units @@ -262,13 +293,16 @@ def off_axis_projection( chunk[ptype, "smoothing_length"].to("code_length").d, bounds, center.to("code_length").d, - width.to("code_length").d, + _width, + periodic, chunk[item].in_units(ounits), buf, mask, normal_vector, north, weight_field=chunk[weight].in_units(wounits), + depth=depth, + kernel_name=kernel_name, ) for chunk in data_source.chunks([], "io"): @@ -281,12 +315,15 @@ def off_axis_projection( chunk[ptype, "smoothing_length"].to("code_length").d, bounds, center.to("code_length").d, - width.to("code_length").d, + _width, + periodic, chunk[weight].to(wounits), weight_buff, mask, normal_vector, north, + depth=depth, + kernel_name=kernel_name, ) normalization_2d_utility(buf, weight_buff) @@ -300,6 +337,7 @@ def off_axis_projection( "north_vector": north_vector, "normal_vector": normal_vector, "width": width, + "depth": depth, "units": funits, "type": "SPH smoothed projection", } @@ -487,7 +525,8 @@ def temp_weightfield(field, data): image *= dl else: mask = image[:, :, 1] == 0 - image[:, :, 0] /= image[:, :, 1] + nmask = np.logical_not(mask) + image[:, :, 0][nmask] /= image[:, :, 1][nmask] image[mask] = 0 return image[:, :, 0] diff --git a/yt/visualization/volume_rendering/old_camera.py b/yt/visualization/volume_rendering/old_camera.py index cca5ff178a..fd171dbc80 100644 --- a/yt/visualization/volume_rendering/old_camera.py +++ b/yt/visualization/volume_rendering/old_camera.py @@ -2439,6 +2439,8 @@ def _render(self, double_check, num_threads, image, sampler, msg): data_object_registry["stereospherical_camera"] = StereoSphericalCamera +# replaced in volume_rendering API by the function of the same name in +# yt/visualization/volume_rendering/off_axis_projection def off_axis_projection( ds, center, diff --git a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py index dc5f4ae13a..57f6f85d52 100644 --- a/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py +++ b/yt/visualization/volume_rendering/tests/test_off_axis_SPH.py @@ -26,6 +26,7 @@ def test_no_rotation(): width = right_edge - left_edge px = ad["all", "particle_position_x"] py = ad["all", "particle_position_y"] + pz = ad["all", "particle_position_y"] hsml = ad["all", "smoothing_length"] quantity_to_smooth = ad["gas", "density"] density = ad["io", "density"] @@ -38,7 +39,7 @@ def test_no_rotation(): ds, center, normal_vector, width, resolution, ("gas", "density") ) pixelize_sph_kernel_projection( - buf2, mask, px, py, hsml, mass, density, quantity_to_smooth, bounds + buf2, mask, px, py, pz, hsml, mass, density, quantity_to_smooth, bounds ) assert_almost_equal(buf1.ndarray_view(), buf2)