Notebook

The Standard Model¶

In this chapter we'll go through all fields in the standard model, explain their meaning and the default functions that calculates them.

In [1]:

from dustpy import Simulation

In [2]:

sim = Simulation()

In [3]:

sim.initialize()

By default the frame object of DustPy consists of four groups for dust, gas, grid, and stellar parameters, a field for the time, which is the integration variable, an integrator object, and a writer object.

In [4]:

sim

Out[4]:

DustPy
------
    dust         : Group (Dust quantities)
    gas          : Group (Gas quantities)
    grid         : Group (Grid quantities)
    star         : Group (Stellar quantities)
  -----
    t            : Field (Time [s]), Integration variable
  -----
    Integrator   : Integrator (Default integrator)
    Writer       : Writer (HDF5 file format using h5py)

Dust¶

In [5]:

sim.dust

Out[5]:

Group (Dust quantities)
-----------------------
    backreaction : Group (Backreaction coefficients)
    boundary     : Group (Boundary conditions)
    coagulation  : Group (Coagulation quantities)
    delta        : Group (Mixing parameters)
    Fi           : Group (Fluxes)
    p            : Group (Probabilities)
    S            : Group (Sources)
    v            : Group (Velocities)
  -----
    a            : Field (Particle size [cm)
    D            : Field (Diffusivity [cm²/s])
    eps          : Field (Dust-to-gas ratio)
    fill         : Field (Filling factor)
    H            : Field (Scale heights [cm])
    kernel       : Field (Collision kernel [cm²/s])
    rho          : Field (Midplane mass density per mass bin [g/cm³])
    rhos         : Field (Solid state density [g/cm³])
    Sigma        : Field (Surface density per mass bin [g/cm²])
    SigmaFloor   : Field (Floor value of surface density [g/cm²])
    St           : Field (Stokes number)
  -----

`Simulation.dust.backreaction`¶

In [6]:

sim.dust.backreaction

Out[6]:

Group (Backreaction coefficients)
---------------------------------
    A            : Field (Pull factor)
    B            : Field (Push factor)
  -----

The backreaction describes the hydrodynamic influence the dust has on the gas. Numerically it consists of two fields A and B of shape (Simulation.grid.Nr,) that describe the pull respectively the push the dust excerts on the gas.

The details of this mechanism are described in Gárate et al. (2019).

Backreaction modifies the radial gas velocity as follows

$v_\mathrm{g} = Av_\mathrm{visc} + 2B\eta v_\mathrm{K}$.

In the standard model we have A=1 and B=0 everywhere, i.e., backreactions is not active.

In [7]:

sim.dust.backreaction.A

Out[7]:

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1.]

In [8]:

sim.dust.backreaction.B

Out[8]:

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0.]

Depending on the type of backreaction that you want to model, you have to provide functions for A and B. Have a look at Appendix A of Gárate et al. (2019) for examples.

`Simulation.dust.boundary`¶

In [9]:

sim.dust.boundary

Out[9]:

Group (Boundary conditions)
---------------------------
    inner        : Boundary
    outer        : Boundary
  -----

By default there are no dust boundary conditions set for the inner and outer grid cells of the dust.

In [10]:

sim.dust.boundary.inner

Out[10]:

No boundary condition set

In [11]:

sim.dust.boundary.outer

Out[11]:

No boundary condition set

The boundary is indirectly defined by the dust fluxes through the inner and outer grid interfaces, which are by default calculated by assuming constant dust velocities at the boundaries.

The boundary conditions can be changed via setcondition().

In [12]:

help(sim.dust.boundary.outer.setcondition)

Help on method setcondition in module dustpy.utils.boundary:

setcondition(condition, value=None) method of dustpy.utils.boundary.Boundary instance
    Function to set boundary condition.
    
    Parameters
    ----------
    condition : string
        Type of boundary conditon:
            - "const_grad" : constant gradient
            - "const_pow" : constant power law
            - "const_val" : constant value
            - "val" : custom value
            - "grad" : custom gradient
            - "pow" : custom power law with set exponent
            - None : Don't impose boundary condition (default)
    value : float or array, optional, default : None
        Value if needed for boundary condition

`Simulation.dust.coagulation`¶

The fields in this group define the behavior of dust growth and are discussed in a separate chapter.

`Simulation.dust.delta`¶

In [13]:

sim.dust.delta

Out[13]:

Group (Mixing parameters)
-------------------------
    rad          : Field (Radial mixing parameter)
    turb         : Field (Turbulent mixing parameter)
    vert         : Field (Vertical mixing parameter)
  -----

The $\delta$ parameters control the mixing of dust particles along vertical and radial directions and turbulent mixing. You can see them similar to the turbulent $\alpha$ parameter. And by default they will have the same value as $\alpha$ as given by Simulation.ini.gas.alpha.

`Simulation.dust.delta.rad`¶

In [14]:

sim.dust.delta.rad

Out[14]:

[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001]

$\delta_\mathrm{rad}$ will be used to calculate the radial, turbulent RMS velocity of the dust.

`Simulation.dust.delta.turb`¶

In [15]:

sim.dust.delta.turb

Out[15]:

[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001]

$\delta_\mathrm{turb}$ will be used to calculate the turbulent collision velocities of the dust particles.

`Simulation.dust.delta.vert`¶

In [16]:

sim.dust.delta.vert

Out[16]:

[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001]

$\delta_\mathrm{vert}$ will be used to calculate the vertical dust scale heights.

