"It's not a mesh; it's a meshoid!" - Alfred B. Einstimes

Meshoid is a Python package for analyzing meshless particle data, performing such operations as density/smoothing length computation, numerical differentiation on unstructured data, and various kernel-weighted projection and averaging operations. It was originally designed to work with simulation data from the GIZMO code, but is totally agnostic to the actual format of the data it works on, so feel free to use it with your favorite meshless or unstructured mesh code!

Full Documentation

Installation

Install from PyPI via pip install meshoid or pull this repo and run pip install . from the repo directory.

Walkthrough

First let's import pylab, Meshoid, and the loadfromsnapshot script for loading GIZMO outputs.

python from matplotlib import pyplot as plt import numpy as np from meshoid import Meshoid from load_from_snapshot import load_from_snapshot %matplotlib inline

Now let's load some of the gas data fields from a snapshot in the public FIRE data release using loadfromsnapshot. In this case we'll perform a density cut at n_H ~ .1 cm^-3 to narrow it down to just the ISM.

```python from os import system from os.path import isfile

download the data - feel free to use your own snapshots!

for i in range(4): if not isfile(f"./snapshot600.{i}.hdf5"): system(f"wget -r -nH --cut-dirs=6 --no-parent --reject='index.html*' -e robots=off https://users.flatironinstitute.org/~chayward/fire2publicrelease/core/m12ires7100/output/snapdir600/snapshot600.{i}.hdf5")

rho = loadfromsnapshot("Density", 0, ".", 600) RHOTONH = 300 # convert to H number density (approx) densitycut = (rho*RHOTONH > .01) pdata = {} fieldstoload = "Masses", "Coordinates", "SmoothingLength", "Velocities", "Density" pdata = {field: loadfromsnapshot(field, 0, ".", 600, particlemask=np.arange(len(rho))[densitycut]) for field in fieldsto_load} ```

Finally, before getting to the meshoid stuff we will also center the coordinates, perform a cut in galactocentric radius at 40kpc, and orient our coordinates to the principal axes of the gas distribution.

```python

do centering

pos = pdata["Coordinates"] center = np.median(pos,axis=0) pos -= center

do radius cut

MAXRADIUS = 40. radiuscut = np.sum(pos*pos,axis=1) < MAXRADIUS * MAXRADIUS pos, mass, hsml, v, rho = pos[radiuscut], pdata["Masses"][radius_cut], pdata["SmoothingLength"][radius_cut], pdata["Velocities"][radius_cut], pdata["Density"][radius_cut] centernew = np.average(pos,axis=0,weights=mass) # now do another re-centering center += centernew pos -= centernew

now get the principal axes - eigenvectors of the second mass moment matrix

covpos = np.cov(pos.T, aweights=mass) w, coordinatebasis = np.linalg.eigh(covpos) coordinatebasis = coordinatebasis[:,w.argsort()[::-1]] # sort so the smallest moment axis is the last = z-axis pos = pos @ coordinatebasis # dot product with each basis vector to get coordinates in new basis ```

OK, now let's start by making a map of gas surface density. We can do so by generating a Meshoid object from the particle masses, coordinates, and smoothing lengths, and then calling the SurfaceDensity method. This function is useful for giving kernel-weighted projected quantities on a Cartesion grid of sighlines.

Meshoid can also adaptively compute particle smoothing lengths on-the-fly provided only the coordinates, but it requires a nearest-neighbor search that takes a while so it's best to provide the smoothing length when you can.

```python import matplotlib.colors as colors from mpltoolkits.axesgrid1 import makeaxeslocatable

M = Meshoid(pos, mass, hsml) rmax = 20 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') fig, ax = plt.subplots(figsize=(6,6)) sigmagasmsunpc2 = M.SurfaceDensity(M.m,center=np.array([0,0,0]),size=40.,res=res)*1e4 p = ax.pcolormesh(X, Y, sigmagasmsunpc2, norm=colors.LogNorm(vmin=.1,vmax=1e3)) divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$\Sigma{\rm gas}$ $(\rm M\odot\,pc^{-2})$",pad=0,cax=cax) ax.setaspect('equal') ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.tightlayout(hpad=1) plt.show() ```

Now let's look at the 3D gas density in a slice through the galaxy, using the Slice method. This will reconstruct the data to a grid of points in a plane slicing through the data. You can chose the order of the reconstruction: 0 will simply give the value of the nearest particle (i.e. reflecting the Voronoi domains), 1 will perform a linear reconstruction from that particle, etc. The best order will depend upon the nature of the data: smooth data will look best with higher-order reconstructions, while messier data will have nasty overshoots and artifacts. Here the density field is quite poorly-resolved, so we will use a 0'th order reconstruction.

```python fig, ax = plt.subplots(figsize=(6,6)) densityslicenHcgs = M.Slice(rho,center=np.array([0,0,0]),size=40.,res=res, order=0) * RHOTONH

alternative to try: default linear reconstruction of log rho to avoid overshoot to negative density

densityslicenHcgs = 10**M.Slice(np.log10(rho),center=np.array([0,0,0]),size=40.,res=res) * RHOTONH

p = ax.pcolormesh(X, Y, densityslicenHcgs,norm=colors.LogNorm(vmin=.01,vmax=1e2)) ax.setaspect('equal') divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$n{\rm H}$ $(\rm cm^{-3})$",cax=cax) ax.setxlabel("X (kpc)") ax.set_ylabel("Y (kpc)") plt.show() ```

Simple Radiative Transfer

Meshoid is also capable of performing radiative transfer with a known emissivity/source function and opacity, neglecting scattering. For example, we can load in the stellar positions and assume a simple constant mass-to-light ratio, and calculate the dust-extincted starlight in the V-band.

