4. The Standard Model
In this chapter we’ll go through all fields in the standard model, explain their meaning and the default functions that calculate them.
[1]:
from dustpy import Simulation
[2]:
sim = Simulation()
[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.
[4]:
sim
[4]:
DustPy
------
dust : Group (Dust quantities)
gas : Group (Gas quantities)
grid : Group (Grid quantities)
star : Group (Stellar quantities)
-----
t : IntVar (Time [s]), Integration variable
-----
Integrator : Integrator (Default integrator)
Writer : Writer (HDF5 file format using h5py)
Dust
[5]:
sim.dust
[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
[6]:
sim.dust.backreaction
[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.
[7]:
sim.dust.backreaction.A
[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.]
[8]:
sim.dust.backreaction.B
[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
[9]:
sim.dust.boundary
[9]:
Group (Boundary conditions)
---------------------------
inner : Boundary
outer : Boundary
-----
By default the inner dust boundary is set to constant gradient and the outer boundary to floor value.
[10]:
sim.dust.boundary.inner
[10]:
Constant gradient
[11]:
sim.dust.boundary.outer
[11]:
Value
The outer boundary therefore stores an array with the floor values of all particle masses.
The boundary conditions can be changed via setcondition().
[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 about dust coagulation.
Simulation.dust.delta
[13]:
sim.dust.delta
[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 look at them as similar to the turbulent \(\alpha\) parameter for the gas. And by default they will have the same value as \(\alpha\) as given by Simulation.ini.gas.alpha.
Simulation.dust.delta.rad
[14]:
sim.dust.delta.rad
[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 dust diffusion.
Simulation.dust.delta.turb
[15]:
sim.dust.delta.turb
[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
[16]:
sim.dust.delta.vert
[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
[17]:
sim.dust.Fi
[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).
Attention: When using implicit dust integration, the interface fluxes are calculated in retrospect after the new dust surface density was found. Changing any Simulation.dust.Fi does not influence the simulation. It is only given for data analysis.
Simulation.dust.Fi.adv
This is the advective flux calculated by \(F_\mathrm{adv} = v_\mathrm{d}\Sigma_\mathrm{d}\).
Attention: When using implicit dust integration, the advective interface fluxes are calculated in retrospect after the new dust surface density was found. Changing Simulation.dust.Fi.adv does not influence the simulation. It is only given for data analysis.
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.
Attention: When using implicit dust integration, the diffusive interface fluxes are calculated in retrospect after the new dust surface density was found. Changing Simulation.dust.Fi.diff does not influence the simulation. It is only given for data analysis.
Simulation.dust.Fi.tot
This is the total flux through the radial grid interfaces \(F_\mathrm{tot} = F_\mathrm{adv} + F_\mathrm{diff}\).
Attention: When using implicit dust integration, the total interface fluxes are calculated in retrospect after the new dust surface density was found. Changing Simulation.dust.Fi.tot does not influence the simulation. It is only given for data analysis.
Simulation.dust.p
[18]:
sim.dust.p
[18]:
Group (Probabilities)
---------------------
frag : Field (Fragmentation probability)
stick : Field (Sticking probability)
-----
(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
[19]:
ir = 30
i = 80
j = 3
sim.dust.p.frag[ir, j, i]
[19]:
0.9948352476770317
Simulation.dust.p.frag
DustPy assumes that the relative collision velocities follow the Maxwell-Boltzmann distribution
\(f\left(\Delta v\right) = \sqrt{\frac{54}{\pi}}\frac{\Delta v^{2}}{v_\mathrm{rms}^3} \exp \left[-\frac{3}{2} \frac{\Delta v^2}{v_\mathrm{rms}^2} \right]\).
The collision kernel for fragmentation is then given by
\(K_\mathrm{f} = \int\limits_{v_\mathrm{frag}}^\infty \sigma_\mathrm{geo} \Delta v f\left(\Delta v\right) \mathrm{d}\Delta v\)
by counting all collisions above the fragmentation velocity.
The fragmentation probability is the given by
\(p_\mathrm{f} = \frac{K_\mathrm{f}}{\sigma_\mathrm{geo} \Delta \bar{v}}\)
with the mean velocity of the Maxwell-Boltzmann distribution
\(\Delta \bar{v} = \sqrt{\frac{8\pi}{2}}v_\mathrm{rms}\).
The root-mean-square velocities \(v_\mathrm{rms}\) are stored in Simulation.dust.v.rel.
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
[20]:
sim.dust.S
[20]:
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 about dust coagulation.
Attention: When using implicit dust integration, the coagulation source terms are calculated in retrospect after the new dust surface density was found. Changing Simulation.dust.S.coag does not influence the simulation. It is only given for data analysis.
Simulation.dust.S.ext
External source terms. These are by default set to zero, i.e. no external sources.
[21]:
sim.dust.S.ext
[21]:
[[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}\).
Attention: When using implicit dust integration, the hydrodynamic source terms are calculated in retrospect after the new dust surface density was found. Changing Simulation.dust.S.hyd does not influence the simulation. It is only given for data analysis.
Simulation.dust.S.tot
These are the total source terms given by \(S_\mathrm{tot} = S_\mathrm{coag} + S_\mathrm{ext} + S_\mathrm{hyd}\).
Simulation.dust.v
[22]:
sim.dust.v
[22]:
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
[23]:
sim.dust.v.rel
[23]:
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])
-----
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}} = \sqrt{ \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
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.
[24]:
sim.dust.v.frag
[24]:
[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}}\mathrm{St} \right) \frac{1}{\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.
[25]:
sim.dust.fill
[25]:
[[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}}{\delta_\mathrm{vert}+\mathrm{St}_i}}\)
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_{ij} = \frac{1}{1+\delta_{ij}} \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.
The collisions rates of equal particle collisions are reduced by a factor of \(\frac{1}{2}\) by the Kronecker \(\delta_{ij}\). In a reservoir of \(N\) equal particles, the first particle can collide with \(N-1\) other particles, the second particle can collide with \(N-2\) other particles without counting the collision with the first particle twice. The third particle can collide with \(N-3\) other particles and so on. The total number of possible collisions is therefore given by
\(\sum\limits_{i=1}^{N-1} i = \frac{1}{2}\left(N\right)\left(N-1\right)\),
which is \(\frac{1}{2}N^2\) in the limit of large \(N\).
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.
[26]:
sim.dust.rhos
[26]:
[[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
[27]:
sim.dust.fill * sim.dust.rhos
[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]]
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 over 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}\)
[28]:
sim.dust.Sigma.sum(-1)
[28]:
[1.11815368e+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 1.24348839e-27 1.08303312e-27
9.43282410e-28 8.21564634e-28 7.15552883e-28 6.23220508e-28
5.42802371e-28 4.72761102e-28 4.11757706e-28 3.58625970e-28
3.12350162e-28 2.72045619e-28 2.36941829e-28 2.06367706e-28
1.79738758e-28 1.56545914e-28 1.36345791e-28 1.18752220e-28
1.03428860e-28 9.00827708e-29 7.84588135e-29 6.83347698e-29
5.95170965e-29 5.18372240e-29 4.51483347e-29 3.93225557e-29]
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 is 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 an 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 in that radial grid cell.
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}\)
Update order
The update order of the dust quantities in the standard model is set to
[29]:
sim.dust.updateorder
[29]:
['delta',
'rhos',
'fill',
'a',
'St',
'H',
'rho',
'backreaction',
'v',
'D',
'eps',
'kernel',
'p',
'S']
[30]:
sim.dust.backreaction.updateorder
[30]:
['A', 'B']
[31]:
sim.dust.delta.updateorder
[31]:
['rad', 'turb', 'vert']
[32]:
sim.dust.Fi.updateorder
[32]:
['adv', 'diff', 'tot']
[33]:
sim.dust.p.updateorder
[33]:
['frag', 'stick']
[34]:
sim.dust.S.updateorder
[34]:
['ext', 'tot']
[35]:
sim.dust.v.updateorder
[35]:
['frag', 'driftmax', 'rel']
[36]:
sim.dust.v.rel.updateorder
[36]:
['azi', 'brown', 'rad', 'turb', 'vert', 'tot']
Note: The quantities that are excluded are calculated in the finalization step of the implicit integrator from the new values of the dust surface densities.
Gas
[37]:
sim.gas
[37]:
Group (Gas quantities)
----------------------
boundary : Group (Boundary conditions)
S : Group (Source terms)
v : Group (Velocities)
-----
alpha : Field (Turbulent alpha parameter)
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
[38]:
sim.gas.boundary
[38]:
Group (Boundary conditions)
---------------------------
inner : Boundary
outer : Boundary
-----
These are the boundary conditions of the gas. By default the inner boundary is set to constant gradient and the outer boundary to floor value.
[39]:
sim.gas.boundary.inner
[39]:
Constant gradient
[40]:
sim.gas.boundary.outer
[40]:
Value
The boundary conditions can be modified with setcondition.
[41]:
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
[42]:
sim.gas.S
[42]:
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.
[43]:
sim.gas.S.ext
[43]:
[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
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
[44]:
sim.gas.v
[44]:
Group (Velocities)
------------------
rad : Field (Radial velocity [cm/s])
visc : Field (Viscous accretion velocity [cm/s])
-----
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}\).
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.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.
[45]:
sim.gas.alpha
[45]:
[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001
0.001 0.001 0.001 0.001 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
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.
[46]:
sim.gas.gamma
[46]:
[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}\).
[47]:
sim.gas.mu
[47]:
[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}\).
[48]:
sim.gas.SigmaFloor
[48]:
[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{1}{2}\frac{0.05\,L_*}{4\,\pi\,r^2\,\sigma_\mathrm{SB}}}\)
Update order
The update order of the gas quantities in the standard model is set to
[49]:
sim.gas.updateorder
[49]:
['gamma',
'mu',
'T',
'alpha',
'cs',
'Hp',
'nu',
'rho',
'n',
'mfp',
'P',
'eta',
'S']
[50]:
sim.gas.S.updateorder
[50]:
['ext', 'tot']
[51]:
sim.gas.v.updateorder
[51]:
['visc', 'rad']
Note: The quantities that are excluded are calculated in the finalization step of the implicit integrator from the new values of the surface densities.
Grid
[52]:
sim.grid
[52]:
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
-----
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.A
Surface area of the annulus a grid cell spans.
\(A_i = \pi \left( r_{i+\frac{1}{2}}^2 - r_{i-\frac{1}{2}}^2 \right)\)
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!
[53]:
sim.grid.m
[53]:
[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.
[54]:
sim.grid.Nm
[54]:
120
Simulation.grid.Nr
Number of radial grid cells.
[55]:
sim.grid.Nr
[55]:
100
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.
[56]:
sim.grid.r == 0.5 * (sim.grid.ri[:-1] + sim.grid.ri[1:])
[56]:
[ 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.
Update order
The update order of the grid quantities in the standard model is set to
[57]:
sim.grid.updateorder
[57]:
['OmegaK']
Star
[58]:
sim.star
[58]:
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.
[59]:
sim.star.M
[59]:
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.
[60]:
sim.star.R
[60]:
139140000000.0
Simulation.star.T
Stellar effective surface temperature. Initially set to the value given in Simulation.ini.star.T.
[61]:
sim.star.T
[61]:
5772.0
Update order
The update order of the stellar quantities in the standard model is set to
[62]:
sim.star.updateorder
[62]:
['M', 'R', 'T', 'L']
Time
Simulation.t is the current time and the integration variable. It starts at zero initially.
[63]:
sim.t
[63]:
0.0
Snapshots are written between \(10^3\) years and \(10^5\) years with 10 snapshots per time decade.
[64]:
sim.t.snapshots
[64]:
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])
The timestep is calculated from the dust and the gas source terms and the current values of the surface densities of dust and gas, while using a safety factor of 0.1, which can be accessed as an attribute of the integration variable.
[65]:
sim.t.cfl
[65]:
0.1
Integrator
The integrator that is used by default as two integration instructions. One for the gas and one for the dust.
[66]:
sim.integrator
[66]:
Integrator (Default integrator)
[67]:
sim.integrator.instructions
[67]:
[Instruction (Dust: implicit 1st-order direct solver),
Instruction (Gas: implicit 1st-order direct solver)]
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
[68]:
sim.gas.Sigma.jacobian() @ sim.gas.Sigma
[68]:
array([ 0.00000000e+00, -2.65695963e-11, -1.51742395e-11, -1.41277753e-11,
-1.31512275e-11, -1.22399343e-11, -1.13895454e-11, -1.05960011e-11,
-9.85551322e-12, -9.16454666e-12, -8.51980274e-12, -7.91820337e-12,
-7.35687637e-12, -6.83314177e-12, -6.34449895e-12, -5.88861478e-12,
-5.46331238e-12, -5.06656081e-12, -4.69646526e-12, -4.35125809e-12,
-4.02929031e-12, -3.72902371e-12, -3.44902350e-12, -3.18795143e-12,
-2.94455937e-12, -2.71768333e-12, -2.50623785e-12, -2.30921082e-12,
-2.12565856e-12, -1.95470127e-12, -1.79551880e-12, -1.64734665e-12,
-1.50947226e-12, -1.38123153e-12, -1.26200553e-12, -1.15121751e-12,
-1.04832998e-12, -9.52842068e-13, -8.64286959e-13, -7.82229559e-13,
-7.06264239e-13, -6.36012738e-13, -5.71122173e-13, -5.11263153e-13,
-4.56128002e-13, -4.05429063e-13, -3.58897093e-13, -3.16279733e-13,
-2.77340052e-13, -2.41855165e-13, -2.09614912e-13, -1.80420610e-13,
-1.54083866e-13, -1.30425466e-13, -1.09274328e-13, -9.04665397e-14,
-7.38444770e-14, -5.92560115e-14, -4.65538218e-14, -3.55948141e-14,
-2.62396618e-14, -1.83524760e-14, -1.18006141e-14, -6.45463314e-15,
-2.18839224e-15, 1.12069694e-15, 3.59127752e-15, 5.33752432e-15,
6.46850987e-15, 7.08748612e-15, 7.29111679e-15, 7.16870159e-15,
6.80143987e-15, 6.26178420e-15, 5.61293508e-15, 4.90852435e-15,
4.19252718e-15, 3.49943057e-15, 2.85467018e-15, 2.27532818e-15,
1.77106488e-15, 1.34523665e-15, 9.96136649e-16, 7.18283408e-16,
5.03678333e-16, 3.42957842e-16, 2.26378068e-16, 1.44589401e-16,
8.91814086e-17, 5.30026964e-17, 3.02814903e-17, 1.65881680e-17,
8.68880524e-18, 4.33878800e-18, 2.05888641e-18, 9.25247153e-19,
3.92320686e-19, 1.56337655e-19, 3.86150120e-20, 0.00000000e+00])
The total source terms, i.e., including external sources and excluding the boundaries are given by
[69]:
sim.gas.Sigma.jacobian() @ sim.gas.Sigma + sim.gas.S.ext
[69]:
[ 0.00000000e+00 -2.65695963e-11 -1.51742395e-11 -1.41277753e-11
-1.31512275e-11 -1.22399343e-11 -1.13895454e-11 -1.05960011e-11
-9.85551322e-12 -9.16454666e-12 -8.51980274e-12 -7.91820337e-12
-7.35687637e-12 -6.83314177e-12 -6.34449895e-12 -5.88861478e-12
-5.46331238e-12 -5.06656081e-12 -4.69646526e-12 -4.35125809e-12
-4.02929031e-12 -3.72902371e-12 -3.44902350e-12 -3.18795143e-12
-2.94455937e-12 -2.71768333e-12 -2.50623785e-12 -2.30921082e-12
-2.12565856e-12 -1.95470127e-12 -1.79551880e-12 -1.64734665e-12
-1.50947226e-12 -1.38123153e-12 -1.26200553e-12 -1.15121751e-12
-1.04832998e-12 -9.52842068e-13 -8.64286959e-13 -7.82229559e-13
-7.06264239e-13 -6.36012738e-13 -5.71122173e-13 -5.11263153e-13
-4.56128002e-13 -4.05429063e-13 -3.58897093e-13 -3.16279733e-13
-2.77340052e-13 -2.41855165e-13 -2.09614912e-13 -1.80420610e-13
-1.54083866e-13 -1.30425466e-13 -1.09274328e-13 -9.04665397e-14
-7.38444770e-14 -5.92560115e-14 -4.65538218e-14 -3.55948141e-14
-2.62396618e-14 -1.83524760e-14 -1.18006141e-14 -6.45463314e-15
-2.18839224e-15 1.12069694e-15 3.59127752e-15 5.33752432e-15
6.46850987e-15 7.08748612e-15 7.29111679e-15 7.16870159e-15
6.80143987e-15 6.26178420e-15 5.61293508e-15 4.90852435e-15
4.19252718e-15 3.49943057e-15 2.85467018e-15 2.27532818e-15
1.77106488e-15 1.34523665e-15 9.96136649e-16 7.18283408e-16
5.03678333e-16 3.42957842e-16 2.26378068e-16 1.44589401e-16
8.91814086e-17 5.30026964e-17 3.02814903e-17 1.65881680e-17
8.68880524e-18 4.33878800e-18 2.05888641e-18 9.25247153e-19
3.92320686e-19 1.56337655e-19 3.86150120e-20 0.00000000e+00]
The dust is integrated with an implicit first-order Euler integration scheme. The Jacobian is calculated from the parameters given in Simulation.dust. 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 hydrodynamic and coagulation source terms are given by
[70]:
(sim.dust.Sigma.jacobian() @ sim.dust.Sigma.ravel()).reshape(sim.dust.Sigma.shape)
[70]:
array([[ 6.69378160e-07, 7.07100367e-07, 7.46948374e-07, ...,
1.97891724e-30, 2.74969659e-30, 3.82069102e-30],
[-6.23483735e-08, -6.59511976e-08, -4.67032489e-08, ...,
-4.64083767e-34, -7.20390542e-34, -1.05662920e-33],
[-4.88452713e-08, -5.16639372e-08, -3.65877154e-08, ...,
-7.22801314e-34, -1.17572957e-33, -1.85635051e-33],
...,
[-8.55226776e-59, -1.06505134e-58, -1.32632479e-58, ...,
-1.18906666e-47, -1.48062233e-47, -1.84366661e-47],
[-2.59893623e-59, -3.23624508e-59, -4.02982080e-59, ...,
-3.61156047e-48, -4.49710457e-48, -5.59978150e-48],
[ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
Note the ravel() and reshape() operations, since the Jacobian works on a one-dimensional state vector.
The total source terms except at the boundaries are given by
[71]:
(sim.dust.Sigma.jacobian() @ sim.dust.Sigma.ravel()).reshape(sim.dust.Sigma.shape) + sim.dust.S.ext
[71]:
[[ 6.69378160e-07 7.07100367e-07 7.46948374e-07 ... 1.97891724e-30
2.74969659e-30 3.82069102e-30]
[-6.23483735e-08 -6.59511976e-08 -4.67032489e-08 ... -4.64083767e-34
-7.20390542e-34 -1.05662920e-33]
[-4.88452713e-08 -5.16639372e-08 -3.65877154e-08 ... -7.22801314e-34
-1.17572957e-33 -1.85635051e-33]
...
[-8.55226776e-59 -1.06505134e-58 -1.32632479e-58 ... -1.18906666e-47
-1.48062233e-47 -1.84366661e-47]
[-2.59893623e-59 -3.23624508e-59 -4.02982080e-59 ... -3.61156047e-48
-4.49710457e-48 -5.59978150e-48]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 ... 0.00000000e+00
0.00000000e+00 0.00000000e+00]]
At the end of a successful integration step, the floor values and boundaries are enforced.
Writer
DustPy uses by default the hdf5writer of simframe.
[72]:
sim.writer
[72]:
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.