`Simulation.dust.Fi`¶

In [17]:

sim.dust.Fi

Out[17]:

Group (Fluxes)
--------------
    adv          : Field (Advective flux [g/cm/s)
    diff         : Field (Diffusive flux [g/cm/s)
    tot          : Field (Total flux [g/cm/s)
  -----

This is a group of fluxes through the radial grid interfaces for every particle mass of shape (Simulation.grid.Nr+1,, Simulation.grid.Nm).

`Simulation.dust.Fi.adv`¶

This is the advective flux calculated by $F_\mathrm{adv} = v_\mathrm{d}\Sigma_\mathrm{d}$.

`Simulation.dust.Fi.diff`¶

This is the diffusive flux calculated by $F_\mathrm{diff} = -D\Sigma_\mathrm{g}\nabla\frac{\Sigma_\mathrm{d}}{\Sigma_\mathrm{gas}}$ for every particle species separately. The diffusive fluxes at the grid boundaries are set to zero to avoid instabilities.

In [18]:

sim.dust.Fi.diff

Out[18]:

[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  0.00000000e+00
   0.00000000e+00  0.00000000e+00]
 [-1.50642471e-14  1.50642471e-14 -1.50642471e-14 ...  1.91213354e-21
   2.57703517e-21  3.42138925e-21]
 [-0.00000000e+00 -0.00000000e+00 -0.00000000e+00 ...  1.67822354e-21
   2.27209437e-21  3.03640156e-21]
 ...
 [-3.57009189e-47 -3.98385408e-47 -4.44555995e-47 ... -1.33086273e-41
  -1.48508841e-41 -1.65718639e-41]
 [-3.89114355e-48 -4.34207055e-48 -4.84525217e-48 ... -1.45046997e-42
  -1.61855623e-42 -1.80612098e-42]
 [ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  0.00000000e+00
   0.00000000e+00  0.00000000e+00]]

`Simulation.dust.Fi.tot`¶

This is the total flux through the radial grid interfaces $F_\mathrm{tot} = F_\mathrm{adv} + F_\mathrm{diff}$.

`Simulation.dust.p`¶

In [19]:

sim.dust.p

Out[19]:

Group (Probabilities)
---------------------
    frag         : Field (Fragmentation probability)
    stick        : Field (Sticking probability)
  -----

These are the fragmentation and sticking probability of certain particle collisions.
Their shape is (Simulation.grid.Nr, Simulation.grid.Nm, Simulation.grid.Nm).

The fragmentation probability of particle i=80, with particle j=3 at radial grid cell ir=30 is given by

In [20]:

ir = 30
i = 80
j = 3
sim.dust.p.frag[ir, j, i]

Out[20]:

1.0

`Simulation.dust.p.frag`¶

The fragmentation probability has a smooth transition between no fragmentation and fragmentation at the fragmentation velocity.

$p_\mathrm{f} = \begin{cases} 1, & \text{if } v_\mathrm{rel}>v_\mathrm{frag}\\ 0, & \text{if } v_\mathrm{rel}<v_\mathrm{frag} - \Delta v\\ \frac{v_\mathrm{rel} - v_\mathrm{frag} + \Delta v}{\Delta v}, & \text{otherwise} \end{cases}$

with $\Delta v = 0.2\,v_\mathrm{frag}$.

`Simulation.dust.p.stick`¶

This is the sticking probability given by $p_\mathrm{s} = 1 - p_\mathrm{f}$.

Bouncing can be easily implemented if $0 \leq p_\mathrm{f} + p_\mathrm{s} < 1$.

`Simulation.dust.S`¶

In [21]:

sim.dust.S

Out[21]:

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])
  -----

These are the source terms of the dust of shape (Simulation.grid.Nr, Simulation.grid.Nm) used to integrate the time evolution of the dust.

`Simulation.dust.S.coag`¶

The source terms from dust coagulation. This is described in detail in a separate chapter.

`Simulation.dust.S.ext`¶

External source terms. These are by default set to zero, i.e. no external sources.

In [22]:

sim.dust.S.ext

Out[22]:

[[0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 ...
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]
 [0. 0. 0. ... 0. 0. 0.]]

`Simulation.dust.S.hyd`¶

These are the hydrodynamic source terms. At grid cell $i$ they are calculated by

$S_{\mathrm{hyd},\,i} = 2\frac{\left( r_{i-\frac{1}{2}}F_{i-\frac{1}{2}}-r_{i+\frac{1}{2}}F_{i+\frac{1}{2}} \right)}{r_{i+\frac{1}{2}}^2 - r_{i-\frac{1}{2}}^2}$.

`Simulation.dust.S.tot`¶

These are the total source terms given by $S_\mathrm{tot} = S_\mathrm{coag} + S_\mathrm{ext} + S_\mathrm{hyd}$.

These are the actual sources the dust integrator uses.

`Simulation.dust.v`¶

In [23]:

sim.dust.v

Out[23]:

Group (Velocities)
------------------
    rel          : Group (Relative velocities)
  -----
    driftmax     : Field (Maximum drift velocity [cm/s])
    frag         : Field (Fragmentation velocity [cm/s])
    rad          : Field (Radial velocity [cm/s])
  -----

These are some dust related velocities the simulation needs for execution.

`Simulation.dust.v.rel`¶

In [24]:

sim.dust.v.rel

Out[24]:

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])
  -----

