How to convert particles' velocities to the gridded velocity field?
WangYun1995 opened this issue · 4 comments
Hello,
I have 1820^3 particles and their spatial velocities ( unit: km/s, shape: (1820^3, 3) ). Then I want to convert these 'scattered' velocity data to the 1024^3 mesh data, such that I can do some further computations on it.
What I know is converting particle masses to the gridded density field by calling to_mesh(), and setting weight='Masses'.
By reading the docs of nbodykit, I also know that I can obtain the mass-weighted velocities, i.e. momentum field per unit mass by setting value = 'velocity_x'. However, what I want is V(x), not the weighted value, i.e. (1+delta)*V(x). So whether does nbodykit offer the corresponding functionality?
We usually prefer to avoid using the velocity mesh because the velocity is undefined when there is no mass in the cell.
If you really do want the velocity mesh, then I guess you already have doodled the solution:
1 wv = (1 + delta) * V(x),
2 w = (1 + delta).
3 v = wv / w.
It may feel a bit quirky at first, but this way the potential problem is more obvious -- "when divide, check for zeros".
For example, you can can say that when w == 0, v == 0, which is a good choice if eventually v is multipled by w again.
Alternatively one can make the resampling window adaptive (e.g. Voronoi tessellation), such that we never run to a zero in w.
We don't have that in nbodykit -- I think it is possible to write one using scipy's qhull interface.
Thank you for your clarification. I still have another question.
I have obtained one momentum field (x-component) from scattered data of IllustrisTNGby calling to_mesh(....., weight='Masses', value='V_x'). However, the maximum value of momentum, i.e. [(1+delta)*V_x] is proportional to 10^7 km/s. After the operation '3 v = wv / w' is done, I found that the maximum of v is proportional to 10^8 km/s greater than light speed. I can't understand why this situation occured.
It could be a unit conversion error too. How does this compare to the velocity value you have from the snapshot file ('V_x'). Are the close?
Could you also check the histogram and standard-dev?
also worth checking if there is an off by 1 error. I don't immediately recall if the field you get after paint is 1+delta or delta. (checking the mean would tell).
It could be a unit conversion error too. How does this compare to the velocity value you have from the snapshot file ('V_x'). Are the close?
Could you also check the histogram and standard-dev?
also worth checking if there is an off by 1 error. I don't immediately recall if the field you get after paint is 1+delta or delta. (checking the mean would tell).
According to your advice, I computed the mean and std of the
- z-component of velocities (from the snapshot file),
- z-component of momentum field (mass-weighted V_z),
- overdensity
- mass_weighted V_z / overdensity.
The results are shown in below.
Obviously, The results of 'Mass weighted V_z / Overdensity' are ridiculous. The code snippet I used to compute 'Mass weighted V_z / Overdensity' is
V_z = np.divide(mass_weighted_vz, overdensity, out=np.zeros(mass_weighted_vz.shape, dtype='float32'), where=overdensity!=0.0)
The histogram of 'z-component of velocities (from the snapshot file)':
The histogram of 'z-component of momentum field (mass-weighted V_z)':