```python from meshoid.radiation import radtransfer, dustabsopacity from astropy import units as u, constants as c kappadustcodeunits = dustabsopacity(0.555) * (u.cm2/u.g).to(u.kpc2/(1e10*c.Msun)) # dust opacity in cgs converted to solar - evaluated at 555nm kappagas = np.repeat(kappadustcodeunits,len(mass)) jgas = np.zeroslike(mass) # assume dust does not emit

have to get the star properties now

xstars = (loadfromsnapshot("Coordinates", 4, ".", 600) - center) @ coordinatebasis mstars = loadfromsnapshot("Masses", 4, ".", 600) hstar = np.repeat(0.1, len(mstars)) # 100pc radii MASSTOLIGHTSOLAR = 1. # emissivity is just the light-to-mass ratio for stars - here assume 1 (old-ish stellar population in V-band) jstar = np.repeat(1e10/(MASSTOLIGHTSOLAR), len(mstars)) # we are assuming a constant emissivity throughout the kernel-smoothed star particles kappastars = np.zeros(len(m_stars))

now combine all emissivities, opacities, masses, kernel lengths

jall = np.atleast2d(np.concatenate([jgas, jstar])).T # 2D because this has shape (numparticles, numbands) (can have an arbitrary number of bands) kappaall = np.atleast2d(np.concatenate([kappagas, kappastars])).T # ditto kappaall = kappaall.clip(1e-100) # we divide by kappa at a certain point so put this to avoid floating-point errors hall = np.concatenate([hsml, hstar]) mall = np.concatenate([mass, mstars]) xall = np.concatenate([pos,xstars],axis=0)

rmax = 10 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') I = radtransfer(jall, mall,kappaall, xall,h_all,res,2*rmax) # actual call to rad transfer solver

screw you I'm not converting this to mag/arcsec^2

fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X,Y, I[:,:,0], norm=colors.LogNorm(vmin=I.max()/1e3, vmax=I.max()), cmap="Greysr") ax.setaspect('equal') divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$I\,\left(L\odot\,\rm kpc^{-2}\,\rm sr^{-1}\right)\, $",cax=cax) ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.show()

```

Differential Operators

Now let's play around with Meshoid's numerical differentiation. Meshoid can take both first (Meshoid.D) and second derivatives (Meshoid.D2) on unstructured data, using a kernel-weighted (or unweighted) least-squares gradient estimator.

As a first sanity check, we can try differentiating the coordinate functions, with respect to those coordinates. That ought to return an identity matrix. Note that you can differentiate scalars, vectors, or even arbitrary tensors that are defined on the meshoid. In general, differentiating a tensor of rank N will return a tensor of rank N+1.

The first time a given differentiation method is called, Meshoid can take a minute to compute the weights that it needs. Hang in there, Meshoid is working diligently and will re-use those weights the next time you need a derivative!

python M.D(pos)

array([[[ 1.00000000e+00, -1.56125113e-17,  1.56125113e-17],
        [ 2.25080371e-16,  1.00000000e+00, -2.12069945e-16],
        [ 1.30104261e-16, -3.03576608e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -5.96514483e-17, -1.77402580e-17],
        [-5.62429877e-17,  1.00000000e+00,  8.52725009e-17],
        [ 2.69437792e-16, -8.81591892e-17,  1.00000000e+00]],

       [[ 1.00000000e+00,  4.98157017e-17,  4.32874494e-16],
        [ 4.84079329e-17,  1.00000000e+00, -7.06018902e-17],
        [ 4.68917440e-18,  4.95344868e-17,  1.00000000e+00]],

       ...,

       [[ 1.00000000e+00, -1.60970141e-17, -2.35271871e-17],
        [-3.77099068e-17,  1.00000000e+00,  3.52535113e-17],
        [ 6.46048970e-17, -2.99888627e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -2.06011965e-16,  1.01874347e-16],
        [-2.15403867e-16,  1.00000000e+00,  9.90893023e-17],
        [ 1.01969214e-16,  5.41396355e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -7.45388994e-19, -2.98697699e-17],
        [-2.60208521e-16,  1.00000000e+00, -1.51354623e-16],
        [ 6.80011603e-16, -1.39645240e-16,  1.00000000e+00]]])

OK now let's look at something physical. Let's calculate the enstrophy, which is just the norm squared of the velocity gradient, and plot it as a slice.

python gradv = M.D(v) enstrophy = np.sum(gradv*gradv, axis=(1,2)) enstrophy_projection = M.Slice(enstrophy,center=np.array([0,0,0]),size=40.,res=res,order=0) fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X, Y, enstrophy_projection*.979**2, norm=colors.LogNorm(vmin=10,vmax=1e7)) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"Enstrophy $(\rm Gyr^{-2})$",cax=cax) ax.set_aspect('equal') ax.set_xlabel("X (kpc)") ax.set_ylabel("Y (kpc)") plt.show()

Walkthrough

First let's import pylab, Meshoid, and the loadfromsnapshot function for loading GIZMO outputs.

python from matplotlib import pyplot as plt import numpy as np from meshoid import Meshoid, load_from_snapshot %matplotlib inline

Now let's load some of the gas data fields from a snapshot in the public FIRE data release using loadfromsnapshot. In this case we'll perform a density cut at n_H ~ .1 cm^-3 to narrow it down to just the ISM.