These are the different sources of relative particle collision velocities used by the coagulation algorithm.
They are used for two different reasons:

To calculated the collision rates
To determine the outcome of a collision

`Simulation.dust.v.rel.azi`¶

Relative collision velocity caused by a differential azimuthal drift of particles with different Stokes numbers calculated by

$v_{\mathrm{rel},\,\mathrm{azi}} = \left| v_{\mathrm{drift},\,\mathrm{max}} \cdot \left( \frac{1}{1+\mathrm{St}_i^2} - \frac{1}{1+\mathrm{St}_j^2} \right) \right|$.

`Simulation.dust.v.brown`¶

Relative collision velocity of particles caused by Brownian motion calculated with

$v_{\mathrm{rel},\,\mathrm{brown}} = \frac{8k_\mathrm{B}T\left(m_i + m_j \right)}{\pi m_i m_j}$.

Since this expression is diverging for very small particle masses, the relative velocity is capped to a maximum value of the sound speed $c_\mathrm{s}$. For very small particles this can still be larger than the fragmentation velocity and can cause unwanted fragmentation in the simple coagulation algorithm implemented in the default model.

`Simulation.dust.v.rel.rad`¶

Relative collision velocity caused by differential radial drift.

$v_{\mathrm{rel},\,\mathrm{rad}} = \left| v_{\mathrm{rad},\,i} - v_{\mathrm{rad},\,j} \right|$.

`Simulation.dust.v.rel.tot`¶

Total relative velocities calculated by using the root mean square of all individual velocity sources.

$v_{\mathrm{rel},\,\mathrm{tot}} = \sqrt{v_{\mathrm{rel},\,\mathrm{azi}}^2 + v_{\mathrm{rel},\,\mathrm{brown}}^2 + v_{\mathrm{rel},\,\mathrm{rad}}^2 + v_{\mathrm{rel},\,\mathrm{turb}}^2 + v_{\mathrm{rel},\,\mathrm{vert}}^2}$.

`Simulation.dust.v.rel.turb`¶

Relative velocities caused by turbulent motion of the particles. This follows the prescription of Ormel & Cuzzi (2007).
It uses Simulation.dust.delta.turb instead of Simulation.gas.alpha to calculate the velocities.

`Simulation.dust.v.rel.vert`¶

Relative collision velocities caused by differential vertical settling of particles.

$v_{\mathrm{rel},\,\mathrm{vert}} = \left| h_i \min \left( \mathrm{St}_i,\,\frac{1}{2}\right) - h_j \min \left( \mathrm{St}_j,\,\frac{1}{2}\right) \right| \cdot \Omega_\mathrm{K}$.

This prescription is taken from Birnstiel et al. (2010) and follows Dullemond & Dominik (2004).

`Simulation.dust.v.driftmax`¶

This is the maximum drift velocity a particle of $\mathrm{St} = 1$ can have.

$v_{\mathrm{drift},\,\mathrm{max}} = \frac{1}{2} B v_\mathrm{visc} - A \eta v_\mathrm{K}$.

See Simulation.dust.backreaction for details.

`Simulation.dust.v.frag`¶

Fragmentation velocities of shape (Simulation.grid.Nr,). By default this is set by the value of Simulation.ini.dust.vfrag.

In [25]:

sim.dust.v.frag

Out[25]:

[100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.
 100. 100.]

`Simulation.dust.v.rad`¶

Radial velocities of the dust particles.

$v_\mathrm{rad} = \left( v_\mathrm{g} + 2 v_{\mathrm{drift}\,\mathrm{max}} \right) \frac{\mathrm{St}}{\mathrm{St}^2+1}$

`Simulation.dust.a`¶

Particle radii calculated with

$a = \sqrt[3]{\frac{3m}{4\pi\rho_\mathrm{s}}}$

`Simulation.dust.D`¶

Dust diffusivity for every particle species of shape (Simulation.grid.Nr, Simulation.grid.Nm) calculated with

$D = \frac{\delta_\mathrm{rad}\,c_\mathrm{s}^2}{\Omega_\mathrm{K}\left( 1 + \mathrm{St}^2 \right)}$.

`Simulation.dust.eps`¶

This is the vertically integrated dust-to-gas ratio of shape (Simulation.grid.Nr,). In the literature this is also often refered to as metallicity $z$. It is calculated via

$\epsilon = \frac{\sum\limits_i \Sigma_{\mathrm{d},\,i}}{\Sigma_\mathrm{g}}$

`Simulation.dust.fill`¶

This describes the filling factor of the dust aggregates. By default this is 1.

In [26]:

sim.dust.fill

Out[26]:

[[1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 ...
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]
 [1. 1. 1. ... 1. 1. 1.]]

`Simulation.dust.H`¶

These are the dust scale heights of shape (Simulation.grid.Nr, Simulation.grid.Nm) calculated with the prescription of Dubrulle et al. (1995)

$H_i = H_\mathrm{P} \cdot \sqrt{\frac{\delta_\mathrm{vert}+\mathrm{St}_i}{\delta_\mathrm{vert}}}$

It is limited to a maximum of the pressure scale height $H_\mathrm{P}$. It uses Simulation.dust.delta.vert as vertical mixing parameter instead of Simulation.gas.alpha.

`Simulation.dust.kernel`¶

These are the vertically integrated collision kernels of shape (Simulation.grid.Nr, Simulation.grid.Nm, Simulation.grid.Nm) calculated with

