Plotting results

Rendering

In SPH we have access to rendering, it can be used in runscripts as follows :

First define the parameters of the plot

pixel_x = 1920
pixel_y = 1080
radius = 5
center = (0.0, 0.0, 0.0)

aspect = pixel_x / pixel_y
pic_range = [-radius * aspect, radius * aspect, -radius, radius]
delta_x = (radius * 2 * aspect, 0.0, 0.0)
delta_y = (0.0, radius * 2, 0.0)

You can do a column integrated plot :

arr_rho = model.render_cartesian_column_integ(
    "rho", "f64", center=(0.0, 0.0, 0.0), delta_x=delta_x, delta_y=delta_y, nx=pixel_x, ny=pixel_y
)

Or a slice :

arr_rho = model.render_cartesian_slice(
    "rho", "f64", center=(0.0, 0.0, 0.0), delta_x=delta_x, delta_y=delta_y, nx=pixel_x, ny=pixel_y
)

Note here that you can save the obtained numpy array using np.save and recover it using np.load

Matplotlib

You can then either do a standard plot like so :

import copy, matplotlib

my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat"))  # copy the default cmap
my_cmap.set_bad(color="black")

plt.figure(figsize=(16 / 2, 9 / 2))
res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range, norm="log", vmin=1e-9)

cbar = plt.colorbar(res, extend="both")
cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")
# or r"$\rho$ [code unit]" for slices

plt.title("t = {:0.3f} [code unit]".format(model.get_time()))
plt.xlabel("x")
plt.ylabel("z")

This should result in something like this :

Or you do the same plot using the splash cinematic way like so :

import copy, matplotlib

my_cmap = copy.copy(matplotlib.colormaps.get_cmap("gist_heat"))  # copy the default cmap
my_cmap.set_bad(color="black")

dpi = 200
plt.figure(dpi=dpi)
plt.gca().set_position((0, 0, 1, 1))
plt.gcf().set_size_inches(pixel_x / dpi, pixel_y / dpi)
plt.axis("off")

res = plt.imshow(arr_rho, cmap=my_cmap, origin="lower", extent=pic_range, norm="log", vmin=1e-9)

axins = plt.gca().inset_axes([0.73, 0.1, 0.25, 0.025])
cbar = plt.colorbar(res, cax=axins, orientation="horizontal", extend="both")
cbar.set_label(r"$\int \rho \, \mathrm{d} z$ [code unit]")

from matplotlib.offsetbox import AnchoredText

anchored_text = AnchoredText("t = {:0.3f} [code unit]".format(model.get_time()), loc=2)
plt.gca().add_artist(anchored_text)

Adding sink particles

You can also add sink particles to the plot like so

output_list = []
for s in model.get_sinks():
    x, y, z = s["pos"]
    output_list.append(plt.Circle((x, z), s["accretion_radius"], color="chartreuse", fill=False))
for circle in output_list:
    plt.gca().add_artist(circle)

Details of the rendering function

In this sample the core function to be used is model.render_cartesian_slice which will render the slice using the GPU trees. It will return a numpy array with the slice.

The arguments to be passed to the function are :

• field name: the field name to render
• field_type: the type of the data to be rendered
• center: the center of the slice
• delta_x and delta_y: the size of the slice
• nx and ny: the number of pixels in the slice

Column density algorithm

The standard SPH kernel is defined as being

\[ W(\mathbf{r},h) = \frac{C_{norm}}{h^3} f\left(q = \vert \frac{\mathbf{r}}{h}\vert \right) \]

If we integrate along the z axis we can define a new kernel

\[ Y(x,y,h) = \frac{C_{norm}}{h^3} \int \mathrm{d}z \, f\left(q = \vert \frac{\mathbf{r}}{h}\vert \right) \]

Which can can transform as follows

\[ Y(x,y,h) = \frac{C_{norm}}{h^2} \int \mathrm{d}q_z \, f\left(\sqrt{q_x^2 + q_y^2 + q_z^2}\right) = \frac{C_{norm}}{h^2} g\left(\sqrt{q_x^2 + q_y^2 }\right) \]

where

\[ g(q_x) = \int \mathrm{d}q_z \, f\left(\sqrt{q_x^2 + q_z^2}\right) \]

We can use a Riemann sum to integrate this function, this gives the following plot

We can see here that already with only 4 points in the riemann sum we approximate the kernel \(g\) fairly well. We can compare the error with a large number of points (1024 here)

indeed only 4 points yield a sufficient precision for the integrated kernel, which make it sufficiently cheap to be usable as is.