```python from os import system from os.path import isfile

download the data - feel free to use your own snapshots!

snapnum = 600 for i in range(4): if not isfile(f"./snapshot{snapnum}.{i}.hdf5"): system(f"wget -r -nH --cut-dirs=6 --no-parent --reject='index.html*' -e robots=off https://users.flatironinstitute.org/~mgrudic/fire2publicrelease/core/m12ires7100/output/snapdir{snapnum}/snapshot{snapnum}.{i}.hdf5")

rho = loadfromsnapshot("Density", 0, ".", snapnum) RHOTONH = 300 # convert to H number density (approx) densitycut = (rho*RHOTONH > .01) pdata = {} fieldstoload = "Masses", "Coordinates", "SmoothingLength", "Velocities", "Density" pdata = {field: loadfromsnapshot(field, 0, ".", snapnum, particlemask=np.arange(len(rho))[densitycut]) for field in fieldsto_load} ```

Finally, before getting to the meshoid stuff we will also center the coordinates, perform a cut in galactocentric radius at 40kpc, and orient our coordinates to the principal axes of the gas distribution.

```python

do centering

pos = pdata["Coordinates"] center = np.median(pos,axis=0) pos -= center

do radius cut

MAXRADIUS = 40. radiuscut = np.sum(pos*pos,axis=1) < MAXRADIUS * MAXRADIUS pos, mass, hsml, v, rho = pos[radiuscut], pdata["Masses"][radius_cut], pdata["SmoothingLength"][radius_cut], pdata["Velocities"][radius_cut], pdata["Density"][radius_cut] centernew = np.average(pos,axis=0,weights=mass) # now do another re-centering center += centernew pos -= centernew

now get the principal axes - eigenvectors of the second mass moment matrix

covpos = np.cov(pos.T, aweights=mass) w, coordinatebasis = np.linalg.eigh(covpos) coordinatebasis = coordinatebasis[:,w.argsort()[::-1]] # sort so the smallest moment axis is the last = z-axis pos = pos @ coordinatebasis # dot product with each basis vector to get coordinates in new basis ```

Projections

OK, now let's start by making a map of gas surface density. We can do so by generating a Meshoid object from the particle masses, coordinates, and smoothing lengths, and then calling the SurfaceDensity method. This function is useful for giving kernel-weighted projected quantities on a Cartesion grid of sighlines.

Meshoid can also adaptively compute particle smoothing lengths on-the-fly provided only the coordinates, but it requires a nearest-neighbor search that takes a while so it's best to provide the smoothing length when you can.

```python import matplotlib.colors as colors from mpltoolkits.axesgrid1 import makeaxeslocatable

M = Meshoid(pos, mass, hsml) rmax = 20 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') fig, ax = plt.subplots(figsize=(6,6)) sigmagasmsunpc2 = M.SurfaceDensity(M.m,center=np.array([0,0,0]),size=40.,res=res)*1e4 p = ax.pcolormesh(X, Y, sigmagasmsunpc2, norm=colors.LogNorm(vmin=.1,vmax=1e3)) divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$\Sigma{\rm gas}$ $(\rm M\odot\,pc^{-2})$",pad=0,cax=cax) ax.setaspect('equal') ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.tightlayout(hpad=1) plt.show() ```

Slices

Now let's look at the 3D gas density in a slice through the galaxy, using the Slice method. This will reconstruct the data to a grid of points in a plane slicing through the data. You can chose the order of the reconstruction: 0 will simply give the value of the nearest particle (i.e. reflecting the Voronoi domains), 1 will perform a linear reconstruction from that particle, etc. The best order will depend upon the nature of the data: smooth data will look best with higher-order reconstructions, while messier data will have nasty overshoots and artifacts. Here the density field is quite poorly-resolved, so a nearest-neighbor reconstruction is probably most appropriate.

```python fig, ax = plt.subplots(1,2,figsize=(6,6), sharex=True, sharey=True, gridspeckw={'widthratios': [1,1]}) densityslicenHcgs = M.Slice(rho,center=np.array([0,0,0]),size=40.,res=res, order=0) * RHOTONH

alternative to try: default linear reconstruction of log rho to avoid overshoot to negative density

M = Meshoid(pos, mass, hsml)

densityslicenHcgs2 = 10**M.Slice(np.log10(rho),center=np.array([0,0,0]),size=40.,res=res) * RHOTONH p = ax[0].pcolormesh(X, Y, densityslicenHcgs,norm=colors.LogNorm(vmin=.01,vmax=1e2)) p2 = ax[1].pcolormesh(X, Y, densityslicenHcgs2,norm=colors.LogNorm(vmin=.01,vmax=1e2))

setaxisargs = {"xlabel": "X (kpc)", "ylabel": "Y (kpc)","aspect": "equal"} for a in ax: a.set(**setaxisargs) divider = makeaxeslocatable(a) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$n{\rm H}$ $(\rm cm^{-3})$",cax=cax) ax[0].settitle(r"Nearest-neighbor $\rho$") ax[1].settitle(r"1st-order reconstruction of $\log \rho$") fig.tight_layout() plt.show() ```

Simple Radiative Transfer

Meshoid is also capable of performing radiative transfer with a known emissivity/source function and opacity, neglecting scattering. This is possible because there is a direct analogy between doing a line-integral through a density field (as in the project above), and solving the time-independent radiative transfer equation.

For example, we can load in the stellar positions and assume a simple constant mass-to-light ratio, and calculate the dust-extincted starlight in the V-band.