$K = \frac{\pi\left( a_i + a_j \right)^2}{\sqrt{2\pi\left( H_i^2 + H_j^2 \right)}} v_\mathrm{rel}$

Please have a look at Birnstiel et al. (2010) for details. The kernel is the geometrical cross section multiplied with the factor described in (A.14). Multiplying the kernel with the densities returns the collision rates.

`Simulation.dust.rho`¶

Midplane mass densities of the dust.

$\rho = \frac{\Sigma}{\sqrt{2\pi}H}$

`Simulation.dust.rhos`¶

Solid state density of the particle material. This is initially set to Simulation.ini.dust.rhoMonomer.

In [27]:

sim.dust.rhos

Out[27]:

[[1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 ...
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]]

To calculate the mass density of a dust aggregate do

In [28]:

sim.dust.fill * sim.dust.rhos

Out[28]:

[[1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 ...
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]
 [1.67 1.67 1.67 ... 1.67 1.67 1.67]]

This is the density that is used to calculate the particle size.

`Simulation.dust.Sigma`¶

This is the dust surface density of every particle species of shape (Simulation.grid.Nr, Simulation.grid.Nm). This is the quantity that is integrated to solve for dust evolution.

The surface densities are integrated of the mass bin, so a numerical summation over the mass dimension returns the total dust surface density.

$\Sigma_\mathrm{d} = \sum\limits_i \Sigma_{\mathrm{d},\,i}$

In [29]:

sim.dust.Sigma.sum(-1)

Out[29]:

[1.11824742e+01 1.04232155e+01 9.71465059e+00 9.05339692e+00
 8.43629793e+00 7.86040780e+00 7.32297742e+00 6.82144128e+00
 6.35340519e+00 5.91663488e+00 5.50904531e+00 5.12869076e+00
 4.77375546e+00 4.44254502e+00 4.13347826e+00 3.84507969e+00
 3.57597248e+00 3.32487186e+00 3.09057900e+00 2.87197527e+00
 2.66801690e+00 2.47773001e+00 2.30020595e+00 2.13459693e+00
 1.98011200e+00 1.83601326e+00 1.70161231e+00 1.57626698e+00
 1.45937826e+00 1.35038739e+00 1.24877321e+00 1.15404966e+00
 1.06576343e+00 9.83491798e-01 9.06840556e-01 8.35442144e-01
 7.68953850e-01 7.07056150e-01 6.49451154e-01 5.95861146e-01
 5.46027226e-01 4.99708026e-01 4.56678518e-01 4.16728886e-01
 3.79663471e-01 3.45299782e-01 3.13467557e-01 2.84007884e-01
 2.56772375e-01 2.31622377e-01 2.08428234e-01 1.87068586e-01
 1.67429708e-01 1.49404882e-01 1.32893806e-01 1.17802035e-01
 1.04040457e-01 9.15247991e-02 8.01751688e-02 6.99156307e-02
 6.06738170e-02 5.23805770e-02 4.49696665e-02 3.83774784e-02
 3.25428186e-02 2.74067268e-02 2.29123459e-02 1.90048396e-02
 1.56313572e-02 1.27410459e-02 1.02851052e-02 8.21687989e-03
 6.49198428e-03 5.06844837e-03 3.90687641e-03 2.97060545e-03
 2.22585132e-03 8.21564634e-26 7.15552883e-26 6.23220508e-26
 5.42802371e-26 4.72761102e-26 4.11757706e-26 3.58625970e-26
 3.12350162e-26 2.72045619e-26 2.36941829e-26 2.06367706e-26
 1.79738758e-26 1.56545914e-26 1.36345791e-26 1.18752220e-26
 1.03428860e-26 9.00827708e-27 7.84588135e-27 6.83347698e-27
 5.95170965e-27 5.18372240e-27 4.51483347e-27 3.93225557e-27]

If you want to plot the dust density distribution you might want to convert the field into a quantity that does not depend on the mass bin width. The plotting script integrated in DustPy if plotting a grid independent density distribution. See the first chapter for details on this.

`Simulation.dust.SigmaFloor`¶

This is the floor value for the dust surface densities. Mass bins that are below their respective floor value will not contribute to coagulation. By default the floor value is the density that corresponds to one physical particle of that mass distributed over and annulus at that radial grid location.

$\Sigma_{\mathrm{d},\,\mathrm{floor}} = \frac{m}{A_\mathrm{annulus}}$

Densities below the floor value therefore correspond to fewer that one physical particle.

`Simulation.dust.St`¶

Stokes numbers of the particles. By default the Epstein and the Stokes I drag regime are considered.

$\mathrm{St} = \begin{cases} \frac{\pi}{2} \frac{a\rho}{\Sigma_\mathrm{g}}, & \text{if } a < \frac{9}{4} \lambda_\mathrm{mfp}\\ \frac{2\pi}{9} \frac{a^2 \rho}{\lambda_\mathrm{mfp} \Sigma_\mathrm{g}}, & \text{else} \end{cases}$

Gas¶

In [30]:

sim.gas

Out[30]:

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])
  -----

`Simulation.gas.boundary`¶

In [31]:

sim.gas.boundary

Out[31]:

Group (Boundary conditions)
---------------------------
    inner        : Boundary
    outer        : Boundary
  -----

These are the boundary conditions of the gas. They are set to constant gradient by default.

In [32]:

sim.gas.boundary.inner

Out[32]:

Constant gradient

In [33]:

sim.gas.boundary.outer

