dmentipl / plonk Goto Github PK
View Code? Open in Web Editor NEWSmoothed particle hydrodynamics analysis and visualization with Python.
License: MIT License
Smoothed particle hydrodynamics analysis and visualization with Python.
License: MIT License
Phantom has integer codes ieos
that determine the equation of state (eos). These codes are written to the snapshot file header. These are required to calculate the temperature or pressure as these variables are not written to file. The meaning of the ieos
codes can be inferred from the source code eos.F90
.
Currently, Plonk supports:
ieos=1
: globally isothermal eos,ieos=2
: adiabatic eos,ieos=3
: vertically isothermal eos (appropriate for discs).We should add support for more Phantom equations of state.
For reference, here is the list of values at the current version:
! This module contains stuff to do with the equation of state
! Current options:
! 1 = isothermal eos
! 2 = adiabatic/polytropic eos
! 3 = eos for a locally isothermal disc as in Lodato & Pringle (2007)
! 4 = GR isothermal
! 6 = eos for a locally isothermal disc as in Lodato & Pringle (2007),
! centered on a sink particle
! 7 = z-dependent locally isothermal eos
! 8 = Barotropic eos
! 9 = Piecewise polytrope
! 10 = MESA EoS
! 11 = isothermal eos with zero pressure
! 12 = ideal gas with radiation pressure
! 14 = locally isothermal prescription from Farris et al. (2014) for binary system
! 15 = Helmholtz free energy eos
! 16 = Shen eos
! 19 = Variable gamma (requires KROME)
I tried different things in order to plot the dust density map in case of a one-fluid dump, without success.
Steps to reproduce the behaviour:
filename = 'mdisc010racc01_00050.h5'
snap = plonk.load_snap(filename)
2a. Gas visualisation (works fine):
viz = plonk.visualize.plot( snap=snap,
quantity='density',
extent=(-250, 250, -250, 250),
cmap='gist_heat')
2b. Dust visualisation (doesn't work):
viz = plonk.visualize.plot( snap=snap,
quantity='dust_density',
extent=(-250, 250, -250, 250),
cmap='gist_heat')
snap.available_arrays()
I was expecting to be able to plot it. Instead, I got this error:ValueError: 2d quantity must be a vector quantity, e.g. "velocity".
For dust quantities, try appending the dust species number,
e.g. "dust_density_001".
I tried also using dust_density_001, without having success.
I checked and apparently density and dust_density are respectively a line array and a column array. Thinking this could be the problem I converted the dust_density into a "line array", but also this doesn't work.
I also tried defining a new snap quantity:
snap['dust_dens']=snap['density']*snap['dust_fraction']
but it is incredibly slow and my notebook crashes.
I would expect to be able to plot the dust density map like the one plotted for the gas, when visualising a one-fluid dump.
Implementing the ability to produce slices though the SPH data across any arbitrary plane. Right now Plonk can only do this for the xy-plane, giving slices with constant z. When this is rotated, the slice stays in the xy plane, even the plane of the disc is no longer in the xy plane (but originally was).
The ability to be able to produce cross-sectional slices across any arbitrary plane, or possibly even a 3D surface.
e.g. being able to slice across a plane inclined with respect to the xy plane.
More advanced:
Being able to define a 3D surface to, for example, take into account the flaring of a disc, so that we can see how quantities change along a constant disc scale height.
Hi,
I'm trying to visualise my phantom simulation but I keep getting the error TypeError: expected dtype object, got 'numpy.dtype[float64]'
I've tried a few different methods of visualising the sim but the TypeError
persists.
Here is a minimum working example using the minor and major phantom dumps.
import plonk
import matplotlib.pyplot as plt
sim = plonk.load_simulation(prefix='disc', directory='hd100453_gas_disc_10mj_20au_low_acc/')
snaps = sim.snaps
units = {'position': 'au', 'density': 'g/cm^3'}
plonk.visualize.animation_images(filename='plots/hd100453_gas_disc_10mj_20au_low_acc.mp4',snaps=snaps,
units=units,
quantity='density',
save_kwargs={'fps': 10, 'dpi': 300},)
This produces the following error.
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
/Volumes/phd_backup/phantom_dumps/animate_models.py in <module>
16 'projection': 'cm'
17 }
---> 18 plonk.animate(
19 filename='animation.mp4',
20 snaps=snaps,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/animation.py in animate(filename, snaps, profiles, quantity, x, y, fig, adaptive_colorbar, adaptive_limits, text, text_kwargs, func_animation_kwargs, save_kwargs, **kwargs)
149 )
150 else:
--> 151 anim = animation_images(
152 filename=filename,
153 snaps=snaps,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/animation.py in animation_images(filename, snaps, quantity, fig, adaptive_colorbar, text, text_kwargs, func_animation_kwargs, save_kwargs, **kwargs)
247
248 if fig is None:
--> 249 ax = snaps[0].image(quantity=quantity, **kwargs)
250 fig = ax.figure
251 image = ax.images[0]
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/visualization.py in image(snap, quantity, x, y, interp, weighted, slice_normal, slice_offset, extent, units, ax, ax_kwargs, colorbar_kwargs, **kwargs)
145 """
146 with snap.context(cache=True):
--> 147 return _interpolation_plot(
148 snap=snap,
149 quantity=quantity,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/visualization.py in _interpolation_plot(snap, quantity, x, y, kind, interp, weighted, slice_normal, slice_offset, extent, units, ax, ax_kwargs, colorbar_kwargs, **kwargs)
339 # Interpolate data to plot
340 num_pixels = _kwargs.pop('num_pixels', None)
--> 341 _data, _extent, _units = _interpolated_data(
342 snap=snap,
343 quantity=quantity,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/visualization.py in _interpolated_data(snap, quantity, x, y, interp, weighted, slice_normal, slice_offset, extent, units, num_pixels)
400
401 # Interpolate in code units
--> 402 interpolated_data = interpolate(
403 snap=snap,
404 quantity=quantity,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/interpolation.py in interpolate(snap, quantity, x, y, interp, weighted, slice_normal, slice_offset, extent, num_pixels)
125
126 if _quantity.ndim == 1:
--> 127 interpolated_data = scalar_interpolation(
128 quantity=_quantity,
129 x_coordinate=x,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/interpolation.py in scalar_interpolation(quantity, x_coordinate, y_coordinate, dist_from_slice, extent, smoothing_length, particle_mass, hfact, weighted, num_pixels)
208 shape (npixx, npixy).
209 """
--> 210 return _interpolate(
211 quantity=quantity,
212 x_coordinate=x_coordinate,
/anaconda3/envs/py38/lib/python3.8/site-packages/plonk/visualize/interpolation.py in _interpolate(quantity, x_coordinate, y_coordinate, dist_from_slice, extent, smoothing_length, particle_mass, hfact, weighted, num_pixels)
352 )
353 else:
--> 354 interpolated_data = interpolate_projection(
355 x=x_coordinate,
356 y=y_coordinate,
TypeError: expected dtype object, got 'numpy.dtype[float64]'
I can share the dumps if you need.
Cheers,
Brodie
The documentation covers loading, and making basic renders, but requires some API exploration to determine how to use render
to, for instance, visualize along a different axis.
Having in the docs examples that show rotating the snap, then render
-ing the snap, would make this process more obvious.
I am using Plonk to make a profile plot and I have a problem with that. After reading the list of profiles and appending them, only the first snap is plotted.
plonk version: 0.7.3
python: 3.7.4
Linux mint 18.1
jupyter lab
import plonk
import matplotlib.pyplot as plt
import numpy as np
filename = 'file_address'
snap = plonk.load_snap(filename)
snap = plonk.load_snap(filename)
snap.add_quantities('disc')
prof = plonk.load_profile(snap, cmin='10 au', cmax='200 au')
snap['timestep'].to('Myr')
snap['x'].to('pc')
sim = plonk.load_simulation(prefix='disc')
sinkFilePath = !ls 'my sink_file_path'
times = sim.properties['time'].to('Myr')
for sinklist in sinkFilePath:
print(sinklist)
profiles = list()
snap = plonk.load_snap(sinklist)
profile = plonk.load_profile(snap, cmin='10 au', cmax='150 au', n_bins=50)
profiles.append(profile)
fig, ax = plt.subplots()
units = {'position': 'pc', 'density': 'g/cm^3'}
for time, profile in zip(times, profiles):
label = f'{time.m:.0f}'
profile.plot(
'radius', 'density', units=units, label=label, ax=ax
)
ax.set_ylabel('Density [g/cm${}^3$]')
ax.legend(title='Time [Myr]', loc='best')
I also need to plot multiple snapshots in a single panel. I used to do it with splash and I want to know if it is possible in plonk. I expect a plot similar to what splash produces.
multiple_plot.pdf
The plonk.animation
function no longer works.
import plonk
sim = plonk.load_sim('disc')
plonk.animation(
filename='density.mp4',
snaps=sim.snaps,
quantity='density',
)
This produces the following error.
---------------------------------------------------------------------------
TypingError Traceback (most recent call last)
<ipython-input-3-d7a778b00f2e> in <module>
----> 1 plonk.animation(
2 filename='density.mp4',
3 snaps=sim.snaps,
4 quantity='density',
5 )
~/repos/plonk/plonk/visualize/animation.py in animation(filename, snaps, quantity, x, y, extent, fig, adaptive_colorbar, text, text_kwargs, func_animation_kwargs, save_kwargs, **kwargs)
204 fig, animate, frames=len(snaps), **func_animation_kwargs
205 )
--> 206 anim.save(filepath, extra_args=['-vcodec', 'libx264'], **save_kwargs)
207 plt.close()
208 if tqdm is not None:
~/conda/lib/python3.8/site-packages/matplotlib/animation.py in save(self, filename, writer, fps, dpi, codec, bitrate, extra_args, metadata, extra_anim, savefig_kwargs, progress_callback)
1118 _is_saving=True, manager=None):
1119 for anim in all_anim:
-> 1120 anim._init_draw() # Clear the initial frame
1121 frame_number = 0
1122 # TODO: Currently only FuncAnimation has a save_count
~/conda/lib/python3.8/site-packages/matplotlib/animation.py in _init_draw(self)
1693 # artists.
1694 if self._init_func is None:
-> 1695 self._draw_frame(next(self.new_frame_seq()))
1696
1697 else:
~/conda/lib/python3.8/site-packages/matplotlib/animation.py in _draw_frame(self, framedata)
1716 # Call the func with framedata and args. If blitting is desired,
1717 # func needs to return a sequence of any artists that were modified.
-> 1718 self._drawn_artists = self._func(framedata, *self._args)
1719 if self._blit:
1720 if self._drawn_artists is None:
~/repos/plonk/plonk/visualize/animation.py in animate(idx)
179 else:
180 _extent = extent
--> 181 interp_data = interpolate(
182 snap=snaps[idx],
183 interp=interp,
~/repos/plonk/plonk/visualize/interpolation.py in interpolate(snap, quantity, x, y, interp, z_slice, extent, **kwargs)
91
92 if _quantity.ndim == 1:
---> 93 interpolated_data = scalar_interpolation(
94 quantity=_quantity,
95 x_coordinate=x,
~/repos/plonk/plonk/visualize/interpolation.py in scalar_interpolation(quantity, x_coordinate, y_coordinate, z_coordinate, extent, smoothing_length, particle_mass, hfact, number_of_pixels, cross_section, density_weighted)
175 shape (npixx, npixy).
176 """
--> 177 return _interpolate(
178 quantity=quantity,
179 x_coordinate=x_coordinate,
~/repos/plonk/plonk/visualize/interpolation.py in _interpolate(quantity, x_coordinate, y_coordinate, z_coordinate, extent, smoothing_length, particle_mass, hfact, number_of_pixels, cross_section, density_weighted)
331 )
332 else:
--> 333 interpolated_data = interpolate_projection(
334 x=x_coordinate,
335 y=y_coordinate,
~/conda/lib/python3.8/site-packages/numba/core/dispatcher.py in _compile_for_args(self, *args, **kws)
413 e.patch_message(msg)
414
--> 415 error_rewrite(e, 'typing')
416 except errors.UnsupportedError as e:
417 # Something unsupported is present in the user code, add help info
~/conda/lib/python3.8/site-packages/numba/core/dispatcher.py in error_rewrite(e, issue_type)
356 raise e
357 else:
--> 358 reraise(type(e), e, None)
359
360 argtypes = []
~/conda/lib/python3.8/site-packages/numba/core/utils.py in reraise(tp, value, tb)
78 value = tp()
79 if value.__traceback__ is not tb:
---> 80 raise value.with_traceback(tb)
81 raise value
82
TypingError: Failed in nopython mode pipeline (step: nopython frontend)
non-precise type pyobject
During: typing of argument at /Users/daniel/repos/plonk/plonk/visualize/splash.py (157)
File "../../../repos/plonk/plonk/visualize/splash.py", line 157:
def interpolate_projection(
<source elided>
"""
coltable = setup_integratedkernel()
^
This error may have been caused by the following argument(s):
- argument 7: cannot determine Numba type of <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>
- argument 8: cannot determine Numba type of <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>
- argument 11: cannot determine Numba type of <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>
- argument 12: cannot determine Numba type of <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>
I expect it to produce a file "density.mp4" with an animation of the rendered density.
I used the dataset from the documentation (https://anaconda.org/dmentipl/plonk_example_data/).
I think the issue is in returning non-dimensional data, and needing to add units.
interp_data = interpolate(
snap=snaps[idx],
interp=interp,
quantity=quantity,
x=x,
y=y,
extent=_extent,
**_kwargs,
)
Currently, the low-level SPH interpolation routines live in splash.py
. These are rewritten from the original Splash Fortran routines using NumPy and Numba for performance.
The package sph-interp contains very similar functions. Plonk should use this package instead of splash.py
. Given that these functions are quite "low-level" it makes sense that they live outside of the Plonk repository. Then other packages can use them.
This function is to replicate the functionality of Phantom's moddump.
A new module/class/function to modify snapshots and write new ones to file.
In line 329 of profile.py, the definition of velocity reads:
vel = self.Snap['vxyz'][self._mask]
and it should be
vel = self.snap['vxyz'][self._mask]
after loading a snapshot, running the following code results in an error
prof = plonk.analysis.Profile(snap, radius_min=0, radius_max=6)
ecc = prof._profile_functions["eccentricity"](prof, gravitational_parameter = 2)
return radial profile of eccentricity
When using the visualization.MultiPlot class, the colorbars are not displayed properly. There are two issues.
I am plotting a Phantom dump obtained with a one fluid multigrain simulation.
This is part of the code I'm using for the multiplot:
totdustfrac = np.zeros(1000000)
for s in range(0,5):
grains = dump.header['grainsize'][s]*dump.header['udist']*10**(4)
dustfrac = dump.particles.arrays['dustfrac'][:,s]
totdustfrac +=dump.particles.arrays['dustfrac'][:,s]
new_units = plonk.Units(length='au', mass='g', time='yr', density=1., velocity=1.0e5)
density = dump.density
density = dump.units.convert_quantity_to_new_units(density, 'density', new_units)
gas_density = density*(1.-totdustfrac)
options = list()
options.append({'render': gas_density,
'render_scale': 'log',
'extent': [-100,100,-100,100],
'inclination': -65,
'position_angle': -21,
'normalize':False,
'render_label': r'$\Sigma_{\rm gas}$ [g cm${}^{-2}$]',
'integrated_z':plonk.constants.au,
})
label_dust = [r'$\Sigma_{{\rm d},1.77\,\mu{\rm m}}$ [g cm${}^{-2}$]',
r'$\Sigma_{{\rm d},5.56\,\mu{\rm m}}$ [g cm${}^{-2}$]',
r'$\Sigma_{{\rm d},0.017\,{\rm mm}}$ [g cm${}^{-2}$]',
r'$\Sigma_{{\rm d},0.054\,{\rm mm}}$ [g cm${}^{-2}$]',
r'$\Sigma_{{\rm d},0.17\,{\rm mm}}$ [g cm${}^{-2}$]',
r'$\Sigma_{{\rm d},0.54\,{\rm mm}}$ [g cm${}^{-2}$]',
r'$\Sigma_{{\rm d},1.7 {\rm mm}}$ [g cm${}^{-2}$]']
for s in range(0,5):
options.append({'render': density*dump.particles.arrays['dustfrac'][:,s],
'render_scale': 'log',
'extent': [-100,100,-100,100],
'inclination': -65,
'position_angle': -21,
'render_label': label_dust[s],
'font_size': 30,
'normalize': False,
'integrated_z':plonk.constants.au,
})
multiplot = plonk.visualization.MultiPlot(
dump,
options,
scale=7,
shape=(3,2),
axes_pad=0.7,
cbar_mode='each',
cbar_location='top',
cbar_pad=0.2
)
NxM subplots, each with a colorbar in log scale.
Here you can see the (1x6) shape image, displaying all the colorbars, but not in log scale (the color rendering is performed properly in log). The (6x1) shape is behaving at the same way.
Here you can see the (2x3) shape image, not displaying properly the last colorbars. The colorbar is plotted only for the 0, the 3rd, 4th and 5th element.
Here you can see the (3x2) shape image, not displaying properly the last colorbars. The colorbar is plotted only for the 0, the 2nd, 4th and 5th element.
Split arrays on SubSnap objects return whole array. I.e. rather than getting only 'x' array you get the whole 'position' array.
Steps to reproduce the behaviour:
import plonk
snap = plonk.load_snap(path_to_snap_file)
subsnap = snap[snap['id'] < 10]
x = subsnap['x']
x.shape
# prints (n_particles, 3)
subsnap['x']
should just return the 'x' array.
Should the 'id' array be 'particle_type'? Then 'id' could be an integer that uniquely identifies a particle. This could be useful to track particles through simulations.
Accessing a particle array on a snapshot is currently lazy, in the sense that the data is only loaded from disc into memory when requested, e.g. with snap['position']
. However, it remains there in the dictionary snap._arrays
.
An alternative is to load the array as a Dask array and use the resulting object as if it were a NumPy array loaded into memory. Then complicated expressions can be written before anything is loaded into memory. The computation is executed with the .compute()
method. Dask also allows for easy parallelization on both local, i.e. a multi-core laptop, and remote hardware, i.e. a supercomputing cluster.
A suggestion is to adjust the __getitem__()
method to return a Dask array, and to not store the array in snap._arrays
.
Column density plots in plonk don't appear to show the same information as in splash? I plotted a broken disc and it looked fuzzier than I thought it should and despite adjustments it didn't improve. Have I missed something? It appears that high density structure in the inner region and the gap between the inner and outer disc are not accurately reproduced in plonk? See the examples attached below.
density_plonk.pdf
density_splash.pdf
Following the example in the wiki as much as possible, I plotted with:
viz = plonk.visualize.render(snap=snap,
quantity='density',
axis=axes,
extent=[-15,15,-15,15],
scalar_options={'cmap': 'gist_heat',
'norm': 'log',
'plot_colorbar': cbar}
)
I'm using:
c816ba88
, gas only simulation.Plonk does not appear to to recover the temperature correctly when used on mcfost+phantom dump, indicating temperatures of a few 1e6K for a protoplanetary disc.
See the data in the dropbox link below, and the following script
"""Calculate vertical profile in discs."""
import matplotlib.pyplot as plt
import numpy as np
import plonk
def vertical_profile(snap, radius, vertical_height, radial_width=1):
"""Vertical profile at particular radius.
Parameters
----------
snap
The Snap.
radius
The radius at which to calculate the vertical profile.
vertical_height
The height over which to compute the vertical profile.
radial_width
The radial width of the subset of particle on which to compute
the vertical profile.
"""
# Load extra quantities, e.g. cylindrical radius
snap.extra_quantities()
# Generate subsnap in a cylindrical ring
subsnap = snap[
(snap['R'] > radius - radial_width / 2)
& (snap['R'] < radius + radial_width / 2)
]
# Return Cartesian z-profile
return plonk.load_profile(
subsnap,
ndim=1,
coordinate='z',
radius_min=-vertical_height / 2,
radius_max=vertical_height / 2,
)
if __name__ == '__main__':
# Load snapshot
snap = plonk.load_snap('hydrostatic_18790.h5')
# Set molecular weight for temperature
snap.set_molecular_weight(2.381)
# Choose radii at which to calculate z-profiles.
radius = np.linspace(25, 150, 5)
# Vertical height
R0 = 10 # reference radius
H0 = 0.034 * R0 # scale height at reference
q = 1 / 4 # sound speed index
scale_height = H0 * (radius / R0) ** (3 / 2 - q)
vertical_height = 10 * scale_height
# Make figure and axes
fig, axs = plt.subplots(ncols=2, figsize=(12, 5))
axs[0].set(xlabel='altitude', ylabel='density')
axs[1].set(xlabel='altitude', ylabel='temperature')
# Loop over all radii and plot density and temperature profiles
for R, dZ in zip(radius, vertical_height):
prof = vertical_profile(snap=snap, radius=R, vertical_height=dZ)
prof.plot('radius', 'density', label=f'radius={R:.0f}', ax=axs[0])
prof.plot('radius', 'temperature', label=f'radius={R:.0f}', ax=axs[1])
axs[1].legend().remove()
plt.show()
I would love to get the correct mcfost temperature ;-)
https://www.dropbox.com/sh/goc9rukjghnca0d/AAABP1n5N-iRlAukFkwUAAmda?dl=0
Sink particles are not rotated with snap.rotate
Steps to reproduce the behaviour:
import plonk
from scipy.spatial.transform.rotation import Rotation
snap = plonk.load_snap(path_to_file)
rotation = Rotation.from_rotvec([1, 1, 0])
print(snap.sinks['position'])
snap.rotate(rotation)
print(snap.sinks['position'])
The two print statements are the same but should show that the snapshot has been rotated.
The sinks should rotate with a snap rotation.
I am trying to make a scatter plot of particles by visualize and I get the following error while I follow the instruction of "Plonk tutorial written for the 3rd Phantom + MCFOST users workshop, Monash University, Feb 24-28, 2020". I also used the sample data of the workshop to be sure that data format is not related in problem.
plonk version: 0.7.3
python: 3.7.4
Linux mint 18.1
jupyter lab
import plonk
import matplotlib.pyplot as plt
import numpy as np
filename = 'path_to_file.h5'
snap = plonk.load_snap(filename)
viz = plonk.visualize.render(snap=snap, quantity='density', extent=[-200, 200, -200, 200])
which raises this error:
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
<ipython-input-57-c81b5b94f9b3> in <module>
----> 1 viz = plonk.visualize.render(snap=snap, quantity='density', extent=[-200, 200, -200, 200])
AttributeError: module 'plonk.visualize' has no attribute 'render'
and when I use dir(plonk.visualize) to return the attributes, render is not in the list.
['__all__',
'__builtins__',
'__cached__',
'__doc__',
'__file__',
'__loader__',
'__name__',
'__package__',
'__path__',
'__spec__',
'animate',
'animation',
'animation_images',
'animation_particles',
'animation_profiles',
'image',
'interpolate',
'interpolation',
'plot',
'plots',
'simulation',
'splash',
'vector',
'visualization',
'visualize_sim']
I also have got an error when using plonk.animate:
FilePath = !ls ./disc_000*.h5
for Filelist in FilePath:
print(Filelist)
snaps = plonk.load_snap(Filelist)
print(snaps)
units = {
'position': 'au',
'density': 'g/cm^3',
'projection': 'cm'
}
plonk.animate(
filename='animation.mp4',
snaps=snaps,
quantity='density',
units=units,
save_kwargs={'fps': 10, 'dpi': 300}
)
and after reading and printing all data, this produces the following error.
---------------------------------------------------------------------------
ZeroDivisionError Traceback (most recent call last)
<ipython-input-83-349ba0d90912> in <module>
13 'projection': 'cm'
14 }
---> 15 plonk.animate(
16 filename='animation.mp4',
17 snaps=snaps,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/animation.py in animate(filename, snaps, profiles, quantity, x, y, fig, adaptive_colorbar, adaptive_limits, text, text_kwargs, func_animation_kwargs, save_kwargs, **kwargs)
149 )
150 else:
--> 151 anim = animation_images(
152 filename=filename,
153 snaps=snaps,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/animation.py in animation_images(filename, snaps, quantity, fig, adaptive_colorbar, text, text_kwargs, func_animation_kwargs, save_kwargs, **kwargs)
247
248 if fig is None:
--> 249 ax = snaps[0].image(quantity=quantity, **kwargs)
250 fig = ax.figure
251 image = ax.images[0]
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/visualization.py in image(snap, quantity, x, y, interp, weighted, slice_normal, slice_offset, extent, units, ax, ax_kwargs, colorbar_kwargs, **kwargs)
145 """
146 with snap.context(cache=True):
--> 147 return _interpolation_plot(
148 snap=snap,
149 quantity=quantity,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/visualization.py in _interpolation_plot(snap, quantity, x, y, kind, interp, weighted, slice_normal, slice_offset, extent, units, ax, ax_kwargs, colorbar_kwargs, **kwargs)
339 # Interpolate data to plot
340 num_pixels = _kwargs.pop('num_pixels', None)
--> 341 _data, _extent, _units = _interpolated_data(
342 snap=snap,
343 quantity=quantity,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/visualization.py in _interpolated_data(snap, quantity, x, y, interp, weighted, slice_normal, slice_offset, extent, units, num_pixels)
400
401 # Interpolate in code units
--> 402 interpolated_data = interpolate(
403 snap=snap,
404 quantity=quantity,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/interpolation.py in interpolate(snap, quantity, x, y, interp, weighted, slice_normal, slice_offset, extent, num_pixels)
125
126 if _quantity.ndim == 1:
--> 127 interpolated_data = scalar_interpolation(
128 quantity=_quantity,
129 x_coordinate=x,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/interpolation.py in scalar_interpolation(quantity, x_coordinate, y_coordinate, dist_from_slice, extent, smoothing_length, particle_mass, hfact, weighted, num_pixels)
208 shape (npixx, npixy).
209 """
--> 210 return _interpolate(
211 quantity=quantity,
212 x_coordinate=x_coordinate,
~/anaconda3/envs/plonk-env/lib/python3.8/site-packages/plonk/visualize/interpolation.py in _interpolate(quantity, x_coordinate, y_coordinate, dist_from_slice, extent, smoothing_length, particle_mass, hfact, weighted, num_pixels)
352 )
353 else:
--> 354 interpolated_data = interpolate_projection(
355 x=x_coordinate,
356 y=y_coordinate,
ZeroDivisionError: division by zero
Using the image(quanity='density') method to produce integrated density plots with a logarithmic colour bar causes artefact in regions of a disc where there are no particles. The minimum density colour does not extend to the edges of the box or fill in the gap around the sink
import plonk
snap = plonk.load_snap('disc_00010.h5')
units = {'position': 'au', 'density': 'g/cm^3', 'projection': 'cm'}
snap.image(
... quantity='density',
... extent=(-120,120,-120,120),
... units=units,
... norm='log',
... vmin=1,
... vmax=3000,
... )
plt.show()
No error is produced, the plot is produced successfully but does not look correct when using the norm='log' option.
Expected the plotted area to be filled with the colour corresponding to vmin
Use the attached .h5 file with plonk (code above) to reproduce problem
I need access to some properties of the simulation itself (number of particles, accretion radius of the sinks, alpha_av etc).
For now, I just did a small function in order to extract this data and have it in a dictionary, but I think it would be useful to have it inside the simulation class, or just have a plonk function to take the .setup and .in file and load them.
Automate building the Conda and pip packages with Travis CI.
Currently, the Conda Forge build recipe requires adding the SHA256 hash from PyPI by hand.
Hi, I'm attempting to shift to the disc centre for both snap images and an animation. In phantom I set a 5mj companion at 10au but the sink position in plonk is significantly far from this.
I rotate and shift the disc with:
snap_file = 'hd100453_bec_run_5mj_10au/disc_00101.h5'
snap = plonk.load_snap(snap_file)
analysis.discs.rotate_face_on(snap)
com = analysis.total.center_of_mass(snap, sinks=False)
snap.translate(-com)
Then setting the units and printing the sink positions:
snap.set_units(position='au')
print(snap.sinks['position'])
The output.
[[1.0657612036467647 -0.7482025165997441 0.08040807544826903] [151.16832435718015 -17.357874795606797 -215.90554676449221] [-14942.404548791676 35342.100422533564 18413.824209907314]] astronomical_unit
If I print the mass with
print(snap.sinks['mass'])
The output:
[3.3820768744152954e+30 3.9782982217579486e+29 9.490655967300002e+27] kilogram
Therefor the third sink is clearly the 5mj companion but the xyz coordinates are obviously not at a 10au radius.
Is this a plonk issue for a phantom issue?
As for the animation, I'm not sure how to centre the disc at all or label the sinks with a marker. Any tips?
I load the entire simulation with:
sim = plonk.load_simulation(prefix='disc', directory='hd100453_bec_run_5mj_10au/')
It appears that I'm shifting correctly to the star sink.
Here's the phantom disc.setup for reference
# input file for disc setup routine
# resolution
np = 1000000 ! number of gas particles
# units
dist_unit = au ! distance unit (e.g. au,pc,kpc,0.1pc)
mass_unit = solarm ! mass unit (e.g. solarm,jupiterm,earthm)
# central object(s)/potential
icentral = 1 ! use sink particles or external potential (0=potential,1=sinks)
nsinks = 2 ! number of sinks
ibinary = 0 ! binary orbit (0=bound,1=unbound [flyby])
# options for binary
m1 = 1.700 ! primary mass
m2 = 0.200 ! secondary mass
binary_a = 207. ! binary semi-major axis
binary_e = 0.320 ! binary eccentricity
binary_i = 49. ! i, inclination (deg)
binary_O = 47. ! Omega, PA of ascending node (deg)
binary_w = 18. ! w, argument of periapsis (deg)
binary_f = 180. ! f, initial true anomaly (deg,180=apastron)
accr1 = 5.000 ! primary accretion radius
accr2 = 5.000 ! secondary accretion radius
# options for multiple discs
use_binarydisc = F ! setup circumbinary disc
use_primarydisc = T ! setup circumprimary disc
use_secondarydisc = F ! setup circumsecondary disc
use_global_iso = F ! globally isothermal or Farris et al. (2014)
# options for circumprimary gas disc
isetgasprimary = 0 ! how to set gas density profile (0=total disc mass,1=mass within annulus,2=surface density normalisation,3=surface density at reference rad\
ius,4=minimum Toomre Q)
itapergasprimary = F ! exponentially taper the outer disc profile
ismoothgasprimary = F ! smooth inner disc
iwarpprimary = F ! warp disc
R_inprimary = 20. ! inner radius
R_refprimary = 20. ! reference radius
R_outprimary = 60. ! outer radius
disc_mprimary = 0.003 ! disc mass
pindexprimary = 1.000 ! p index
qindexprimary = 0.250 ! q index
posanglprimary = 183.5 ! position angle (deg)
inclprimary = 14.9 ! inclination (deg)
H_Rprimary = 0.050 ! H/R at R=R_ref
alphaSS = 0.005 ! desired alphaSS
# set planets
setplanets = 1 ! add planets? (0=no,1=yes)
nplanets = 1 ! number of planets
# planet:1
mplanet1 = 5.000 ! planet mass (in Jupiter mass)
rplanet1 = 10. ! planet distance from star
inclplanet1 = 14.9 ! planet orbital inclination (deg)
accrplanet1 = 0.250 ! planet accretion radius (in Hill radius)
# timestepping
norbits = 10 ! maximum number of outer planet orbits
deltat = 0.100 ! output interval as fraction of orbital period
Hi Daniel,
I hope you are well. I'm trying to track individual particles from snapshot to snapshot through the simulation and am having some issues. I thought of using the particle ID numbers but these don't seem to be allocated to the same particles in one snapshot as they are in the next.
Is there something I am missing that enables the user to do something like this? Find the ID number of a particle at x1,y1,z1 at t = t_1 and then, using that ID number, find x2,y2,z2 of the particle at t = t2?
It would be nice to have a method, say snap.physical_units()
that sets the snap to give arrays and other values in physical units rather than having to multiply by the appropriate combination of snap.properties['utime']
, snap.properties['udist']
, snap.properties['umass']
.
Implement a method on the snap like snap.physical_units()
with arguments to set the time, mass, distance units, such that when accessing arrays like snap['position']
they return an array with physical units.
This solution could return a NumPy array in the correct units. Or, perhaps a better option is to return a Pint array which contains the units along with the array.
The alternative is to continue to do this manually with the units provided in snap.properties
.
We could use unyt or Astropy units instead.
The docs have a getting started section, but in order to play along, we need the sample data. The results of each step should be provided as well so we can verify the library is working as intended.
sample dataset referenced in the docs should be included
A clear and concise description of any alternative solutions or features you've considered.
Add any other context or screenshots about the feature request here.
I am trying to change the position of the colorbar to the left hand side of the plot. I think I should be able to do this by setting colorbar_kwargs={'position': 'left'}
, but this raises the error
AttributeError: 'AxesImage' object has no property 'colorbar_kawrgs'
It seems that colorbar_kwargs
is getting passed to the plots.imshow
function, which raises this error.
The code below reproduces the error
dumpfile = 'internal/e40_q2_gtd/rescale_hires_00092.h5'
snap= plonk.load_snap(dumpfile)
RENDER = 'density'
dens = plonk.visualize.plot(snap=snap, quantity=RENDER,
show_colorbar=True, cmap='gist_heat', norm='log', colorbar_kwargs={'position': 'left'})
Also another thing I noticed in the documentation for plonk.visualize.plot
, it states it returns a matplotlib axes object:
Returns
-------
ax
The matplotlib Axes object.
However, this does not seem to be the case:
>>> print(dens)
<plonk.Visualization "rescale_hires_00092.h5">
I am also wondering if it is possible to make the plot_object
object available to the user somehow in the above object, or make it accessible some other way. It would be useful if one wants to have more control over other matplotlib features themselves. For example, if I wanted to customise the colorbar after making the plot, I need plot_object
to do so.
Would you be interested in adding support for other snapshot formats, e.g. the AMUSE HDF5 format?
Currently the only Phantom header quantities read in are 'time', 'udist', 'utime', 'umass', 'hfact', 'ieos', 'gamma', 'polyk', 'qfacdisc'.
Other quantities may be required.
This could use matplotlib animations, https://matplotlib.org/api/animation_api.html.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.