```python from meshoid.radiation import radtransfer, dustabsopacity from astropy import units as u, constants as c

kappadustcodeunits = dustabsopacity(0.555).to(u.kpc*2/(1e10c.Msun)).value # dust opacity in cgs converted to solar - evaluated at 555nm kappagas = np.repeat(kappadustcodeunits,len(mass)) jgas = np.zeroslike(mass) # assume dust does not emit

have to get the star properties now

xstars = (loadfromsnapshot("Coordinates", 4, ".", 600) - center) @ coordinatebasis mstars = loadfromsnapshot("Masses", 4, ".", 600) hstar = np.repeat(0.1, len(mstars)) # 100pc radii MASSTOLIGHTSOLAR = 1. # emissivity is just the light-to-mass ratio for stars - here assume 1 (old-ish stellar population in V-band) jstar = np.repeat(1e10/(MASSTOLIGHTSOLAR), len(mstars)) # we are assuming a constant emissivity throughout the kernel-smoothed star particles kappastars = np.zeros(len(m_stars))

now combine all emissivities, opacities, masses, kernel lengths

jall = np.atleast2d(np.concatenate([jgas, jstar])).T # 2D because this has shape (numparticles, numbands) (can have an arbitrary number of bands) kappaall = np.atleast2d(np.concatenate([kappagas, kappastars])).T # ditto kappaall = kappaall.clip(1e-100) # we divide by kappa at a certain point so put this to avoid floating-point errors hall = np.concatenate([hsml, hstar]) mall = np.concatenate([mass, mstars]) xall = np.concatenate([pos,xstars],axis=0)

rmax = 10 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') I = radtransfer(jall, mall,kappaall, xall,h_all,res,2*rmax) # actual call to rad transfer solver

screw you I'm not converting this to mag/arcsec^2

fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X,Y, I[:,:,0], norm=colors.LogNorm(vmin=I.max()/1e3, vmax=I.max()), cmap="Greysr") ax.setaspect('equal') divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$I\,\left(L\odot\,\rm kpc^{-2}\,\rm sr^{-1}\right)\, $",cax=cax) ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.show()

```

Differential Operators

Now let's play around with Meshoid's numerical differentiation. Meshoid can take both first (Meshoid.D) and second derivatives (Meshoid.D2) on unstructured data, using a kernel-weighted (or unweighted) least-squares gradient estimator.

As a first sanity check, we can try differentiating the coordinate functions, with respect to those coordinates. That ought to return an identity matrix. Note that you can differentiate scalars, vectors, or even arbitrary tensors that are defined on the meshoid. In general, differentiating a tensor of rank N will return a tensor of rank N+1.

The first time a given differentiation method is called, Meshoid can take a minute to compute the weights that it needs. Hang in there, Meshoid is working diligently and will re-use those weights the next time you need a derivative!

python M.D(pos)

array([[[ 1.00000000e+00, -1.56125113e-17,  1.56125113e-17],
        [ 2.25080371e-16,  1.00000000e+00, -2.12069945e-16],
        [ 1.30104261e-16, -3.03576608e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -5.96514483e-17, -1.77402580e-17],
        [-5.62429877e-17,  1.00000000e+00,  8.52725009e-17],
        [ 2.69437792e-16, -8.81591892e-17,  1.00000000e+00]],

       [[ 1.00000000e+00,  4.98157017e-17,  4.32874494e-16],
        [ 4.84079329e-17,  1.00000000e+00, -7.06018902e-17],
        [ 4.68917440e-18,  4.95344868e-17,  1.00000000e+00]],

       ...,

       [[ 1.00000000e+00, -1.60970141e-17, -2.35271871e-17],
        [-3.77099068e-17,  1.00000000e+00,  3.52535113e-17],
        [ 6.46048970e-17, -2.99888627e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -2.06011965e-16,  1.01874347e-16],
        [-2.15403867e-16,  1.00000000e+00,  9.90893023e-17],
        [ 1.01969214e-16,  5.41396355e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -7.45388994e-19, -2.98697699e-17],
        [-2.60208521e-16,  1.00000000e+00, -1.51354623e-16],
        [ 6.80011603e-16, -1.39645240e-16,  1.00000000e+00]]],
      shape=(2378703, 3, 3))

OK now let's look at something physical. Let's calculate the enstrophy, which is just the norm squared of the velocity gradient, and plot it as a slice.

python gradv = M.D(v) enstrophy = np.sum(gradv*gradv, axis=(1,2)) enstrophy_projection = M.Slice(enstrophy,center=np.array([0,0,0]),size=40.,res=res,order=0) fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X, Y, enstrophy_projection*.979**2, norm=colors.LogNorm(vmin=10,vmax=1e7)) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"Enstrophy $(\rm Gyr^{-2})$",cax=cax) ax.set_aspect('equal') ax.set_xlabel("X (kpc)") ax.set_ylabel("Y (kpc)") plt.show()

python M.particle_mask

array([      0,       1,       2, ..., 2378700, 2378701, 2378702],
      shape=(2378703,))

```python

```

Walkthrough

First let's import pylab, Meshoid, and the loadfromsnapshot function for loading GIZMO outputs.

python from matplotlib import pyplot as plt import numpy as np from meshoid import Meshoid, load_from_snapshot %matplotlib inline

Now let's load some of the gas data fields from a snapshot in the public FIRE data release using loadfromsnapshot. In this case we'll perform a density cut at n_H ~ .1 cm^-3 to narrow it down to just the ISM.

```python from os import system from os.path import isfile

download the data - feel free to use your own snapshots!

snapnum = 600 for i in range(4): if not isfile(f"./snapshot{snapnum}.{i}.hdf5"): system(f"wget -r -nH --cut-dirs=6 --no-parent --reject='index.html*' -e robots=off https://users.flatironinstitute.org/~mgrudic/fire2publicrelease/core/m12ires7100/output/snapdir{snapnum}/snapshot{snapnum}.{i}.hdf5")

rho = loadfromsnapshot("Density", 0, ".", snapnum) RHOTONH = 300 # convert to H number density (approx) densitycut = (rho*RHOTONH > .01) pdata = {} fieldstoload = "Masses", "Coordinates", "SmoothingLength", "Velocities", "Density" pdata = {field: loadfromsnapshot(field, 0, ".", snapnum, particlemask=np.arange(len(rho))[densitycut]) for field in fieldsto_load} ```

---------------------------------------------------------------------------