Out[33]:

Value (1e-100)

The boundary conditions can be modified with setcondition.

In [34]:

help(sim.gas.boundary.inner.setcondition)

Help on method setcondition in module dustpy.utils.boundary:

setcondition(condition, value=None) method of dustpy.utils.boundary.Boundary instance
    Function to set boundary condition.
    
    Parameters
    ----------
    condition : string
        Type of boundary conditon:
            - "const_grad" : constant gradient
            - "const_pow" : constant power law
            - "const_val" : constant value
            - "val" : custom value
            - "grad" : custom gradient
            - "pow" : custom power law with set exponent
            - None : Don't impose boundary condition (default)
    value : float or array, optional, default : None
        Value if needed for boundary condition

If the gas surface density follows a power law $\propto R^{-1}$ the constant gradient boundary condition should work fine. Other values can lead to deviations at the inner boundary. See the chapter about gas evolution tests for details.

`Simulation.gas.S`¶

In [35]:

sim.gas.S

Out[35]:

Group (Source terms)
--------------------
    ext          : Field (External sources [g/cm²/s])
    hyd          : Field (Hydrodynamic sources [g/cm²/s])
    tot          : Field (Total sources [g/cm²/s])
  -----

These are the source terms of the gas.

`Simulation.gas.S.ext`¶

These are the external source terms for gas evolution, e.g. infall. By default these are set to zero.

In [36]:

sim.gas.S.ext

Out[36]:

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0.]

`Simulation.gas.S.hyd`¶

These are the hydrodynamic source terms of the gas evolution.
Attention: Since the gas evolution is calculated implicitly, the hydrodynamic source terms are calculated in retrospect after the new gas surface density was found. Changing Simulation.gas.S.hyd does not influence the simulation. It is only given for data analysis.

`Simulation.gas.S.tot`¶

These are the total source terms of gas evolution.

$S_\mathrm{tot} = S_\mathrm{ext} + S_\mathrm{hyd}$

Attention: Since the gas evolution is calculated implicitly, the total source terms are calculated in retrospect after the new gas surface density was found. Changing Simulation.gas.S.tot does not influence the simulation. It is only given for data analysis.

`Simulation.gas.v`¶

In [37]:

sim.gas.v

Out[37]:

Group (Velocities)
------------------
    rad          : Field (Radial velocity [cm/s])
    visc         : Field (Viscous accretion velocity [cm/s])
  -----

These are velocities that are relevant for the gas evolution.
Attention: Since the gas evolution is calculated implicitly, the velocities are calculated in retrospect after the new gas surface density was found. Changing anything in Simulation.gas.v does not influence the simulation. It is only given for data analysis.

`Simulation.gas.v.rad`¶

This is the radial gas velocity. It is given by

$v_\mathrm{g} = Av_\mathrm{visc} + 2B\eta v_\mathrm{K}$.

See Simulation.dust.backreaction for details. If backreaction is turned off, i.e., $A=1$ and $B=0$, the radial velocity is identical to the viscous velocity.
Attention: Since the gas evolution is calculated implicitly, the velocities are calculated in retrospect after the new gas surface density was found. Changing anything in Simulation.gas.v.rad does not influence the simulation. It is only given for data analysis.

`Simulation.gas.v.visc`¶

This is the radial viscous gas velocity

$v_\mathrm{visc} = -\frac{3}{\Sigma_\mathrm{g}\sqrt{R}} \frac{\partial}{\partial R} \left( \Sigma_\mathrm{g} \nu \sqrt{R} \right)$

Attention: Since the gas evolution is calculated implicitly, the velocities are calculated in retrospect after the new gas surface density was found. Changing anything in Simulation.gas.v.visc does not influence the simulation. It is only given for data analysis.

`Simulation.gas.alpha`¶

This is the turbulent viscosity parameter according Shakura & Sunyaev (1973). It is initially set to the value in Simulation.ini.gas.alpha.

In [38]:

sim.gas.alpha

Out[38]:

[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
 0.001 0.001 0.001 0.001]

`Simulation.gas.cs`¶

This is the adiabatic sound speed in the midplane of the disk

$c_\mathrm{s} = \sqrt{\frac{\gamma k_\mathrm{B} T}{\mu}}$.

For isothermal simulations set $\gamma=1$.

`Simulation.gas.eta`¶

This is the midplane pressure gradient parameter $\eta$ given by

$\eta = -\frac{1}{2} \left( \frac{H_\mathrm{P}}{r} \right)^2 \frac{\partial \log P}{\partial \log r}$

It describes the degree of "sub-Keplerity" of the disk

$v_\phi^2 = \left( 1-2\eta \right) v_\mathrm{K}^2$.

`Simulation.gas.Fi`¶

These are the mass fluxes of gas through the grid cell interfaces.
Attention: Since the gas evolution is calculated implicitly, the fluxes are calculated in retrospect after the new gas surface density was found. Changing anything in Simulation.gas.Fi does not influence the simulation. It is only given for data analysis.

`Simulation.gas.gamma`¶

This is the ratio of specific heats. It is initially set to the value given by Simulation.ini.gas.gamma.

In [39]:

sim.gas.gamma

Out[39]:

[1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4
 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4
 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4
 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4
 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4
 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4]

For isothermal simulations set this to $1$.

`Simulation.gas.Hp`¶

Pressure scale height of the gas given by the ratio of the isothermal sound speed to the Keplerian frequency.

