#!/usr/bin/env python
# coding: utf-8

# # 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


# ## Dust

# In[5]:


sim.dust


# ### `Simulation.dust.backreaction`

# In[6]:


sim.dust.backreaction


# 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)](https://doi.org/10.3847/1538-4357/aaf4fc).
# 
# 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


# In[8]:


sim.dust.backreaction.B


# 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)](https://doi.org/10.3847/1538-4357/aaf4fc) for examples.

# ### `Simulation.dust.boundary`

# In[9]:


sim.dust.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


# In[11]:


sim.dust.boundary.outer


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


# ### `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


# 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


# $\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


# $\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


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

# ### `Simulation.dust.Fi`

# In[17]:


sim.dust.Fi


# 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


# #### `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


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


# #### `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


# 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


# #### `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


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

# #### `Simulation.dust.v.rel`

# In[24]:


sim.dust.v.rel


# These are the different sources of relative particle collision velocities used by the coagulation algorithm.  
# They are used for two different reasons:
# 1. To calculated the collision rates
# 2. 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)](https://doi.org/10.1051/0004-6361:20066899).  
# 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)](https://doi.org/10.1051/0004-6361/200913731) and follows [Dullemond & Dominik (2004)](https://doi.org/10.1051/0004-6361:20040284).

# #### `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


# #### `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


# ### `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)](https://doi.org/10.1006/icar.1995.1058)
# 
# $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)](https://doi.org/10.1051/0004-6361/200913731) 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


# To calculate the mass density of a dust aggregate do

# In[28]:


sim.dust.fill * sim.dust.rhos


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


# 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


# ### `Simulation.gas.boundary`

# In[31]:


sim.gas.boundary


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

# In[32]:


sim.gas.boundary.inner


# In[33]:


sim.gas.boundary.outer


# The boundary conditions can be modified with `setcondition`.

# In[34]:


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


# 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


# 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


# #### `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


# 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)](https://ui.adsabs.harvard.edu/abs/1973A%26A....24..337S/abstract). It is initially set to the value in `Simulation.ini.gas.alpha`.

# In[38]:


sim.gas.alpha


# ### `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


# 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


# ### `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


# 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


# 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


# ### `Simulation.grid.Nm`
# 
# Number of mass bins.

# In[44]:


sim.grid.Nm


# ### `Simulation.grid.Nr`
# 
# Number of radial grid cells.

# In[45]:


sim.grid.Nr


# ### `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:])


# ### `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


# ### `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


# ### `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


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

# In[50]:


sim.star.T


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

# In[51]:


sim.t


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

# In[52]:


sim.t.snapshots


# 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


# In[54]:


sim.t.stepsize


# ## 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


# In[56]:


sim.integrator.instructions


# 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


# 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


# 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


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