KeyboardInterrupt                         Traceback (most recent call last)

Cell In[2], line 15
     13 pdata = {}
     14 fields_to_load = "Masses", "Coordinates", "SmoothingLength", "Velocities", "Density"
---> 15 pdata = {field: load_from_snapshot(field, 0, ".", snapnum, particle_mask=np.arange(len(rho))[density_cut]) for field in fields_to_load}


Cell In[2], line 15, in <dictcomp>(.0)
     13 pdata = {}
     14 fields_to_load = "Masses", "Coordinates", "SmoothingLength", "Velocities", "Density"
---> 15 pdata = {field: load_from_snapshot(field, 0, ".", snapnum, particle_mask=np.arange(len(rho))[density_cut]) for field in fields_to_load}


File ~/code/meshoid/examples/load_from_snapshot.py:187, in load_from_snapshot(value, ptype, sdir, snum, particle_mask, axis_mask, units_to_physical, four_char, snapshot_name, snapdir_name, extension)
    185     q_t = numpy.array(file['PartType'+str(ptype)+'/'+value+'/']).take(axis_mask,axis=1)
    186 else:
--> 187     q_t = numpy.array(file['PartType'+str(ptype)+'/'+value+'/'])                    
    188 # check data has non-zero size
    189 if(q_t.size > 0): 
    190     # if this is the first time we are actually reading it, parse it and determine the shape of the vector, to build the data container


File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()


File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()


File ~/python_work/lib/python3.11/site-packages/h5py/_hl/dataset.py:1112, in Dataset.__array__(self, dtype, copy)
   1109 if self.size == 0:
   1110     return arr
-> 1112 self.read_direct(arr)
   1113 return arr


File ~/python_work/lib/python3.11/site-packages/h5py/_hl/dataset.py:1068, in Dataset.read_direct(self, dest, source_sel, dest_sel)
   1065     dest_sel = sel.select(dest.shape, dest_sel)
   1067 for mspace in dest_sel.broadcast(source_sel.array_shape):
-> 1068     self.id.read(mspace, fspace, dest, dxpl=self._dxpl)


KeyboardInterrupt: 

Finally, before getting to the meshoid stuff we will also center the coordinates, perform a cut in galactocentric radius at 40kpc, and orient our coordinates to the principal axes of the gas distribution.

```python

do centering

pos = pdata["Coordinates"] center = np.median(pos,axis=0) pos -= center

do radius cut

MAXRADIUS = 40. radiuscut = np.sum(pos*pos,axis=1) < MAXRADIUS * MAXRADIUS pos, mass, hsml, v, rho = pos[radiuscut], pdata["Masses"][radius_cut], pdata["SmoothingLength"][radius_cut], pdata["Velocities"][radius_cut], pdata["Density"][radius_cut] centernew = np.average(pos,axis=0,weights=mass) # now do another re-centering center += centernew pos -= centernew

now get the principal axes - eigenvectors of the second mass moment matrix

covpos = np.cov(pos.T, aweights=mass) w, coordinatebasis = np.linalg.eigh(covpos) coordinatebasis = coordinatebasis[:,w.argsort()[::-1]] # sort so the smallest moment axis is the last = z-axis pos = pos @ coordinatebasis # dot product with each basis vector to get coordinates in new basis ```

Projections

OK, now let's start by making a map of gas surface density. We can do so by generating a Meshoid object from the particle masses, coordinates, and smoothing lengths, and then calling the SurfaceDensity method. This function is useful for giving kernel-weighted projected quantities on a Cartesion grid of sighlines.

Meshoid can also adaptively compute particle smoothing lengths on-the-fly provided only the coordinates, but it requires a nearest-neighbor search that takes a while so it's best to provide the smoothing length when you can.

```python import matplotlib.colors as colors from mpltoolkits.axesgrid1 import makeaxeslocatable

M = Meshoid(pos, mass, hsml) rmax = 20 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') fig, ax = plt.subplots(figsize=(6,6)) sigmagasmsunpc2 = M.SurfaceDensity(M.m,center=np.array([0,0,0]),size=40.,res=res)*1e4 p = ax.pcolormesh(X, Y, sigmagasmsunpc2, norm=colors.LogNorm(vmin=.1,vmax=1e3)) divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$\Sigma{\rm gas}$ $(\rm M\odot\,pc^{-2})$",pad=0,cax=cax) ax.setaspect('equal') ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.tightlayout(hpad=1) plt.show() ```

Slices

Now let's look at the 3D gas density in a slice through the galaxy, using the Slice method. This will reconstruct the data to a grid of points in a plane slicing through the data. You can chose the order of the reconstruction: 0 will simply give the value of the nearest particle (i.e. reflecting the Voronoi domains), 1 will perform a linear reconstruction from that particle, etc. The best order will depend upon the nature of the data: smooth data will look best with higher-order reconstructions, while messier data will have nasty overshoots and artifacts. Here the density field is quite poorly-resolved, so a nearest-neighbor reconstruction is probably most appropriate.

```python fig, ax = plt.subplots(1,2,figsize=(6,6), sharex=True, sharey=True, gridspeckw={'widthratios': [1,1]}) densityslicenHcgs = M.Slice(rho,center=np.array([0,0,0]),size=40.,res=res, order=0) * RHOTONH

alternative to try: default linear reconstruction of log rho to avoid overshoot to negative density

M = Meshoid(pos, mass, hsml)

densityslicenHcgs2 = 10**M.Slice(np.log10(rho),center=np.array([0,0,0]),size=40.,res=res) * RHOTONH p = ax[0].pcolormesh(X, Y, densityslicenHcgs,norm=colors.LogNorm(vmin=.01,vmax=1e2)) p2 = ax[1].pcolormesh(X, Y, densityslicenHcgs2,norm=colors.LogNorm(vmin=.01,vmax=1e2))

setaxisargs = {"xlabel": "X (kpc)", "ylabel": "Y (kpc)","aspect": "equal"} for a in ax: a.set(**setaxisargs) divider = makeaxeslocatable(a) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$n{\rm H}$ $(\rm cm^{-3})$",cax=cax) ax[0].settitle(r"Nearest-neighbor $\rho$") ax[1].settitle(r"1st-order reconstruction of $\log \rho$") fig.tight_layout() plt.show() ```

