[yt-4.0] sph deposition normalization

[AshKelly]

Bug report

Bug summary

Robin Kooistra reported in the slack that the density from the mean SPH gather/scatter of gas was way of what it should be e.g. critical density.

The reason for this is that we over-normalize the results of the deposition after we have calculated the density.

Namely, the yt depositon for gas density follows Eq. 3 of the SPLASH method paper, whereas most codes follow Eq. 9. The choice of Eq. 3 leads to two issues:

1) The normalization is done automatically. However, the accuracy of the normalization depends m/h^3 == rho which is not always true for SPH simulations, and in fact, rarely is!

2) As Eq. 3 already returns the normalized field, when we normalize we are normalizing an already normalized field.

SPLASH

I propose we change to the typical choice of Eq. 9 which will involve a few minor changes in yt.utilities.lib.pixelization_routines.pyx. I plan to work on this in the #pullrequest-triage on May 8th.

Code for reproduction

import numpy as np
import h5py
from astropy.cosmology import FlatLambdaCDM
from astropy import units as u
import yt
sim = 'TNG100-3' #path to the simulation files
snap = 33 #snapshot number
#extracting cosmology, etc.
header = h5py.File(sim+'/snap_0%d.0.hdf5'%snap,'r')['Header'].attrs
rs = header['Redshift']
maxchunks = header['NumFilesPerSnapshot']
L = header['BoxSize']
cosmo = FlatLambdaCDM(H0=header['HubbleParam']*100.,Om0=header['Omega0'],Ob0=header['OmegaBaryon'],Tcmb0=2.725)
rhocrit_bar = cosmo.critical_density(rs)*cosmo.Ob(rs)  #critical_density in cgs units
rhocrit_bar_simunits = rhocrit_bar.to(1.e10*u.M_sun/cosmo.h/((u.kpc/cosmo.h)**3)) #critical density in simulation units
#loading necessary parameters from hdf5 files into memory
posx = []
posy = []
posz = []
dens = []
m = []
smolen = []
for chunk in range(maxchunks):
  data = h5py.File(sim+'/snap_0%d.%d.hdf5'%(snap,chunk),'r')  
  coords = np.array(data['PartType0/Coordinates'])
  density = np.array(data['PartType0/Density'])
  mass = np.array(data['PartType0/Masses'])
  posx.append(coords[:,0])
  posy.append(coords[:,1])
  posz.append(coords[:,2])
  dens.append(density)
  m.append(mass)
  if (sim == 'Illustris-1') or (sim == 'Illustris-3'):
    smolen.append(np.array(data['PartType0/SmoothingLength']))
  elif (sim == 'TNG100-1') or (sim == 'TNG100-3'):
    smolen.append((3./4./np.pi*mass/density)**(1./3.))
posx = np.hstack(posx).astype('float64')
posy = np.hstack(posy).astype('float64')
posz = np.hstack(posz).astype('float64')
dens = np.hstack(dens).astype('float64')
m = np.hstack(m).astype('float64')
smolen = np.hstack(smolen).astype('float64')
bbox = np.array([[0.,L],[0.,L],[0.,L]])
data = dict(density = dens,
            particle_position_x = posx,
            particle_position_y = posy,
            particle_position_z = posz,
            particle_mass = m,
            smoothing_length = smolen
            )
ds = yt.load_particles(data,bbox=bbox)
ds.num_neighbors = 16
ds.sph_smoothing_style = 'gather'
N = 600 #number of grid-cells per side
ag = ds.arbitrary_grid(ds.domain_left_edge,ds.domain_right_edge,dims=[N,N,N])
mesh_dens = ag['io','density']

Version Information

• yt version: yt-4.0

[AshKelly]

@matthewturk @JBorrow feel free to pitch in ideas or suggestions for this!

[munkm]

I plan to work on this in the #pullrequest-triage on May 8th.

Excellent!!!!

[jzuhone]

@AshKelly did we ever resolve this? I don't think that we did.

[JBorrow]

No, we didn't. Though me and Ash did write a paper about this at some point... https://ui.adsabs.harvard.edu/abs/2021arXiv210605281B/abstract

[jzuhone]

@JBorrow thanks for this--I think on some reasonable time scale we should try implementing this in yt. I skimmed the paper but do you think it would be very difficult or can you think of any potential hangups? I'm happy to take a whack at it myself but I'll probably bug you with questions.

[JBorrow]

It's very simple, all you're doing is taking the first order integration of the kernel to the grid (e.g. W(r_center_of_pixel, h) being the contribution to that pixel) and turning it into something either higher order, or like I did here just creating more resolution elements to do this first order integration over. It boils down to a case of:

if (pixel_overlap_of_kernel < threshold): create_fake_refined_grid() project_onto_fake_refined_grid() derefine_grid_and_add_to(main_grid)

The code for swiftismio is super simple too, as it's all numba accelerated:

Base case (no refinement): https://github.com/SWIFTSIM/swiftsimio/blob/master/swiftsimio/visualisation/projection_backends/fast.py Subsampled case (with refinement): https://github.com/SWIFTSIM/swiftsimio/blob/master/swiftsimio/visualisation/projection_backends/subsampled.py