Common Mesh Operations

In this section, we detail some common operations for manipulating data defined on a mesh via a MeshSource object.

Previewing the Mesh

The MeshSource.preview() function allows users to preview a low-resolution mesh by resampling the mesh and gathering the mesh to all ranks. It can also optionally project the mesh across multiple axes, which enables visualizations of the projected density field for quick data inspection by the user.

For example, below we initialize a LinearMesh object on a \(128^3\) mesh and preview the mesh on a \(64^3\) mesh after projecting the field along two axes:

In [1]: from nbodykit.lab import LinearMesh, cosmology

In [2]: from matplotlib import pyplot as plt

In [3]: cosmo = cosmology.Planck15

In [4]: Plin = cosmology.EHPower(cosmo, redshift=0)

In [5]: mesh = LinearMesh(Plin, Nmesh=128, BoxSize=1380, seed=42)

In [6]: density = mesh.preview(Nmesh=64, axes=(0,1))

In [7]: plt.imshow(density)
Out[7]: <matplotlib.image.AxesImage at 0x7f2fb964e5f8>
../_images/density-preview.png

Note

The previewed mesh result is broadcast to all ranks, so each rank allocates \(\mathrm{Nmesh}^3\) in memory.

Saving and Loading a Mesh

The MeshSource.save() function paints the mesh via a call to MeshSource.paint() and saves the mesh using a bigfile format. The output mesh saved to file can be in either configuration space or Fourier space, by specifying mode as either real or complex.

Below, we save our LinearMesh to a bigfile file:

# save the RealField
In [8]: mesh.save('linear-mesh-real.bigfile', mode='real', dataset='Field')

# save the ComplexField
In [9]: mesh.save('linear-mesh-complex.bigfile', mode='real', dataset='Field')

The saved mesh can be loaded from disk using the BigFileMesh class:

In [10]: from nbodykit.lab import BigFileMesh

In [11]: import numpy

# load the mesh in the form of a RealField
In [12]: real_mesh = BigFileMesh('linear-mesh-real.bigfile', 'Field')

# return the RealField via paint
In [13]: rfield = real_mesh.paint(mode='real')

# load the mesh in the form of a ComplexField
In [14]: complex_mesh = BigFileMesh('linear-mesh-complex.bigfile', 'Field')

# FFT to get the ComplexField as a RealField
In [15]: rfield2 = complex_mesh.paint(mode='real')

# the two RealFields must be the same!
In [16]: numpy.allclose(rfield.value, rfield2.value)
Out[16]: True

Here, we load our meshes in configuration space and Fourier space and then paint both with mode=real and verify that the results are the same.

Applying Functions to the Mesh

nbodykit supports performing transformations to the mesh data by applying arbitrary functions in either configuration space or Fourier space. Users can use the MeshSource.apply() function to apply these transformations. The function applied to the mesh should take two arguments, x and v:

  1. The x argument provides a list of length three holding the coordinate arrays that define the mesh. These arrays broadcast to the full shape of the mesh, i.e., they have shapes \((N_x,1,1)\), \((1,N_y,1)\), and \((1,1,N_z)\) if the mesh has shape \((N_x, N_y, N_z)\).
  2. The v argument is the array holding the value of the mesh field at the coordinate arrays in x

The units of the x coordinate arrays depend upon the values of the kind and mode keywords passed to the apply() function. The various cases are:

mode kind range of the “x” argument
real relative \([-L/2, L/2)\)
real index \([0, N)\)
complex wavenumber \([- \pi N/L, \pi N / L)\)
complex circular \([-\pi, \pi)\)
complex index \([0, N)\)

Here, \(L\) is the size of the box and N is the number of cells per mesh side.

One common use of the MeshSource.apply() functionality is applying compensation function to the mesh to correct for the interpolation window. The table of built-in compensation functions in the Compensation: Deconvolving the Window Kernel section provide examples of the syntax needed to apply functions to the mesh.

In the example below, we apply a filter function in Fourier space that divides the mesh by the squared norm of the wavenumber k on the mesh, and then print out the first few mesh cells of the filtered mesh to verify the function was applied properly.

In [17]: def filter(k, v):
   ....:     kk = sum(ki ** 2 for ki in k) # k^2 on the mesh
   ....:     kk[kk == 0] = 1
   ....:     return v / kk # divide the mesh by k^2
   ....: 

# apply the filter and get a new mesh
In [18]: filtered_mesh = mesh.apply(filter, mode='complex', kind='wavenumber')

# get the filtered RealField object
In [19]: filtered_rfield = filtered_mesh.paint(mode='real')

In [20]: print("head of filtered Realfield = ",  filtered_rfield[:10,0,0])
head of filtered Realfield =  [  849.83032227   836.96740723   859.5670166    887.95947266   920.81628418
   985.9119873   1063.7590332   1141.55944824  1226.26428223  1310.90649414]

In [21]: print("head of original RealField = ",  rfield[:10,0,0])
head of original RealField =  [ 0.10047328  0.83333337  1.72612357  2.29414558  1.68662632  1.89175129
  1.43605208  0.86345029  0.87090874  1.63552475]

Resampling a Mesh

Users can resample a mesh by specifying the Nmesh keyword to the MeshSource.paint() function. For example, below we resample a LinearMesh object, changing the mesh resolution from Nmesh=128 to Nmesh=32.

In [22]: from nbodykit.lab import LinearMesh, cosmology

# linear mesh
In [23]: Plin = cosmology.EHPower(cosmology.Planck15, redshift=0.55)

In [24]: source = LinearMesh(Plin, Nmesh=64, BoxSize=512, seed=42)

# paint, re-sampling to Nmesh=32
In [25]: real = source.paint(mode='real', Nmesh=32)

In [26]: print("original Nmesh = ", source.attrs['Nmesh'])
original Nmesh =  [64 64 64]

In [27]: print("resampled Nmesh = ", real.Nmesh)
resampled Nmesh =  [32 32 32]

In [28]: print("shape of resampled density field = ", real.cshape)
shape of resampled density field =  [32 32 32]