Simple Radiative Transfer

Meshoid is also capable of performing radiative transfer with a known emissivity/source function and opacity, neglecting scattering. This is possible because there is a direct analogy between doing a line-integral through a density field (as in the project above), and solving the time-independent radiative transfer equation.

For example, we can load in the stellar positions and assume a simple constant mass-to-light ratio, and calculate the dust-extincted starlight in the V-band.

```python from meshoid.radiation import radtransfer, dustabsopacity from astropy import units as u, constants as c

kappadustcodeunits = dustabsopacity(0.555).to(u.kpc*2/(1e10c.Msun)).value # dust opacity in cgs converted to solar - evaluated at 555nm kappagas = np.repeat(kappadustcodeunits,len(mass)) jgas = np.zeroslike(mass) # assume dust does not emit

have to get the star properties now

xstars = (loadfromsnapshot("Coordinates", 4, ".", 600) - center) @ coordinatebasis mstars = loadfromsnapshot("Masses", 4, ".", 600) hstar = np.repeat(0.1, len(mstars)) # 100pc radii MASSTOLIGHTSOLAR = 1. # emissivity is just the light-to-mass ratio for stars - here assume 1 (old-ish stellar population in V-band) jstar = np.repeat(1e10/(MASSTOLIGHTSOLAR), len(mstars)) # we are assuming a constant emissivity throughout the kernel-smoothed star particles kappastars = np.zeros(len(m_stars))

now combine all emissivities, opacities, masses, kernel lengths

jall = np.atleast2d(np.concatenate([jgas, jstar])).T # 2D because this has shape (numparticles, numbands) (can have an arbitrary number of bands) kappaall = np.atleast2d(np.concatenate([kappagas, kappastars])).T # ditto kappaall = kappaall.clip(1e-100) # we divide by kappa at a certain point so put this to avoid floating-point errors hall = np.concatenate([hsml, hstar]) mall = np.concatenate([mass, mstars]) xall = np.concatenate([pos,xstars],axis=0)

rmax = 10 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') I = radtransfer(jall, mall,kappaall, xall,h_all,res,2*rmax) # actual call to rad transfer solver

screw you I'm not converting this to mag/arcsec^2

fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X,Y, I[:,:,0], norm=colors.LogNorm(vmin=I.max()/1e3, vmax=I.max()), cmap="Greysr") ax.setaspect('equal') divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$I\,\left(L\odot\,\rm kpc^{-2}\,\rm sr^{-1}\right)\, $",cax=cax) ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.show()

```

Differential Operators

Now let's play around with Meshoid's numerical differentiation. Meshoid can take both first (Meshoid.D) and second derivatives (Meshoid.D2) on unstructured data, using a kernel-weighted (or unweighted) least-squares gradient estimator.

As a first sanity check, we can try differentiating the coordinate functions, with respect to those coordinates. That ought to return an identity matrix. Note that you can differentiate scalars, vectors, or even arbitrary tensors that are defined on the meshoid. In general, differentiating a tensor of rank N will return a tensor of rank N+1.

The first time a given differentiation method is called, Meshoid can take a minute to compute the weights that it needs. Hang in there, Meshoid is working diligently and will re-use those weights the next time you need a derivative!

python M.D(pos)

array([[[ 1.00000000e+00, -1.56125113e-17,  1.56125113e-17],
        [ 2.25080371e-16,  1.00000000e+00, -2.12069945e-16],
        [ 1.30104261e-16, -3.03576608e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -5.96514483e-17, -1.77402580e-17],
        [-5.62429877e-17,  1.00000000e+00,  8.52725009e-17],
        [ 2.69437792e-16, -8.81591892e-17,  1.00000000e+00]],

       [[ 1.00000000e+00,  4.98157017e-17,  4.32874494e-16],
        [ 4.84079329e-17,  1.00000000e+00, -7.06018902e-17],
        [ 4.68917440e-18,  4.95344868e-17,  1.00000000e+00]],

       ...,

       [[ 1.00000000e+00, -1.60970141e-17, -2.35271871e-17],
        [-3.77099068e-17,  1.00000000e+00,  3.52535113e-17],
        [ 6.46048970e-17, -2.99888627e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -2.06011965e-16,  1.01874347e-16],
        [-2.15403867e-16,  1.00000000e+00,  9.90893023e-17],
        [ 1.01969214e-16,  5.41396355e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -7.45388994e-19, -2.98697699e-17],
        [-2.60208521e-16,  1.00000000e+00, -1.51354623e-16],
        [ 6.80011603e-16, -1.39645240e-16,  1.00000000e+00]]],
      shape=(2378703, 3, 3))

OK now let's look at something physical. Let's calculate the enstrophy, which is just the norm squared of the velocity gradient, and plot it as a slice.

python gradv = M.D(v) enstrophy = np.sum(gradv*gradv, axis=(1,2)) enstrophy_projection = M.Slice(enstrophy,center=np.array([0,0,0]),size=40.,res=res,order=0) fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X, Y, enstrophy_projection*.979**2, norm=colors.LogNorm(vmin=10,vmax=1e7)) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"Enstrophy $(\rm Gyr^{-2})$",cax=cax) ax.set_aspect('equal') ax.set_xlabel("X (kpc)") ax.set_ylabel("Y (kpc)") plt.show()