$H_\mathrm{P} = \frac{c_{\mathrm{s},\,\mathrm{iso}}}{\Omega_\mathrm{K}}$.

`Simulation.gas.mfp`¶

Mean free path of the gas in the midplane of the disk

$\lambda_\mathrm{mfp} = \frac{1}{\sqrt{2}\,n\,\sigma_\mathrm{H_2}}$

`Simulation.gas.mu`¶

Mean molecular weight of the gas. This is initially set by the value given in Simulation.ini.gas.mu and is equal to $2.3\,m_\mathrm{P}$.

In [40]:

sim.gas.mu

Out[40]:

[3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24
 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24]

`Simulation.gas.n`¶

Midplane number density of the gas given by

$n = \frac{\rho}{\mu}$

`Simulation.gas.nu`¶

Kinematic viscosity of the gas given by

$\nu = \alpha c_\mathrm{s} H_\mathrm{P}$.

`Simulation.gas.P`¶

Midplane gas pressure given by

$P = \frac{\rho\,c_\mathrm{s}^2}{\gamma}$.

`Simulation.gas.rho`¶

Midplane gas mass density given by

$\rho = \frac{\Sigma_\mathrm{g}}{\sqrt{2\pi}H_\mathrm{P}}$.

`Simulation.gas.Sigma`¶

Gas surface density. This is the quantity that is integrated with an implicit Euler first-order scheme.

`Simulation.gas.SigmaFloor`¶

This is the floor value of the gas surface density. By default it is $10^{-100}$.

In [41]:

sim.gas.SigmaFloor

Out[41]:

[1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100 1.e-100
 1.e-100]

If the gas surface density is at any point below it's floor value it will be automatically set to the floor value at the end of a time step.

`Simulation.gas.T`¶

This is the midplane gas temperature. It is calculated by assuming a passively irradiated disk with a constant irradiation angle of $0.05$.

$T\left( r \right) = \sqrt[4]{\frac{0.05\,L_*}{4\,\pi\,r^2\,\sigma_\mathrm{SB}}}$

Grid¶

In [42]:

sim.grid

Out[42]:

Group (Grid quantities)
-----------------------
    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
  -----

These are all quantities that define the radial and the mass grid. Once they are defined they are constant and should not be changed. Additionally, the Keplerian frequency is located here.

`Simulation.grid.m`¶

The mass grid. It has to be strictly logarithmic. Please only use Simulation.ini.dust.mmin, Simulation.ini.dust.mmax, and Simulation.ini.dust.Nmbpd to set the mass grid and do not set it manually!

In [43]:

sim.grid.m

Out[43]:

[1.00000000e-12 1.38949549e-12 1.93069773e-12 2.68269580e-12
 3.72759372e-12 5.17947468e-12 7.19685673e-12 1.00000000e-11
 1.38949549e-11 1.93069773e-11 2.68269580e-11 3.72759372e-11
 5.17947468e-11 7.19685673e-11 1.00000000e-10 1.38949549e-10
 1.93069773e-10 2.68269580e-10 3.72759372e-10 5.17947468e-10
 7.19685673e-10 1.00000000e-09 1.38949549e-09 1.93069773e-09
 2.68269580e-09 3.72759372e-09 5.17947468e-09 7.19685673e-09
 1.00000000e-08 1.38949549e-08 1.93069773e-08 2.68269580e-08
 3.72759372e-08 5.17947468e-08 7.19685673e-08 1.00000000e-07
 1.38949549e-07 1.93069773e-07 2.68269580e-07 3.72759372e-07
 5.17947468e-07 7.19685673e-07 1.00000000e-06 1.38949549e-06
 1.93069773e-06 2.68269580e-06 3.72759372e-06 5.17947468e-06
 7.19685673e-06 1.00000000e-05 1.38949549e-05 1.93069773e-05
 2.68269580e-05 3.72759372e-05 5.17947468e-05 7.19685673e-05
 1.00000000e-04 1.38949549e-04 1.93069773e-04 2.68269580e-04
 3.72759372e-04 5.17947468e-04 7.19685673e-04 1.00000000e-03
 1.38949549e-03 1.93069773e-03 2.68269580e-03 3.72759372e-03
 5.17947468e-03 7.19685673e-03 1.00000000e-02 1.38949549e-02
 1.93069773e-02 2.68269580e-02 3.72759372e-02 5.17947468e-02
 7.19685673e-02 1.00000000e-01 1.38949549e-01 1.93069773e-01
 2.68269580e-01 3.72759372e-01 5.17947468e-01 7.19685673e-01
 1.00000000e+00 1.38949549e+00 1.93069773e+00 2.68269580e+00
 3.72759372e+00 5.17947468e+00 7.19685673e+00 1.00000000e+01
 1.38949549e+01 1.93069773e+01 2.68269580e+01 3.72759372e+01
 5.17947468e+01 7.19685673e+01 1.00000000e+02 1.38949549e+02
 1.93069773e+02 2.68269580e+02 3.72759372e+02 5.17947468e+02
 7.19685673e+02 1.00000000e+03 1.38949549e+03 1.93069773e+03
 2.68269580e+03 3.72759372e+03 5.17947468e+03 7.19685673e+03
 1.00000000e+04 1.38949549e+04 1.93069773e+04 2.68269580e+04
 3.72759372e+04 5.17947468e+04 7.19685673e+04 1.00000000e+05]

`Simulation.grid.Nm`¶

Number of mass bins.

In [44]:

sim.grid.Nm

Out[44]:

`Simulation.grid.Nr`¶

Number of radial grid cells.

In [45]:

sim.grid.Nr

Out[45]:

`Simulation.grid.OmegaK`¶

Keplerian frequency given by

$\Omega_\mathrm{K} = \sqrt{\frac{G\,M_*}{r^3}}$.

`Simulation.grid.r`¶

Radial grid cell centers. The radial grid cell centers are exactly in the middle between the grid cell interfaces.

In [46]:

sim.grid.r == 0.5 * (sim.grid.ri[:-1] + sim.grid.ri[1:])

Out[46]:

[ 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]

`Simulation.grid.ri`¶

Locations of the grid cell interfaces. Simulation.grid.ri[0] and Simulation.grid.ri[-1] are the inner and outer grid boundaries.

Star¶

In [47]:

sim.star

Out[47]:

Group (Stellar quantities)
--------------------------
    L            : Field (Luminosity [erg/s])
    M            : Field (Mass [g])
    R            : Field (Radius [cm])
    T            : Field (Effective temperature [K])
  -----

`Simulation.star.L`¶

Stellar luminosity given by

$L = 4\pi\,R_*^2\,\sigma_\mathrm{SB}\,T_*^4$

`Simulation.star.M`¶

Stellar mass. Initially set to the value given in Simulation.ini.star.M, which corresponds to one Solar mass.

In [48]:

sim.star.M

Out[48]:

1.988409870698051e+33

`Simulation.star.R`¶

Stellar radius. Initially set to the value given in Simulation.ini.star.R, which corresponds to 2 Solar radii.

In [49]:

sim.star.R

Out[49]:

139140000000.0

`Simulation.star.T`¶

Stellar effective surface temperature. Initially set to the value given in Simulation.ini.star.T.

In [50]:

sim.star.T

Out[50]:

5772.0

Time¶

Simulation.t is the current time and the integration variable. It starts at zero initially.

In [51]:

sim.t

Out[51]:

0.0

Snapshots are written between $10^3$ years and $10^5$ years with 10 snapshots per time decade.

In [52]:

sim.t.snapshots

Out[52]:

array([3.15576000e+10, 3.97286646e+10, 5.00154254e+10, 6.29656900e+10,
       7.92691073e+10, 9.97938935e+10, 1.25633068e+11, 1.58162662e+11,
       1.99114995e+11, 2.50670927e+11, 3.15576000e+11, 3.97286646e+11,
       5.00154254e+11, 6.29656900e+11, 7.92691073e+11, 9.97938935e+11,
       1.25633068e+12, 1.58162662e+12, 1.99114995e+12, 2.50670927e+12,
       3.15576000e+12])

Since the dust is integrated with an adaptive stepsize scheme, the timestep function will return the suggested stepsize of the dust integrator, which is initially set to one year.

In [53]:

sim.t.suggested

Out[53]:

31557600.0

In [54]:

sim.t.stepsize

Out[54]:

31557600.0

Integrator¶

The integrator that is used by default as two integration instructions. One for the gas and one for the dust.

In [55]:

sim.integrator

Out[55]:

Integrator (Default integrator)

In [56]:

sim.integrator.instructions

Out[56]:

[Instruction (Dust: explicit 5th-order adaptive Cash-Karp method),
 Instruction (Gas: implicit 1st-order Euler method)]

The gas is integrated with an implicit first-order Euler integration scheme. The Jacobian is calculated from the parameters given in Simulation.gas. Parameters like gas velocities, fluxes and source terms are calculated once the new values of the gas surface density have been found.

Except for the values at the boundaries, the advective source terms are given by

In [57]:

sim.gas.Sigma.jacobian() @ sim.gas.Sigma

Out[57]:

[ 3.54124798e-05 -1.93788620e-11 -1.80453135e-11 -1.68008509e-11
 -1.56395333e-11 -1.45558169e-11 -1.35445284e-11 -1.26008399e-11
 -1.17202464e-11 -1.08985441e-11 -1.01318100e-11 -9.41638378e-12
 -8.74884973e-12 -8.12602081e-12 -7.54492330e-12 -7.00278259e-12
 -6.49700996e-12 -6.02519016e-12 -5.58506990e-12 -5.17454708e-12
 -4.79166070e-12 -4.43458153e-12 -4.10160329e-12 -3.79113452e-12
 -3.50169095e-12 -3.23188835e-12 -2.98043588e-12 -2.74612994e-12
 -2.52784828e-12 -2.32454465e-12 -2.13524373e-12 -1.95903636e-12
 -1.79507516e-12 -1.64257036e-12 -1.50078596e-12 -1.36903605e-12
 -1.24668147e-12 -1.13312657e-12 -1.02781620e-12 -9.30232957e-13
 -8.39894458e-13 -7.56350873e-13 -6.79182551e-13 -6.07997779e-13
 -5.42430665e-13 -4.82139127e-13 -4.26802977e-13 -3.76122108e-13
 -3.29814763e-13 -2.87615883e-13 -2.49275545e-13 -2.14557473e-13
 -1.83237630e-13 -1.55102892e-13 -1.29949808e-13 -1.07583453e-13
 -8.78163775e-14 -7.04676704e-14 -5.53621362e-14 -4.23296062e-14
 -3.12043925e-14 -2.18248951e-14 -1.40333743e-14 -7.67589565e-15
 -2.60245162e-15  1.33274078e-15  4.27077277e-15  6.34742190e-15
  7.69239796e-15  8.42848892e-15  8.67064796e-15  8.52507094e-15
  8.08832069e-15  7.44655832e-15  6.67494233e-15  5.83725208e-15
  4.98578315e-15  4.16154774e-15  3.39479409e-15  2.70583646e-15
  2.10616296e-15  1.59976499e-15  1.18461279e-15  8.54187739e-16
  5.98977857e-16  4.07847906e-16  2.69210409e-16  1.71946744e-16
  1.06055166e-16  6.30311837e-17  3.60109638e-17  1.97267674e-17
  1.03327890e-17  5.15971756e-18  2.44844236e-18  1.10031050e-18
  4.66550552e-19  1.85917852e-19  6.93330477e-20  0.00000000e+00]

