Basic usage¶
DustPy is using on the simframe package for running scientific simulations. For a detailed description of the usage of simframe please have a look at the Simframe Documentation.
In this notebook demonstrates how to run the most simple DustPy model, i.e., the model that is run by default, how to plot data, and how to resume simulations from dump files.
The Simulation Frame¶
To set up a model we have to import the Simulation class from the DustPy package.
from dustpy import Simulation
We can now create an instance of this class.
sim = Simulation()
sim is now an empty simulation object that controls our simulation.
sim
DustPy
------
dust : Group (Dust quantities)
gas : Group (Gas quantities)
grid : Group (Grid quantities)
star : Group (Stellar quantities)
-----
t : NoneType
-----
Integrator : not specified
Writer : not specified
All the fields are initialized with None. All attributes can be easiliy addressed via
sim.gas
Group (Gas quantities)
----------------------
boundary : Group (Boundary conditions)
S : Group (Source terms)
v : Group (Velocities)
-----
alpha : NoneType
cs : NoneType
eta : NoneType
Fi : NoneType
gamma : NoneType
Hp : NoneType
mfp : NoneType
mu : NoneType
n : NoneType
nu : NoneType
P : NoneType
rho : NoneType
Sigma : NoneType
SigmaFloor : NoneType
T : NoneType
-----
Initializing¶
We can now initialize the sim object with Simulation.initialize(). DustPy will then fill all the fields with default values.
sim.initialize()
As you can see, the sim object has now values assigned to the fields. All quantities are in cgs units.
sim.gas
Group (Gas quantities)
----------------------
boundary : Group (Boundary conditions)
S : Group (Source terms)
v : Group (Velocities)
-----
alpha : Field (Turbulent alpha parameters)
cs : Field (Sound speed [cm/s])
eta : Field (Pressure gradient parameter)
Fi : Field (Gas flux interfaces [g/cm/s])
gamma : Field (Adiabatic index)
Hp : Field (Pressure scale height [cm])
mfp : Field (Midplane mean free path [cm])
mu : Field (Mean molecular weight [g])
n : Field (Miplane number density [1/cm³])
nu : Field (Kinematic viscosity [cm²/s])
P : Field (Midplane pressure [g/cm/s²])
rho : Field (Miplane mass density [g/cm³])
Sigma : Field (Surface density [g/cm²])
SigmaFloor : Field (Floor value of surface density [g/cm²])
T : Field (Temperature [K])
-----
You can also display the full table of contents of the sim object.
sim.toc
DustPy - dust: Group (Dust quantities) - a: Field (Particle size [cm) - backreaction: Group (Backreaction coefficients) - A: Field (Pull factor) - B: Field (Push factor) - boundary: Group (Boundary conditions) - inner: Constant gradient - outer: Value - coagulation: Group (Coagulation quantities) - A: Field (Fragment normalization factors), constant - eps: Field (Remnant mass distribution), constant - lf_ind: Field (Index of largest fragment), constant - phi: Field (Fragment distribution), constant - rm_ind: Field (Smaller index of remnant), constant - stick: Field (Sticking matrix), constant - stick_ind: Field (Non-zero elements of sticking matrix), constant - D: Field (Diffusivity [cm²/s]) - delta: Group (Mixing parameters) - rad: Field (Radial mixing parameter) - turb: Field (Turbulent mixing parameter) - vert: Field (Vertical mixing parameter) - eps: Field (Dust-to-gas ratio) - Fi: Group (Fluxes) - adv: Field (Advective flux [g/cm/s) - diff: Field (Diffusive flux [g/cm/s) - tot: Field (Total flux [g/cm/s) - fill: Field (Filling factor) - H: Field (Scale heights [cm]) - kernel: Field (Collision kernel [cm²/s]) - p: Group (Probabilities) - frag: Field (Fragmentation probability) - stick: Field (Sticking probability) - rho: Field (Midplane mass density per mass bin [g/cm³]) - rhos: Field (Solid state density [g/cm³]) - S: Group (Sources) - coag: Field (Coagulation sources [g/cm²/s]) - ext: Field (External sources [g/cm²/s]) - hyd: Field (Hydrodynamic sources [g/cm²/s]) - tot: Field (Tot sources [g/cm²/s]) - Sigma: Field (Surface density per mass bin [g/cm²]) - SigmaFloor: Field (Floor value of surface density [g/cm²]) - St: Field (Stokes number) - v: Group (Velocities) - driftmax: Field (Maximum drift velocity [cm/s]) - frag: Field (Fragmentation velocity [cm/s]) - rad: Field (Radial velocity [cm/s]) - rel: Group (Relative velocities) - azi: Field (Relative azimuthal velocity [cm/s]) - brown: Field (Relative Brownian motion velocity [cm/s]) - rad: Field (Relative radial velocity [cm/s]) - tot: Field (Total relative velocity [cm/s]) - turb: Field (Relative turbulent velocity [cm/s]) - vert: Field (Relative vertical settling velocity [cm/s]) - gas: Group (Gas quantities) - alpha: Field (Turbulent alpha parameters) - boundary: Group (Boundary conditions) - inner: Constant gradient - outer: Value - cs: Field (Sound speed [cm/s]) - eta: Field (Pressure gradient parameter) - Fi: Field (Gas flux interfaces [g/cm/s]) - gamma: Field (Adiabatic index) - Hp: Field (Pressure scale height [cm]) - mfp: Field (Midplane mean free path [cm]) - mu: Field (Mean molecular weight [g]) - n: Field (Miplane number density [1/cm³]) - nu: Field (Kinematic viscosity [cm²/s]) - P: Field (Midplane pressure [g/cm/s²]) - rho: Field (Miplane mass density [g/cm³]) - S: Group (Source terms) - ext: Field (External sources [g/cm²/s]) - hyd: Field (Hydrodynamic sources [g/cm²/s]) - tot: Field (Total sources [g/cm²/s]) - Sigma: Field (Surface density [g/cm²]) - SigmaFloor: Field (Floor value of surface density [g/cm²]) - T: Field (Temperature [K]) - v: Group (Velocities) - rad: Field (Radial velocity [cm/s]) - visc: Field (Viscous accretion velocity [cm/s]) - grid: Group (Grid quantities) - A: Field (Radial grid annulus area [cm²]), constant - m: Field (Mass grid [g]), constant - Nm: Field (# of mass bins), constant - Nr: Field (# of radial grid cells), constant - OmegaK: Field (Keplerian frequency [1/s]) - r: Field (Radial grid cell centers [cm]), constant - ri: Field (Radial grid cell interfaces [cm]), constant - star: Group (Stellar quantities) - L: Field (Luminosity [erg/s]) - M: Field (Mass [g]) - R: Field (Radius [cm]) - T: Field (Effective temperature [K]) - t: Field (Time [s]), Integration variable
Running a Simulation¶
The simulation is now ready to go. This might take a few minutes. The default simulation is running for 100,000 years.
sim.run()
DustPy v0.5.1 Documentation: https://stammler.github.io/dustpy/ PyPI: https://pypi.org/project/dustpy/ GitHub: https://github.com/stammler/dustpy/ Please read README.md on the GitHub repository for information about the Terms of Usage. Checking for mass conservation... - Sticking: max. rel. error: 2.81e-14 for particle collision m[114] = 1.93e+04 g with m[116] = 3.73e+04 g - Full fragmentation: max. rel. error: 3.33e-16 for particle collision m[4] = 3.73e-12 g with m[6] = 7.20e-12 g - Cratering: max. rel. error: 1.78e-15 for particle collision m[110] = 5.18e+03 g with m[118] = 7.20e+04 g Creating data directory 'data'. Writing file data/data0000.hdf5 Writing dump file data/frame.dmp Writing file data/data0001.hdf5 Writing dump file data/frame.dmp Writing file data/data0002.hdf5 Writing dump file data/frame.dmp Writing file data/data0003.hdf5 Writing dump file data/frame.dmp Writing file data/data0004.hdf5 Writing dump file data/frame.dmp Writing file data/data0005.hdf5 Writing dump file data/frame.dmp Writing file data/data0006.hdf5 Writing dump file data/frame.dmp Writing file data/data0007.hdf5 Writing dump file data/frame.dmp Writing file data/data0008.hdf5 Writing dump file data/frame.dmp Writing file data/data0009.hdf5 Writing dump file data/frame.dmp Writing file data/data0010.hdf5 Writing dump file data/frame.dmp Writing file data/data0011.hdf5 Writing dump file data/frame.dmp Writing file data/data0012.hdf5 Writing dump file data/frame.dmp Writing file data/data0013.hdf5 Writing dump file data/frame.dmp Writing file data/data0014.hdf5 Writing dump file data/frame.dmp Writing file data/data0015.hdf5 Writing dump file data/frame.dmp Writing file data/data0016.hdf5 Writing dump file data/frame.dmp Writing file data/data0017.hdf5 Writing dump file data/frame.dmp Writing file data/data0018.hdf5 Writing dump file data/frame.dmp Writing file data/data0019.hdf5 Writing dump file data/frame.dmp Writing file data/data0020.hdf5 Writing dump file data/frame.dmp Writing file data/data0021.hdf5 Writing dump file data/frame.dmp Execution time: 0:25:54
By default DustPy has written output files into the data/ directory.
Plotting¶
DustPy is coming with a simple plotting script that can be used to check the status of a simulation.
from dustpy import plot
The plotting script does either take the simulation object as argument or a data directory.
If the argument is a simulation object the script is only plotting the current state.
plot.panel(sim)
The blue and the green lines in the top left panel are analytical estimates for the fragmentation and drift barrier taken from Birnstiel et al. (2012).
If you pass the data directory as argument, you also have access to the time evolution.
Furthermore, some plots can be addressed by specifying the time it, radial ir, or mass im index.
plot.panel("data", it=10, ir=10, im=5)
It's furthermore possible to use an interactive plotting script.
plot.ipanel("data")
The top middle and bottom left panels show the dust density distribution along the gray dashed lines in the top left panel.
Code Units¶
The dust density distribution plotted here is given by
$\Sigma_\mathrm{d} \left( r \right) = \int\limits_0^\infty \sigma \left(r, m \right) \mathrm{d} \log m$.
In this way the distribution is independent of the mass grid.
The code units DustPy uses are $\Sigma_{\mathrm{d},\,i} \left( r \right) \equiv \Sigma_\mathrm{d} \left(r, m_i \right)$ with
$\Sigma_\mathrm{d} \left( r \right) = \sum\limits_i \Sigma_{\mathrm{d},\,i} \left( r \right)$,
meaning the numerical sum over the mass dimension is the dust surface density.
sim.dust.Sigma.sum(-1)
[1.09467116e+01 1.01999127e+01 9.50180414e+00 8.85574147e+00 8.25193309e+00 7.68811165e+00 7.16111068e+00 6.66886490e+00 6.20952668e+00 5.78043584e+00 5.37983291e+00 5.00620617e+00 4.65734633e+00 4.33177619e+00 4.02824798e+00 3.74494515e+00 3.48065183e+00 3.23431318e+00 3.00445835e+00 2.79011526e+00 2.59034797e+00 2.40400028e+00 2.23032528e+00 2.06843796e+00 1.91746960e+00 1.77682330e+00 1.64573458e+00 1.52354221e+00 1.40971224e+00 1.30366897e+00 1.20490326e+00 1.11286820e+00 1.02721065e+00 9.47430890e-01 8.73155655e-01 8.04088985e-01 7.39755386e-01 6.79974636e-01 6.24338045e-01 5.72630768e-01 5.24619663e-01 4.79986454e-01 4.38580691e-01 4.00148755e-01 3.64525312e-01 3.31532873e-01 3.00968028e-01 2.72726922e-01 2.46597677e-01 2.22502736e-01 2.00270404e-01 1.79812878e-01 1.61002856e-01 1.43743896e-01 1.27930015e-01 1.13480732e-01 1.00296441e-01 8.83097336e-02 7.74274145e-02 6.75813725e-02 5.87121699e-02 5.07514056e-02 4.36518950e-02 3.73744694e-02 3.18521444e-02 2.70040286e-02 2.27480354e-02 1.90186800e-02 1.57674283e-02 1.29545972e-02 1.05430583e-02 8.49458255e-03 6.76826065e-03 5.31433734e-03 4.06259933e-03 2.91403465e-03 1.79320922e-03 8.10323703e-04 2.53491048e-04 5.35555727e-05 7.63359898e-06 8.05766160e-07 6.46162501e-08 3.98729596e-09 1.89679589e-10 6.88597892e-12 1.86172607e-13 3.59701757e-15 4.68098893e-17 3.80976200e-19 1.78887787e-21 4.50726625e-24 5.91422204e-27 9.37143393e-29 7.84597776e-29 6.83347699e-29 5.95170965e-29 5.18372240e-29 4.51483347e-29 3.93225557e-29]
To convert this into $\sigma$ we have to divide the code density distribution by $B = \frac{\Delta m}{m}$, where $\Delta m$ is the width of the mass bin centered on $m$.
Since the mass grid is strictly logarithmic the following relation holds
$m_{i+1} = A \cdot m_i$.
The grid constant $A$ can be easily calculated.
import numpy as np
A = np.mean(sim.grid.m[1:]/sim.grid.m[:-1])
A
1.3894954943731377
We further assume that the mass bin center is exactly in the middle between the mass bin interfaces
$m_i = \frac{1}{2} \left( m_{i-\frac{1}{2}} + m_{i+\frac{1}{2}} \right)$.
Solving for the interfaces yields
$\begin{split} m_{i-\frac{1}{2}} &= \frac{2}{A+1} m_i \\ m_{i+\frac{1}{2}} &= \frac{2}{A+1} A\cdot m_i. \end{split}$
We therefore have
$\Delta m_i = m_{i+\frac{1}{2}} - m_{i-\frac{1}{2}} = 2\frac{A-1}{A+1}m_i$
and
$B = \frac{\Delta m_i}{m_i} = 2\frac{A-1}{A+1}$.
B = 2 * (A-1) / (A+1)
sim.dust.Sigma / B
[[1.74750727e-01 1.42894613e-01 1.62086549e-01 ... 1.52526930e-23 2.11935481e-23 2.94483396e-23] [1.67458492e-01 1.37111787e-01 1.54857088e-01 ... 1.32845402e-23 1.84588088e-23 2.56484316e-23] [1.60965637e-01 1.31983478e-01 1.48366855e-01 ... 1.15703508e-23 1.60769503e-23 2.23388501e-23] ... [4.45718657e-46 6.19324066e-46 8.60547999e-46 ... 2.30858850e-29 3.20777332e-29 4.45718657e-29] [3.88204722e-46 5.39408712e-46 7.49505975e-46 ... 2.01069653e-29 2.79385376e-29 3.88204722e-29] [3.38112178e-46 4.69805348e-46 6.52792414e-46 ... 1.75124347e-29 2.43334490e-29 3.38112178e-29]]
Reading data files¶
If you want to read data files you can use the read/writer module provided by simframe, that is used to write the data.
from dustpy import hdf5writer
Make sure that the correct data directory is assigned to the writer.
hdf5writer
Writer (HDF5 file format using h5py)
------------------------------------
Data directory : data
File names : data/data0000.hdf5
Overwrite : False
Dumping : True
Options : {'com': 'lzf', 'comopts': None}
Verbosity : 1
You can now read a single data file with
data0003 = hdf5writer.read.output(3)
This function returns a namespace and the data can simply be accessed in the same way as for the Simulation object.
data0003.gas.Sigma
array([1.11372257e+003, 1.03888346e+003, 9.69003565e+002, 9.03557194e+002,
8.42299854e+002, 7.85001698e+002, 7.31441889e+002, 6.81404999e+002,
6.34679896e+002, 5.91060674e+002, 5.50348583e+002, 5.12353945e+002,
4.76897425e+002, 4.43810557e+002, 4.12935662e+002, 3.84125429e+002,
3.57242338e+002, 3.32158049e+002, 3.08752805e+002, 2.86914856e+002,
2.66539935e+002, 2.47530752e+002, 2.29796535e+002, 2.13252592e+002,
1.97819910e+002, 1.83424774e+002, 1.69998414e+002, 1.57476681e+002,
1.45799732e+002, 1.34911751e+002, 1.24760673e+002, 1.15297944e+002,
1.06478278e+002, 9.82594451e+001, 9.06020661e+001, 8.34694219e+001,
7.68272763e+001, 7.06437094e+001, 6.48889625e+001, 5.95352925e+001,
5.45568354e+001, 4.99294793e+001, 4.56307440e+001, 4.16396695e+001,
3.79367099e+001, 3.45036346e+001, 3.13234351e+001, 2.83802366e+001,
2.56592155e+001, 2.31465211e+001, 2.08292014e+001, 1.86951334e+001,
1.67329568e+001, 1.49320113e+001, 1.32822780e+001, 1.17743230e+001,
1.03992453e+001, 9.14862750e+000, 8.01448992e+000, 6.98924831e+000,
6.06567494e+000, 5.23686358e+000, 4.49619841e+000, 3.83732714e+000,
3.25413853e+000, 2.74074452e+000, 2.29146712e+000, 1.90083010e+000,
1.56355549e+000, 1.27456471e+000, 1.02898398e+000, 8.22153596e-001,
6.49640253e-001, 5.07251662e-001, 3.91052358e-001, 2.97379529e-001,
2.22857617e-001, 1.64410412e-001, 1.19269449e-001, 8.49776980e-002,
5.93877692e-002, 4.06542358e-002, 2.72200485e-002, 1.77974559e-002,
1.13442374e-002, 7.03640331e-003, 4.23874564e-003, 2.47473017e-003,
1.39717242e-003, 7.60960500e-004, 3.98793478e-004, 2.00545033e-004,
9.64877928e-005, 4.42747967e-005, 1.93105045e-005, 7.97641217e-006,
3.10820798e-006, 1.13786173e-006, 3.88322832e-007, 1.00000000e-100])
You can also read the full data directory with
data = hdf5writer.read.all()
The data has now an additional dimension for time.
data.gas.Sigma.shape
(22, 100)
Data files can be quite large and reading the entire data set can take some time. Is is also possible to only read a single field from the data files instead of the entire data set.
SigmaGas = hdf5writer.read.sequence("gas.Sigma")
SigmaGas.shape
(22, 100)
It is also possible to exclude certain fields from being written into the data files to save memory by setting their save attribute to False.
sim.dust.kernel.save = False
Reading Dump Files¶
The data files are only containing the pure data, but no information about the operations DustPy has to perform, like customized functions. Therefore, it's not possible to directly restart a simulation from data files.
simframe is saving a by default the most recent dump file, from which a simulation can be restarted.
Attention: Malware can be injected with dump files, which are pickled objects. Only read dump files that you have created yourself or from sources that you trust! Dump files have to be read with the same version of DustPy as they were written. Otherwise, it is not guaranteed to work.
from dustpy import readdump
sim_restart = readdump("data/frame.dmp")
This is now a simulation frame that should be identical to our previous object.
sim.gas.Sigma == sim_restart.gas.Sigma
[ True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True True]
We can now, for example, add more snapshots and restart the simulation. Here we just want to extend the run by one year.
from dustpy import constants as c
sim_restart.t.snapshots = np.concatenate((sim_restart.t.snapshots, [100001.*c.year]))
The current time is
sim_restart.t / c.year
100000.0
We can now restart the simulation for another year.
sim_restart.run()
DustPy v0.5.1 Documentation: https://stammler.github.io/dustpy/ PyPI: https://pypi.org/project/dustpy/ GitHub: https://github.com/stammler/dustpy/ Please read README.md on the GitHub repository for information about the Terms of Usage. Checking for mass conservation... - Sticking: max. rel. error: 2.81e-14 for particle collision m[114] = 1.93e+04 g with m[116] = 3.73e+04 g - Full fragmentation: max. rel. error: 3.33e-16 for particle collision m[4] = 3.73e-12 g with m[6] = 7.20e-12 g - Cratering: max. rel. error: 1.78e-15 for particle collision m[110] = 5.18e+03 g with m[118] = 7.20e+04 g Writing file data/data0022.hdf5 Writing dump file data/frame.dmp Execution time: 0:00:01
Another file was written and the current time is now
sim_restart.t / c.year
100001.0