Walkthrough

First let's import pylab, Meshoid, and the loadfromsnapshot function for loading GIZMO outputs.

python from matplotlib import pyplot as plt import numpy as np from meshoid import Meshoid, load_from_snapshot %matplotlib inline

Now let's load some of the gas data fields from a snapshot in the public FIRE data release using loadfromsnapshot. In this case we'll perform a density cut at n_H ~ .1 cm^-3 to narrow it down to just the ISM.

```python from os import system from os.path import isfile

download the data - feel free to use your own snapshots!

snapnum = 600 for i in range(4): if not isfile(f"./snapshot{snapnum}.{i}.hdf5"): system(f"wget -r -nH --cut-dirs=6 --no-parent --reject='index.html*' -e robots=off https://users.flatironinstitute.org/~mgrudic/fire2publicrelease/core/m12ires7100/output/snapdir{snapnum}/snapshot{snapnum}.{i}.hdf5")

rho = loadfromsnapshot("Density", 0, ".", snapnum) RHOTONH = 300 # convert to H number density (approx) densitycut = (rho*RHOTONH > .01) pdata = {} fieldstoload = "Masses", "Coordinates", "SmoothingLength", "Velocities", "Density" pdata = {field: loadfromsnapshot(field, 0, ".", snapnum, particlemask=np.arange(len(rho))[densitycut]) for field in fieldsto_load} ```

Finally, before getting to the meshoid stuff we will also center the coordinates, perform a cut in galactocentric radius at 40kpc, and orient our coordinates to the principal axes of the gas distribution.

```python

do centering

pos = pdata["Coordinates"] center = np.median(pos,axis=0) pos -= center

do radius cut

MAXRADIUS = 40. radiuscut = np.sum(pos*pos,axis=1) < MAXRADIUS * MAXRADIUS pos, mass, hsml, v, rho = pos[radiuscut], pdata["Masses"][radius_cut], pdata["SmoothingLength"][radius_cut], pdata["Velocities"][radius_cut], pdata["Density"][radius_cut] centernew = np.average(pos,axis=0,weights=mass) # now do another re-centering center += centernew pos -= centernew

now get the principal axes - eigenvectors of the second mass moment matrix

covpos = np.cov(pos.T, aweights=mass) w, coordinatebasis = np.linalg.eigh(covpos) coordinatebasis = coordinatebasis[:,w.argsort()[::-1]] # sort so the smallest moment axis is the last = z-axis pos = pos @ coordinatebasis # dot product with each basis vector to get coordinates in new basis ```

Projections

OK, now let's start by making a map of gas surface density. We can do so by generating a Meshoid object from the particle masses, coordinates, and smoothing lengths, and then calling the SurfaceDensity method. This function is useful for giving kernel-weighted projected quantities on a Cartesion grid of sighlines.

Meshoid can also adaptively compute particle smoothing lengths on-the-fly provided only the coordinates, but it requires a nearest-neighbor search that takes a while so it's best to provide the smoothing length when you can.

```python import matplotlib.colors as colors from mpltoolkits.axesgrid1 import makeaxeslocatable

M = Meshoid(pos, mass, hsml) rmax = 20 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') fig, ax = plt.subplots(figsize=(6,6)) sigmagasmsunpc2 = M.SurfaceDensity(M.m,center=np.array([0,0,0]),size=40.,res=res)*1e4 p = ax.pcolormesh(X, Y, sigmagasmsunpc2, norm=colors.LogNorm(vmin=.1,vmax=1e3)) divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$\Sigma{\rm gas}$ $(\rm M\odot\,pc^{-2})$",pad=0,cax=cax) ax.setaspect('equal') ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.tightlayout(hpad=1) plt.show() ```

Slices

Now let's look at the 3D gas density in a slice through the galaxy, using the Slice method. This will reconstruct the data to a grid of points in a plane slicing through the data. You can chose the order of the reconstruction: 0 will simply give the value of the nearest particle (i.e. reflecting the Voronoi domains), 1 will perform a linear reconstruction from that particle, etc. The best order will depend upon the nature of the data: smooth data will look best with higher-order reconstructions, while messier data will have nasty overshoots and artifacts. Here the density field is quite poorly-resolved, so a nearest-neighbor reconstruction is probably most appropriate.

```python fig, ax = plt.subplots(1,2,figsize=(6,6), sharex=True, sharey=True, gridspeckw={'widthratios': [1,1]}) densityslicenHcgs = M.Slice(rho,center=np.array([0,0,0]),size=40.,res=res, order=0) * RHOTONH

alternative to try: default linear reconstruction of log rho to avoid overshoot to negative density

M = Meshoid(pos, mass, hsml)

densityslicenHcgs2 = 10**M.Slice(np.log10(rho),center=np.array([0,0,0]),size=40.,res=res) * RHOTONH p = ax[0].pcolormesh(X, Y, densityslicenHcgs,norm=colors.LogNorm(vmin=.01,vmax=1e2)) p2 = ax[1].pcolormesh(X, Y, densityslicenHcgs2,norm=colors.LogNorm(vmin=.01,vmax=1e2))

setaxisargs = {"xlabel": "X (kpc)", "ylabel": "Y (kpc)","aspect": "equal"} for a in ax: a.set(**setaxisargs) divider = makeaxeslocatable(a) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$n{\rm H}$ $(\rm cm^{-3})$",cax=cax) ax[0].settitle(r"Nearest-neighbor $\rho$") ax[1].settitle(r"1st-order reconstruction of $\log \rho$") fig.tight_layout() plt.show() ```

Simple Radiative Transfer