The total source terms, i.e., including external sources and excluding the boundaries are given by

In [58]:

sim.gas.Sigma.jacobian() @ sim.gas.Sigma + sim.gas.S.ext

Out[58]:

[ 3.54124798e-05 -1.93788620e-11 -1.80453135e-11 -1.68008509e-11
 -1.56395333e-11 -1.45558169e-11 -1.35445284e-11 -1.26008399e-11
 -1.17202464e-11 -1.08985441e-11 -1.01318100e-11 -9.41638378e-12
 -8.74884973e-12 -8.12602081e-12 -7.54492330e-12 -7.00278259e-12
 -6.49700996e-12 -6.02519016e-12 -5.58506990e-12 -5.17454708e-12
 -4.79166070e-12 -4.43458153e-12 -4.10160329e-12 -3.79113452e-12
 -3.50169095e-12 -3.23188835e-12 -2.98043588e-12 -2.74612994e-12
 -2.52784828e-12 -2.32454465e-12 -2.13524373e-12 -1.95903636e-12
 -1.79507516e-12 -1.64257036e-12 -1.50078596e-12 -1.36903605e-12
 -1.24668147e-12 -1.13312657e-12 -1.02781620e-12 -9.30232957e-13
 -8.39894458e-13 -7.56350873e-13 -6.79182551e-13 -6.07997779e-13
 -5.42430665e-13 -4.82139127e-13 -4.26802977e-13 -3.76122108e-13
 -3.29814763e-13 -2.87615883e-13 -2.49275545e-13 -2.14557473e-13
 -1.83237630e-13 -1.55102892e-13 -1.29949808e-13 -1.07583453e-13
 -8.78163775e-14 -7.04676704e-14 -5.53621362e-14 -4.23296062e-14
 -3.12043925e-14 -2.18248951e-14 -1.40333743e-14 -7.67589565e-15
 -2.60245162e-15  1.33274078e-15  4.27077277e-15  6.34742190e-15
  7.69239796e-15  8.42848892e-15  8.67064796e-15  8.52507094e-15
  8.08832069e-15  7.44655832e-15  6.67494233e-15  5.83725208e-15
  4.98578315e-15  4.16154774e-15  3.39479409e-15  2.70583646e-15
  2.10616296e-15  1.59976499e-15  1.18461279e-15  8.54187739e-16
  5.98977857e-16  4.07847906e-16  2.69210409e-16  1.71946744e-16
  1.06055166e-16  6.30311837e-17  3.60109638e-17  1.97267674e-17
  1.03327890e-17  5.15971756e-18  2.44844236e-18  1.10031050e-18
  4.66550552e-19  1.85917852e-19  6.93330477e-20  0.00000000e+00]

The dust is integrated with a 5th-order Cash-Karp adaptive stepsize scheme using Simulation.dust.S.tot as source terms.

At the end of a successful integration step, the floor values and boundaries are enforced.

Writer¶

DustPy uses by default the hdf5writer of simframe.

In [59]:

sim.writer

Out[59]:

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

It is by default writing dump files and prevents overwriting of already existing data files.

The Standard Model¶

Dust¶

Simulation.dust.backreaction¶

Simulation.dust.boundary¶

Simulation.dust.coagulation¶

Simulation.dust.delta¶

Simulation.dust.delta.rad¶

Simulation.dust.delta.turb¶

Simulation.dust.delta.vert¶

Simulation.dust.Fi¶

Simulation.dust.Fi.adv¶

Simulation.dust.Fi.diff¶

Simulation.dust.Fi.tot¶

Simulation.dust.p¶

Simulation.dust.p.frag¶

Simulation.dust.p.stick¶

Simulation.dust.S¶

Simulation.dust.S.coag¶

Simulation.dust.S.ext¶

Simulation.dust.S.hyd¶

Simulation.dust.S.tot¶

Simulation.dust.v¶

Simulation.dust.v.rel¶

Simulation.dust.v.rel.azi¶

Simulation.dust.v.brown¶

Simulation.dust.v.rel.rad¶

Simulation.dust.v.rel.tot¶

Simulation.dust.v.rel.turb¶

Simulation.dust.v.rel.vert¶

Simulation.dust.v.driftmax¶

Simulation.dust.v.frag¶

Simulation.dust.v.rad¶

Simulation.dust.a¶

Simulation.dust.D¶

Simulation.dust.eps¶

Simulation.dust.fill¶

Simulation.dust.H¶

Simulation.dust.kernel¶

Simulation.dust.rho¶

Simulation.dust.rhos¶

Simulation.dust.Sigma¶

Simulation.dust.SigmaFloor¶

Simulation.dust.St¶