Meshoid is also capable of performing radiative transfer with a known emissivity/source function and opacity, neglecting scattering. This is possible because there is a direct analogy between doing a line-integral through a density field (as in the project above), and solving the time-independent radiative transfer equation.

For example, we can load in the stellar positions and assume a simple constant mass-to-light ratio, and calculate the dust-extincted starlight in the V-band.

```python from meshoid.radiation import radtransfer, dustabsopacity from astropy import units as u, constants as c

kappadustcodeunits = dustabsopacity(0.555).to(u.kpc*2/(1e10c.Msun)).value # dust opacity in cgs converted to solar - evaluated at 555nm kappagas = np.repeat(kappadustcodeunits,len(mass)) jgas = np.zeroslike(mass) # assume dust does not emit

have to get the star properties now

xstars = (loadfromsnapshot("Coordinates", 4, ".", 600) - center) @ coordinatebasis mstars = loadfromsnapshot("Masses", 4, ".", 600) hstar = np.repeat(0.1, len(mstars)) # 100pc radii MASSTOLIGHTSOLAR = 1. # emissivity is just the light-to-mass ratio for stars - here assume 1 (old-ish stellar population in V-band) jstar = np.repeat(1e10/(MASSTOLIGHTSOLAR), len(mstars)) # we are assuming a constant emissivity throughout the kernel-smoothed star particles kappastars = np.zeros(len(m_stars))

now combine all emissivities, opacities, masses, kernel lengths

jall = np.atleast2d(np.concatenate([jgas, jstar])).T # 2D because this has shape (numparticles, numbands) (can have an arbitrary number of bands) kappaall = np.atleast2d(np.concatenate([kappagas, kappastars])).T # ditto kappaall = kappaall.clip(1e-100) # we divide by kappa at a certain point so put this to avoid floating-point errors hall = np.concatenate([hsml, hstar]) mall = np.concatenate([mass, mstars]) xall = np.concatenate([pos,xstars],axis=0)

rmax = 10 res = 1024 X = Y = np.linspace(-rmax, rmax, res) X, Y = np.meshgrid(X, Y, indexing='ij') I = radtransfer(jall, mall,kappaall, xall,h_all,res,2*rmax) # actual call to rad transfer solver

screw you I'm not converting this to mag/arcsec^2

fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X,Y, I[:,:,0], norm=colors.LogNorm(vmin=I.max()/1e3, vmax=I.max()), cmap="Greysr") ax.setaspect('equal') divider = makeaxeslocatable(ax) cax = divider.appendaxes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"$I\,\left(L\odot\,\rm kpc^{-2}\,\rm sr^{-1}\right)\, $",cax=cax) ax.setxlabel("X (kpc)") ax.setylabel("Y (kpc)") plt.show()

```

Differential Operators

Now let's play around with Meshoid's numerical differentiation. Meshoid can take both first (Meshoid.D) and second derivatives (Meshoid.D2) on unstructured data, using a kernel-weighted (or unweighted) least-squares gradient estimator.

As a first sanity check, we can try differentiating the coordinate functions, with respect to those coordinates. That ought to return an identity matrix. Note that you can differentiate scalars, vectors, or even arbitrary tensors that are defined on the meshoid. In general, differentiating a tensor of rank N will return a tensor of rank N+1.

The first time a given differentiation method is called, Meshoid can take a minute to compute the weights that it needs. Hang in there, Meshoid is working diligently and will re-use those weights the next time you need a derivative!

python M.D(pos)

array([[[ 1.00000000e+00, -1.56125113e-17,  1.56125113e-17],
        [ 2.25080371e-16,  1.00000000e+00, -2.12069945e-16],
        [ 1.30104261e-16, -3.03576608e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -5.96514483e-17, -1.77402580e-17],
        [-5.62429877e-17,  1.00000000e+00,  8.52725009e-17],
        [ 2.69437792e-16, -8.81591892e-17,  1.00000000e+00]],

       [[ 1.00000000e+00,  4.98157017e-17,  4.32874494e-16],
        [ 4.84079329e-17,  1.00000000e+00, -7.06018902e-17],
        [ 4.68917440e-18,  4.95344868e-17,  1.00000000e+00]],

       ...,

       [[ 1.00000000e+00, -1.60970141e-17, -2.35271871e-17],
        [-3.77099068e-17,  1.00000000e+00,  3.52535113e-17],
        [ 6.46048970e-17, -2.99888627e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -2.06011965e-16,  1.01874347e-16],
        [-2.15403867e-16,  1.00000000e+00,  9.90893023e-17],
        [ 1.01969214e-16,  5.41396355e-16,  1.00000000e+00]],

       [[ 1.00000000e+00, -7.45388994e-19, -2.98697699e-17],
        [-2.60208521e-16,  1.00000000e+00, -1.51354623e-16],
        [ 6.80011603e-16, -1.39645240e-16,  1.00000000e+00]]],
      shape=(2378703, 3, 3))

OK now let's look at something physical. Let's calculate the enstrophy, which is just the norm squared of the velocity gradient, and plot it as a slice.

python gradv = M.D(v) enstrophy = np.sum(gradv*gradv, axis=(1,2)) enstrophy_projection = M.Slice(enstrophy,center=np.array([0,0,0]),size=40.,res=res,order=0) fig, ax = plt.subplots(figsize=(6,6)) p = ax.pcolormesh(X, Y, enstrophy_projection*.979**2, norm=colors.LogNorm(vmin=10,vmax=1e7)) divider = make_axes_locatable(ax) cax = divider.append_axes("right", size="5%", pad=0.0) fig.colorbar(p,label=r"Enstrophy $(\rm Gyr^{-2})$",cax=cax) ax.set_aspect('equal') ax.set_xlabel("X (kpc)") ax.set_ylabel("Y (kpc)") plt.show()