{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Standard Model\n",
    "\n",
    "In this chapter we'll go through all fields in the standard model, explain their meaning and the default functions that calculates them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dustpy import Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = Simulation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim.initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DustPy\n",
       "------\n",
       "    dust         : Group (Dust quantities)\n",
       "    gas          : Group (Gas quantities)\n",
       "    grid         : Group (Grid quantities)\n",
       "    star         : Group (Stellar quantities)\n",
       "  -----\n",
       "    t            : Field (Time [s]), \u001b[95mIntegration variable\u001b[0m\n",
       "  -----\n",
       "    Integrator   : Integrator (Default integrator)\n",
       "    Writer       : Writer (HDF5 file format using h5py)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dust"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Dust quantities)\n",
       "-----------------------\n",
       "    backreaction : Group (Backreaction coefficients)\n",
       "    boundary     : Group (Boundary conditions)\n",
       "    coagulation  : Group (Coagulation quantities)\n",
       "    delta        : Group (Mixing parameters)\n",
       "    Fi           : Group (Fluxes)\n",
       "    p            : Group (Probabilities)\n",
       "    S            : Group (Sources)\n",
       "    v            : Group (Velocities)\n",
       "  -----\n",
       "    a            : Field (Particle size [cm)\n",
       "    D            : Field (Diffusivity [cm²/s])\n",
       "    eps          : Field (Dust-to-gas ratio)\n",
       "    fill         : Field (Filling factor)\n",
       "    H            : Field (Scale heights [cm])\n",
       "    kernel       : Field (Collision kernel [cm²/s])\n",
       "    rho          : Field (Midplane mass density per mass bin [g/cm³])\n",
       "    rhos         : Field (Solid state density [g/cm³])\n",
       "    Sigma        : Field (Surface density per mass bin [g/cm²])\n",
       "    SigmaFloor   : Field (Floor value of surface density [g/cm²])\n",
       "    St           : Field (Stokes number)\n",
       "  -----"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.backreaction`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Backreaction coefficients)\n",
       "---------------------------------\n",
       "    A            : Field (Pull factor)\n",
       "    B            : Field (Push factor)\n",
       "  -----"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.backreaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "The details of this mechanism are described in [Gárate et al. (2019)](https://doi.org/10.3847/1538-4357/aaf4fc).\n",
    "\n",
    "Backreaction modifies the radial gas velocity as follows\n",
    "\n",
    "$v_\\mathrm{g} = Av_\\mathrm{visc} + 2B\\eta v_\\mathrm{K}$.\n",
    "\n",
    "In the standard model we have `A=1` and `B=0` everywhere, i.e., backreactions is not active."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
       " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
       " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
       " 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.\n",
       " 1. 1. 1. 1.]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.backreaction.A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0.]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.backreaction.B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.boundary`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Boundary conditions)\n",
       "---------------------------\n",
       "    inner        : Boundary\n",
       "    outer        : Boundary\n",
       "  -----"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default the inner dust boundary is set to constant gradient and the outer boundary to floor value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Constant gradient"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.boundary.inner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Value"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.boundary.outer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The outer boundary therefore stores an array with the floor values of all particle masses.\n",
    "\n",
    "The boundary conditions can be changed via `setcondition()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method setcondition in module dustpy.utils.boundary:\n",
      "\n",
      "setcondition(condition, value=None) method of dustpy.utils.boundary.Boundary instance\n",
      "    Function to set boundary condition.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    condition : string\n",
      "        Type of boundary conditon:\n",
      "            - \"const_grad\" : constant gradient\n",
      "            - \"const_pow\" : constant power law\n",
      "            - \"const_val\" : constant value\n",
      "            - \"val\" : custom value\n",
      "            - \"grad\" : custom gradient\n",
      "            - \"pow\" : custom power law with set exponent\n",
      "            - None : Don't impose boundary condition (default)\n",
      "    value : float or array, optional, default : None\n",
      "        Value if needed for boundary condition\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(sim.dust.boundary.outer.setcondition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.coagulation`\n",
    "\n",
    "The fields in this group define the behavior of dust growth and are discussed in a separate chapter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.delta`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Mixing parameters)\n",
       "-------------------------\n",
       "    rad          : Field (Radial mixing parameter)\n",
       "    turb         : Field (Turbulent mixing parameter)\n",
       "    vert         : Field (Vertical mixing parameter)\n",
       "  -----"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.delta.rad`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.delta.rad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\delta_\\mathrm{rad}$ will be used to calculate the radial, turbulent RMS velocity of the dust."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.delta.turb`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.delta.turb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\delta_\\mathrm{turb}$ will be used to calculate the turbulent collision velocities of the dust particles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.delta.vert`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.delta.vert"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\delta_\\mathrm{vert}$ will be used to calculate the vertical dust scale heights."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.Fi`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Fluxes)\n",
       "--------------\n",
       "    adv          : Field (Advective flux [g/cm/s)\n",
       "    diff         : Field (Diffusive flux [g/cm/s)\n",
       "    tot          : Field (Total flux [g/cm/s)\n",
       "  -----"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.Fi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a group of fluxes through the radial grid interfaces for every particle mass of shape `(Simulation.grid.Nr+1,, Simulation.grid.Nm)`.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.Fi.adv`\n",
    "\n",
    "This is the advective flux calculated by $F_\\mathrm{adv} = v_\\mathrm{d}\\Sigma_\\mathrm{d}$.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.Fi.diff`\n",
    "\n",
    "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.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  0.00000000e+00\n",
       "   0.00000000e+00  0.00000000e+00]\n",
       " [-0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  0.00000000e+00\n",
       "   0.00000000e+00  0.00000000e+00]\n",
       " [-0.00000000e+00 -0.00000000e+00 -0.00000000e+00 ...  8.42570438e-24\n",
       "   1.14312063e-23  1.53229589e-23]\n",
       " ...\n",
       " [-3.21920964e-49 -3.59231024e-49 -4.00864247e-49 ... -1.20006815e-43\n",
       "  -1.33913684e-43 -1.49432136e-43]\n",
       " [-0.00000000e+00 -0.00000000e+00 -0.00000000e+00 ... -0.00000000e+00\n",
       "  -0.00000000e+00 -0.00000000e+00]\n",
       " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  0.00000000e+00\n",
       "   0.00000000e+00  0.00000000e+00]]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.Fi.diff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.Fi.tot`\n",
    "\n",
    "This is the total flux through the radial grid interfaces $F_\\mathrm{tot} = F_\\mathrm{adv} + F_\\mathrm{diff}$.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.p`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Probabilities)\n",
       "---------------------\n",
       "    frag         : Field (Fragmentation probability)\n",
       "    stick        : Field (Sticking probability)\n",
       "  -----"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the fragmentation and sticking probability of certain particle collisions.  \n",
    "Their shape is `(Simulation.grid.Nr, Simulation.grid.Nm, Simulation.grid.Nm)`.\n",
    "\n",
    "The fragmentation probability of particle `i=80`, with particle `j=3` at radial grid cell `ir=30` is given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ir = 30\n",
    "i = 80\n",
    "j = 3\n",
    "sim.dust.p.frag[ir, j, i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.p.frag`\n",
    "\n",
    "The fragmentation probability has a smooth transition between no fragmentation and fragmentation at the fragmentation velocity.\n",
    "\n",
    "$p_\\mathrm{f} = \n",
    "\\begin{cases}\n",
    "    1, & \\text{if } v_\\mathrm{rel}>v_\\mathrm{frag}\\\\\n",
    "    0, & \\text{if } v_\\mathrm{rel}<v_\\mathrm{frag} - \\Delta v\\\\\n",
    "    \\frac{v_\\mathrm{rel} - v_\\mathrm{frag} + \\Delta v}{\\Delta v}, & \\text{otherwise}\n",
    "\\end{cases}$\n",
    "\n",
    "with $\\Delta v = 0.2\\,v_\\mathrm{frag}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.p.stick`\n",
    "\n",
    "This is the sticking probability given by $p_\\mathrm{s} = 1 - p_\\mathrm{f}$.\n",
    "\n",
    "Bouncing can be easily implemented if $0 \\leq p_\\mathrm{f} + p_\\mathrm{s} < 1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.S`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Sources)\n",
       "---------------\n",
       "    coag         : Field (Coagulation sources [g/cm²/s])\n",
       "    ext          : Field (External sources [g/cm²/s])\n",
       "    hyd          : Field (Hydrodynamic sources [g/cm²/s])\n",
       "    tot          : Field (Tot sources [g/cm²/s])\n",
       "  -----"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.S.coag`\n",
    "\n",
    "The source terms from dust coagulation. This is described in detail in a separate chapter.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.S.ext`\n",
    "\n",
    "External source terms. These are by default set to zero, i.e. no external sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[0. 0. 0. ... 0. 0. 0.]\n",
       " [0. 0. 0. ... 0. 0. 0.]\n",
       " [0. 0. 0. ... 0. 0. 0.]\n",
       " ...\n",
       " [0. 0. 0. ... 0. 0. 0.]\n",
       " [0. 0. 0. ... 0. 0. 0.]\n",
       " [0. 0. 0. ... 0. 0. 0.]]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.S.ext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.S.hyd`\n",
    "\n",
    "These are the hydrodynamic source terms. At grid cell $i$ they are calculated by\n",
    "\n",
    "$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}$.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.S.tot`\n",
    "\n",
    "These are the total source terms given by $S_\\mathrm{tot} = S_\\mathrm{coag} + S_\\mathrm{ext} + S_\\mathrm{hyd}$.\n",
    "\n",
    "These are the actual sources the dust integrator uses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.v`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Velocities)\n",
       "------------------\n",
       "    rel          : Group (Relative velocities)\n",
       "  -----\n",
       "    driftmax     : Field (Maximum drift velocity [cm/s])\n",
       "    frag         : Field (Fragmentation velocity [cm/s])\n",
       "    rad          : Field (Radial velocity [cm/s])\n",
       "  -----"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are some dust related velocities the simulation needs for execution. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.v.rel`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Relative velocities)\n",
       "---------------------------\n",
       "    azi          : Field (Relative azimuthal velocity [cm/s])\n",
       "    brown        : Field (Relative Brownian motion velocity [cm/s])\n",
       "    rad          : Field (Relative radial velocity [cm/s])\n",
       "    tot          : Field (Total relative velocity [cm/s])\n",
       "    turb         : Field (Relative turbulent velocity [cm/s])\n",
       "    vert         : Field (Relative vertical settling velocity [cm/s])\n",
       "  -----"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.v.rel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the different sources of relative particle collision velocities used by the coagulation algorithm.  \n",
    "They are used for two different reasons:\n",
    "1. To calculated the collision rates\n",
    "2. To determine the outcome of a collision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `Simulation.dust.v.rel.azi`\n",
    "\n",
    "Relative collision velocity caused by a differential azimuthal drift of particles with different Stokes numbers calculated by\n",
    "\n",
    "$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|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `Simulation.dust.v.brown`\n",
    "\n",
    "Relative collision velocity of particles caused by Brownian motion calculated with\n",
    "\n",
    "$v_{\\mathrm{rel},\\,\\mathrm{brown}} = \\frac{8k_\\mathrm{B}T\\left(m_i + m_j \\right)}{\\pi m_i m_j}$.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `Simulation.dust.v.rel.rad`\n",
    "\n",
    "Relative collision velocity caused by differential radial drift.\n",
    "\n",
    "$v_{\\mathrm{rel},\\,\\mathrm{rad}} = \\left| v_{\\mathrm{rad},\\,i} - v_{\\mathrm{rad},\\,j} \\right|$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `Simulation.dust.v.rel.tot`\n",
    "\n",
    "Total relative velocities calculated by using the root mean square of all individual velocity sources.\n",
    "\n",
    "$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}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `Simulation.dust.v.rel.turb`\n",
    "\n",
    "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).  \n",
    "It uses `Simulation.dust.delta.turb` instead of `Simulation.gas.alpha` to calculate the velocities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### `Simulation.dust.v.rel.vert`\n",
    "\n",
    "Relative collision velocities caused by differential vertical settling of particles.\n",
    "\n",
    "$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}$.\n",
    "\n",
    "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)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.v.driftmax`\n",
    "\n",
    "This is the maximum drift velocity a particle of $\\mathrm{St} = 1$ can have.\n",
    "\n",
    "$v_{\\mathrm{drift},\\,\\mathrm{max}} = \\frac{1}{2} B v_\\mathrm{visc} - A \\eta v_\\mathrm{K}$.\n",
    "\n",
    "See `Simulation.dust.backreaction` for details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.v.frag`\n",
    "\n",
    "Fragmentation velocities of shape `(Simulation.grid.Nr,)`. By default this is set by the value of `Simulation.ini.dust.vfrag`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100. 100.\n",
       " 100. 100.]"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.v.frag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.dust.v.rad`\n",
    "\n",
    "Radial velocities of the dust particles.\n",
    "\n",
    "$v_\\mathrm{rad} = \\left( v_\\mathrm{g} + 2 v_{\\mathrm{drift}\\,\\mathrm{max}} \\right) \\frac{\\mathrm{St}}{\\mathrm{St}^2+1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.a`\n",
    "\n",
    "Particle radii calculated with\n",
    "\n",
    "$a = \\sqrt[3]{\\frac{3m}{4\\pi\\rho_\\mathrm{s}}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.D`\n",
    "\n",
    "Dust diffusivity for every particle species of shape `(Simulation.grid.Nr, Simulation.grid.Nm)` calculated with\n",
    "\n",
    "$D = \\frac{\\delta_\\mathrm{rad}\\,c_\\mathrm{s}^2}{\\Omega_\\mathrm{K}\\left( 1 + \\mathrm{St}^2 \\right)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.eps`\n",
    "\n",
    "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\n",
    "\n",
    "$\\epsilon = \\frac{\\sum\\limits_i \\Sigma_{\\mathrm{d},\\,i}}{\\Sigma_\\mathrm{g}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.fill`\n",
    "\n",
    "This describes the filling factor of the dust aggregates. By default this is 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1. 1. 1. ... 1. 1. 1.]\n",
       " [1. 1. 1. ... 1. 1. 1.]\n",
       " [1. 1. 1. ... 1. 1. 1.]\n",
       " ...\n",
       " [1. 1. 1. ... 1. 1. 1.]\n",
       " [1. 1. 1. ... 1. 1. 1.]\n",
       " [1. 1. 1. ... 1. 1. 1.]]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.fill"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.H`\n",
    "\n",
    "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)\n",
    "\n",
    "$H_i = H_\\mathrm{P} \\cdot \\sqrt{\\frac{\\delta_\\mathrm{vert}}{\\delta_\\mathrm{vert}+\\mathrm{St}_i}}$\n",
    "\n",
    "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`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.kernel`\n",
    "\n",
    "These are the vertically integrated collision kernels of shape `(Simulation.grid.Nr, Simulation.grid.Nm, Simulation.grid.Nm)` calculated with\n",
    "\n",
    "$K = \\frac{\\pi\\left( a_i + a_j \\right)^2}{\\sqrt{2\\pi\\left( H_i^2 + H_j^2 \\right)}} v_\\mathrm{rel}$\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.rho`\n",
    "\n",
    "Midplane mass densities of the dust.\n",
    "\n",
    "$\\rho = \\frac{\\Sigma}{\\sqrt{2\\pi}H}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.rhos`\n",
    "\n",
    "Solid state density of the particle material. This is initially set to `Simulation.ini.dust.rhoMonomer`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " ...\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]]"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.rhos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate the mass density of a dust aggregate do"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " ...\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]\n",
       " [1.67 1.67 1.67 ... 1.67 1.67 1.67]]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.fill * sim.dust.rhos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the density that is used to calculate the particle size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.Sigma`\n",
    "\n",
    "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.\n",
    "\n",
    "The surface densities are integrated of the mass bin, so a numerical summation over the mass dimension returns the total dust surface density.\n",
    "\n",
    "$\\Sigma_\\mathrm{d} = \\sum\\limits_i \\Sigma_{\\mathrm{d},\\,i}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.11815368e+01 1.04232155e+01 9.71465059e+00 9.05339692e+00\n",
       " 8.43629793e+00 7.86040780e+00 7.32297742e+00 6.82144128e+00\n",
       " 6.35340519e+00 5.91663488e+00 5.50904531e+00 5.12869076e+00\n",
       " 4.77375546e+00 4.44254502e+00 4.13347826e+00 3.84507969e+00\n",
       " 3.57597248e+00 3.32487186e+00 3.09057900e+00 2.87197527e+00\n",
       " 2.66801690e+00 2.47773001e+00 2.30020595e+00 2.13459693e+00\n",
       " 1.98011200e+00 1.83601326e+00 1.70161231e+00 1.57626698e+00\n",
       " 1.45937826e+00 1.35038739e+00 1.24877321e+00 1.15404966e+00\n",
       " 1.06576343e+00 9.83491798e-01 9.06840556e-01 8.35442144e-01\n",
       " 7.68953850e-01 7.07056150e-01 6.49451154e-01 5.95861146e-01\n",
       " 5.46027226e-01 4.99708026e-01 4.56678518e-01 4.16728886e-01\n",
       " 3.79663471e-01 3.45299782e-01 3.13467557e-01 2.84007884e-01\n",
       " 2.56772375e-01 2.31622377e-01 2.08428234e-01 1.87068586e-01\n",
       " 1.67429708e-01 1.49404882e-01 1.32893806e-01 1.17802035e-01\n",
       " 1.04040457e-01 9.15247991e-02 8.01751688e-02 6.99156307e-02\n",
       " 6.06738170e-02 5.23805770e-02 4.49696665e-02 3.83774784e-02\n",
       " 3.25428186e-02 2.74067268e-02 2.29123459e-02 1.90048396e-02\n",
       " 1.56313572e-02 1.27410459e-02 1.02851052e-02 8.21687989e-03\n",
       " 6.49198428e-03 5.06844837e-03 3.90687641e-03 2.97060545e-03\n",
       " 2.22585132e-03 8.21564634e-28 7.15552883e-28 6.23220508e-28\n",
       " 5.42802371e-28 4.72761102e-28 4.11757706e-28 3.58625970e-28\n",
       " 3.12350162e-28 2.72045619e-28 2.36941829e-28 2.06367706e-28\n",
       " 1.79738758e-28 1.56545914e-28 1.36345791e-28 1.18752220e-28\n",
       " 1.03428860e-28 9.00827708e-29 7.84588135e-29 6.83347698e-29\n",
       " 5.95170965e-29 5.18372240e-29 4.51483347e-29 3.93225557e-29]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.Sigma.sum(-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.SigmaFloor`\n",
    "\n",
    "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.\n",
    "\n",
    "$\\Sigma_{\\mathrm{d},\\,\\mathrm{floor}} = \\frac{m}{A_\\mathrm{annulus}}$\n",
    "\n",
    "Densities below the floor value therefore correspond to fewer that one physical particle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.dust.St`\n",
    "\n",
    "Stokes numbers of the particles. By default the Epstein and the Stokes I drag regime are considered.\n",
    "\n",
    "$\\mathrm{St} =\n",
    "\\begin{cases}\n",
    "\\frac{\\pi}{2} \\frac{a\\rho}{\\Sigma_\\mathrm{g}}, & \\text{if } a < \\frac{9}{4} \\lambda_\\mathrm{mfp}\\\\\n",
    "\\frac{2\\pi}{9} \\frac{a^2 \\rho}{\\lambda_\\mathrm{mfp} \\Sigma_\\mathrm{g}}, & \\text{else}\n",
    "\\end{cases}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update order\n",
    "\n",
    "The update order of the dust quantities in the standard model is set to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['delta',\n",
       " 'rhos',\n",
       " 'fill',\n",
       " 'a',\n",
       " 'St',\n",
       " 'H',\n",
       " 'rho',\n",
       " 'backreaction',\n",
       " 'v',\n",
       " 'D',\n",
       " 'eps',\n",
       " 'kernel',\n",
       " 'p',\n",
       " 'S']"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['A', 'B']"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.backreaction.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['rad', 'turb', 'vert']"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.delta.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['adv', 'diff', 'tot']"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.Fi.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['frag', 'stick']"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.p.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ext', 'tot']"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.S.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['frag', 'driftmax', 'rel']"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.v.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['azi', 'brown', 'rad', 'turb', 'vert', 'tot']"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.dust.v.rel.updateorder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: The quantities that are excluded are calculated in the finalization step of the implicit integrator from the new values of the surface densities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Gas quantities)\n",
       "----------------------\n",
       "    boundary     : Group (Boundary conditions)\n",
       "    S            : Group (Source terms)\n",
       "    v            : Group (Velocities)\n",
       "  -----\n",
       "    alpha        : Field (Turbulent alpha parameters)\n",
       "    cs           : Field (Sound speed [cm/s])\n",
       "    eta          : Field (Pressure gradient parameter)\n",
       "    Fi           : Field (Gas flux interfaces [g/cm/s])\n",
       "    gamma        : Field (Adiabatic index)\n",
       "    Hp           : Field (Pressure scale height [cm])\n",
       "    mfp          : Field (Midplane mean free path [cm])\n",
       "    mu           : Field (Mean molecular weight [g])\n",
       "    n            : Field (Miplane number density [1/cm³])\n",
       "    nu           : Field (Kinematic viscosity [cm²/s])\n",
       "    P            : Field (Midplane pressure [g/cm/s²])\n",
       "    rho          : Field (Miplane mass density [g/cm³])\n",
       "    Sigma        : Field (Surface density [g/cm²])\n",
       "    SigmaFloor   : Field (Floor value of surface density [g/cm²])\n",
       "    T            : Field (Temperature [K])\n",
       "  -----"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.boundary`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Boundary conditions)\n",
       "---------------------------\n",
       "    inner        : Boundary\n",
       "    outer        : Boundary\n",
       "  -----"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the boundary conditions of the gas. They are set to constant gradient by default."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Constant gradient"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.boundary.inner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Value"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.boundary.outer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary conditions can be modified with `setcondition`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method setcondition in module dustpy.utils.boundary:\n",
      "\n",
      "setcondition(condition, value=None) method of dustpy.utils.boundary.Boundary instance\n",
      "    Function to set boundary condition.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    condition : string\n",
      "        Type of boundary conditon:\n",
      "            - \"const_grad\" : constant gradient\n",
      "            - \"const_pow\" : constant power law\n",
      "            - \"const_val\" : constant value\n",
      "            - \"val\" : custom value\n",
      "            - \"grad\" : custom gradient\n",
      "            - \"pow\" : custom power law with set exponent\n",
      "            - None : Don't impose boundary condition (default)\n",
      "    value : float or array, optional, default : None\n",
      "        Value if needed for boundary condition\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(sim.gas.boundary.inner.setcondition)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.S`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Source terms)\n",
       "--------------------\n",
       "    ext          : Field (External sources [g/cm²/s])\n",
       "    hyd          : Field (Hydrodynamic sources [g/cm²/s])\n",
       "    tot          : Field (Total sources [g/cm²/s])\n",
       "  -----"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.S"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the source terms of the gas.\n",
    "\n",
    "`Simulation.gas.S` has by default a distole that calculates the implicit fluxes through the boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.gas.S.ext`\n",
    "\n",
    "These are the external source terms for gas evolution, e.g. infall. By default these are set to zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.\n",
       " 0. 0. 0. 0.]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.S.ext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.gas.S.hyd`\n",
    "\n",
    "These are the hydrodynamic source terms of the gas evolution.  \n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.gas.S.tot`\n",
    "\n",
    "These are the total source terms of gas evolution.\n",
    "\n",
    "$S_\\mathrm{tot} = S_\\mathrm{ext} + S_\\mathrm{hyd}$\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.v`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Velocities)\n",
       "------------------\n",
       "    rad          : Field (Radial velocity [cm/s])\n",
       "    visc         : Field (Viscous accretion velocity [cm/s])\n",
       "  -----"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are velocities that are relevant for the gas evolution.  \n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.gas.v.rad`\n",
    "\n",
    "This is the radial gas velocity. It is given by\n",
    "\n",
    "$v_\\mathrm{g} = Av_\\mathrm{visc} + 2B\\eta v_\\mathrm{K}$.\n",
    "\n",
    "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.  \n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `Simulation.gas.v.visc`\n",
    "\n",
    "This is the radial viscous gas velocity\n",
    "\n",
    "$v_\\mathrm{visc} = -\\frac{3}{\\Sigma_\\mathrm{g}\\sqrt{R}} \\frac{\\partial}{\\partial R} \\left( \\Sigma_\\mathrm{g} \\nu \\sqrt{R} \\right)$\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.alpha`\n",
    "\n",
    "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`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001\n",
       " 0.001 0.001 0.001 0.001]"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.cs`\n",
    "\n",
    "This is the adiabatic sound speed in the midplane of the disk\n",
    "\n",
    "$c_\\mathrm{s} = \\sqrt{\\frac{\\gamma k_\\mathrm{B} T}{\\mu}}$.\n",
    "\n",
    "For isothermal simulations set $\\gamma=1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.eta`\n",
    "\n",
    "This is the midplane pressure gradient parameter $\\eta$ given by\n",
    "\n",
    "$\\eta = -\\frac{1}{2} \\left( \\frac{H_\\mathrm{P}}{r} \\right)^2 \\frac{\\partial \\log P}{\\partial \\log r}$\n",
    "\n",
    "It describes the degree of \"sub-Keplerity\" of the disk\n",
    "\n",
    "$v_\\phi^2 = \\left( 1-2\\eta \\right) v_\\mathrm{K}^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.Fi`\n",
    "\n",
    "These are the mass fluxes of gas through the grid cell interfaces.  \n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.gamma`\n",
    "\n",
    "This is the ratio of specific heats. It is initially set to the value given by `Simulation.ini.gas.gamma`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4 1.4]"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.gamma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For isothermal simulations set this to $1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.Hp`\n",
    "\n",
    "Pressure scale height of the gas given by the ratio of the isothermal sound speed to the Keplerian frequency.\n",
    "\n",
    "$H_\\mathrm{P} = \\frac{c_{\\mathrm{s},\\,\\mathrm{iso}}}{\\Omega_\\mathrm{K}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.mfp`\n",
    "\n",
    "Mean free path of the gas in the midplane of the disk\n",
    "\n",
    "$\\lambda_\\mathrm{mfp} = \\frac{1}{\\sqrt{2}\\,n\\,\\sigma_\\mathrm{H_2}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.mu`\n",
    "\n",
    "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}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24\n",
       " 3.84703042e-24 3.84703042e-24 3.84703042e-24 3.84703042e-24]"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.n`\n",
    "\n",
    "Midplane number density of the gas given by\n",
    "\n",
    "$n = \\frac{\\rho}{\\mu}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.nu`\n",
    "\n",
    "Kinematic viscosity of the gas given by\n",
    "\n",
    "$\\nu = \\alpha c_\\mathrm{s} H_\\mathrm{P}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.P`\n",
    "\n",
    "Midplane gas pressure given by\n",
    "\n",
    "$P = \\frac{\\rho\\,c_\\mathrm{s}^2}{\\gamma}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.rho`\n",
    "\n",
    "Midplane gas mass density given by\n",
    "\n",
    "$\\rho = \\frac{\\Sigma_\\mathrm{g}}{\\sqrt{2\\pi}H_\\mathrm{P}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.Sigma`\n",
    "\n",
    "Gas surface density. This is the quantity that is integrated with an implicit Euler first-order scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.SigmaFloor`\n",
    "\n",
    "This is the floor value of the gas surface density. By default it is $10^{-100}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 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\n",
       " 1.e-100]"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.SigmaFloor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.gas.T`\n",
    "\n",
    "This is the midplane gas temperature. It is calculated by assuming a passively irradiated disk with a constant irradiation angle of $0.05$.\n",
    "\n",
    "$T\\left( r \\right) = \\sqrt[4]{\\frac{0.05\\,L_*}{4\\,\\pi\\,r^2\\,\\sigma_\\mathrm{SB}}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update order\n",
    "\n",
    "The update order of the gas quantities in the standard model is set to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['gamma',\n",
       " 'mu',\n",
       " 'T',\n",
       " 'alpha',\n",
       " 'cs',\n",
       " 'Hp',\n",
       " 'nu',\n",
       " 'rho',\n",
       " 'n',\n",
       " 'mfp',\n",
       " 'P',\n",
       " 'eta',\n",
       " 'S']"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['ext', 'tot']"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.S.updateorder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['visc', 'rad']"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.v.updateorder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note**: The quantities that are excluded are calculated in the finalization step of the implicit integrator from the new values of the surface densities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Grid quantities)\n",
       "-----------------------\n",
       "    A            : Field (Radial grid annulus area [cm²]), \u001b[95mconstant\u001b[0m\n",
       "    m            : Field (Mass grid [g]), \u001b[95mconstant\u001b[0m\n",
       "    Nm           : Field (# of mass bins), \u001b[95mconstant\u001b[0m\n",
       "    Nr           : Field (# of radial grid cells), \u001b[95mconstant\u001b[0m\n",
       "    OmegaK       : Field (Keplerian frequency [1/s])\n",
       "    r            : Field (Radial grid cell centers [cm]), \u001b[95mconstant\u001b[0m\n",
       "    ri           : Field (Radial grid cell interfaces [cm]), \u001b[95mconstant\u001b[0m\n",
       "  -----"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.grid.m`\n",
    "\n",
    "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!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.00000000e-12 1.38949549e-12 1.93069773e-12 2.68269580e-12\n",
       " 3.72759372e-12 5.17947468e-12 7.19685673e-12 1.00000000e-11\n",
       " 1.38949549e-11 1.93069773e-11 2.68269580e-11 3.72759372e-11\n",
       " 5.17947468e-11 7.19685673e-11 1.00000000e-10 1.38949549e-10\n",
       " 1.93069773e-10 2.68269580e-10 3.72759372e-10 5.17947468e-10\n",
       " 7.19685673e-10 1.00000000e-09 1.38949549e-09 1.93069773e-09\n",
       " 2.68269580e-09 3.72759372e-09 5.17947468e-09 7.19685673e-09\n",
       " 1.00000000e-08 1.38949549e-08 1.93069773e-08 2.68269580e-08\n",
       " 3.72759372e-08 5.17947468e-08 7.19685673e-08 1.00000000e-07\n",
       " 1.38949549e-07 1.93069773e-07 2.68269580e-07 3.72759372e-07\n",
       " 5.17947468e-07 7.19685673e-07 1.00000000e-06 1.38949549e-06\n",
       " 1.93069773e-06 2.68269580e-06 3.72759372e-06 5.17947468e-06\n",
       " 7.19685673e-06 1.00000000e-05 1.38949549e-05 1.93069773e-05\n",
       " 2.68269580e-05 3.72759372e-05 5.17947468e-05 7.19685673e-05\n",
       " 1.00000000e-04 1.38949549e-04 1.93069773e-04 2.68269580e-04\n",
       " 3.72759372e-04 5.17947468e-04 7.19685673e-04 1.00000000e-03\n",
       " 1.38949549e-03 1.93069773e-03 2.68269580e-03 3.72759372e-03\n",
       " 5.17947468e-03 7.19685673e-03 1.00000000e-02 1.38949549e-02\n",
       " 1.93069773e-02 2.68269580e-02 3.72759372e-02 5.17947468e-02\n",
       " 7.19685673e-02 1.00000000e-01 1.38949549e-01 1.93069773e-01\n",
       " 2.68269580e-01 3.72759372e-01 5.17947468e-01 7.19685673e-01\n",
       " 1.00000000e+00 1.38949549e+00 1.93069773e+00 2.68269580e+00\n",
       " 3.72759372e+00 5.17947468e+00 7.19685673e+00 1.00000000e+01\n",
       " 1.38949549e+01 1.93069773e+01 2.68269580e+01 3.72759372e+01\n",
       " 5.17947468e+01 7.19685673e+01 1.00000000e+02 1.38949549e+02\n",
       " 1.93069773e+02 2.68269580e+02 3.72759372e+02 5.17947468e+02\n",
       " 7.19685673e+02 1.00000000e+03 1.38949549e+03 1.93069773e+03\n",
       " 2.68269580e+03 3.72759372e+03 5.17947468e+03 7.19685673e+03\n",
       " 1.00000000e+04 1.38949549e+04 1.93069773e+04 2.68269580e+04\n",
       " 3.72759372e+04 5.17947468e+04 7.19685673e+04 1.00000000e+05]"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.grid.m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.grid.Nm`\n",
    "\n",
    "Number of mass bins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "120"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.grid.Nm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.grid.Nr`\n",
    "\n",
    "Number of radial grid cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "100"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.grid.Nr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.grid.OmegaK`\n",
    "\n",
    "Keplerian frequency given by\n",
    "\n",
    "$\\Omega_\\mathrm{K} = \\sqrt{\\frac{G\\,M_*}{r^3}}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.grid.r`\n",
    "\n",
    "Radial grid cell centers. The radial grid cell centers are exactly in the middle between the grid cell interfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True  True  True  True  True  True  True  True  True\n",
       "  True  True  True  True]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.grid.r == 0.5 * (sim.grid.ri[:-1] + sim.grid.ri[1:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.grid.ri`\n",
    "\n",
    "Locations of the grid cell interfaces. `Simulation.grid.ri[0]` and `Simulation.grid.ri[-1]` are the inner and outer grid boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update order\n",
    "\n",
    "The update order of the grid quantities in the standard model is set to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['OmegaK']"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.grid.updateorder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Group (Stellar quantities)\n",
       "--------------------------\n",
       "    L            : Field (Luminosity [erg/s])\n",
       "    M            : Field (Mass [g])\n",
       "    R            : Field (Radius [cm])\n",
       "    T            : Field (Effective temperature [K])\n",
       "  -----"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.star"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.star.L`\n",
    "\n",
    "Stellar luminosity given by\n",
    "\n",
    "$L = 4\\pi\\,R_*^2\\,\\sigma_\\mathrm{SB}\\,T_*^4$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.star.M`\n",
    "\n",
    "Stellar mass. Initially set to the value given in `Simulation.ini.star.M`, which corresponds to one Solar mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.988409870698051e+33"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.star.M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.star.R`\n",
    "\n",
    "Stellar radius. Initially set to the value given in `Simulation.ini.star.R`, which corresponds to 2 Solar radii."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "139140000000.0"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.star.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `Simulation.star.T`\n",
    "\n",
    "Stellar effective surface temperature. Initially set to the value given in `Simulation.ini.star.T`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5772.0"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.star.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update order\n",
    "\n",
    "The update order of the stellar quantities in the standard model is set to"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['M', 'R', 'T', 'L']"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.star.updateorder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time\n",
    "\n",
    "`Simulation.t` is the current time and the integration variable. It starts at zero initially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Snapshots are written between $10^3$ years and $10^5$ years with 10 snapshots per time decade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.15576000e+10, 3.97286646e+10, 5.00154254e+10, 6.29656900e+10,\n",
       "       7.92691073e+10, 9.97938935e+10, 1.25633068e+11, 1.58162662e+11,\n",
       "       1.99114995e+11, 2.50670927e+11, 3.15576000e+11, 3.97286646e+11,\n",
       "       5.00154254e+11, 6.29656900e+11, 7.92691073e+11, 9.97938935e+11,\n",
       "       1.25633068e+12, 1.58162662e+12, 1.99114995e+12, 2.50670927e+12,\n",
       "       3.15576000e+12])"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.t.snapshots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.1"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.t.cfl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrator\n",
    "\n",
    "The integrator that is used by default as two integration instructions. One for the gas and one for the dust."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Integrator (Default integrator)"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.integrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Instruction (Dust: implicit 1st-order direct solver),\n",
       " Instruction (Gas: implicit 1st-order direct solver)]"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.integrator.instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "Except for the values at the boundaries, the advective source terms are given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.85143464e-04, -3.15967530e-11, -1.80453135e-11, -1.68008509e-11,\n",
       "       -1.56395333e-11, -1.45558169e-11, -1.35445284e-11, -1.26008399e-11,\n",
       "       -1.17202464e-11, -1.08985441e-11, -1.01318100e-11, -9.41638378e-12,\n",
       "       -8.74884973e-12, -8.12602081e-12, -7.54492330e-12, -7.00278259e-12,\n",
       "       -6.49700996e-12, -6.02519016e-12, -5.58506990e-12, -5.17454708e-12,\n",
       "       -4.79166070e-12, -4.43458153e-12, -4.10160329e-12, -3.79113452e-12,\n",
       "       -3.50169095e-12, -3.23188835e-12, -2.98043588e-12, -2.74612994e-12,\n",
       "       -2.52784828e-12, -2.32454465e-12, -2.13524373e-12, -1.95903636e-12,\n",
       "       -1.79507516e-12, -1.64257036e-12, -1.50078596e-12, -1.36903605e-12,\n",
       "       -1.24668147e-12, -1.13312657e-12, -1.02781620e-12, -9.30232957e-13,\n",
       "       -8.39894458e-13, -7.56350873e-13, -6.79182551e-13, -6.07997779e-13,\n",
       "       -5.42430665e-13, -4.82139127e-13, -4.26802977e-13, -3.76122108e-13,\n",
       "       -3.29814763e-13, -2.87615883e-13, -2.49275545e-13, -2.14557473e-13,\n",
       "       -1.83237630e-13, -1.55102892e-13, -1.29949808e-13, -1.07583453e-13,\n",
       "       -8.78163775e-14, -7.04676704e-14, -5.53621362e-14, -4.23296062e-14,\n",
       "       -3.12043925e-14, -2.18248951e-14, -1.40333743e-14, -7.67589565e-15,\n",
       "       -2.60245162e-15,  1.33274078e-15,  4.27077277e-15,  6.34742190e-15,\n",
       "        7.69239796e-15,  8.42848892e-15,  8.67064796e-15,  8.52507094e-15,\n",
       "        8.08832069e-15,  7.44655832e-15,  6.67494233e-15,  5.83725208e-15,\n",
       "        4.98578315e-15,  4.16154774e-15,  3.39479409e-15,  2.70583646e-15,\n",
       "        2.10616296e-15,  1.59976499e-15,  1.18461279e-15,  8.54187739e-16,\n",
       "        5.98977857e-16,  4.07847906e-16,  2.69210409e-16,  1.71946744e-16,\n",
       "        1.06055166e-16,  6.30311837e-17,  3.60109638e-17,  1.97267674e-17,\n",
       "        1.03327890e-17,  5.15971756e-18,  2.44844236e-18,  1.10031050e-18,\n",
       "        4.66550552e-19,  1.85917852e-19,  4.59212470e-20,  0.00000000e+00])"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.Sigma.jacobian() @ sim.gas.Sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The total source terms, i.e., including external sources and excluding the boundaries are given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ 5.85143464e-04 -3.15967530e-11 -1.80453135e-11 -1.68008509e-11\n",
       " -1.56395333e-11 -1.45558169e-11 -1.35445284e-11 -1.26008399e-11\n",
       " -1.17202464e-11 -1.08985441e-11 -1.01318100e-11 -9.41638378e-12\n",
       " -8.74884973e-12 -8.12602081e-12 -7.54492330e-12 -7.00278259e-12\n",
       " -6.49700996e-12 -6.02519016e-12 -5.58506990e-12 -5.17454708e-12\n",
       " -4.79166070e-12 -4.43458153e-12 -4.10160329e-12 -3.79113452e-12\n",
       " -3.50169095e-12 -3.23188835e-12 -2.98043588e-12 -2.74612994e-12\n",
       " -2.52784828e-12 -2.32454465e-12 -2.13524373e-12 -1.95903636e-12\n",
       " -1.79507516e-12 -1.64257036e-12 -1.50078596e-12 -1.36903605e-12\n",
       " -1.24668147e-12 -1.13312657e-12 -1.02781620e-12 -9.30232957e-13\n",
       " -8.39894458e-13 -7.56350873e-13 -6.79182551e-13 -6.07997779e-13\n",
       " -5.42430665e-13 -4.82139127e-13 -4.26802977e-13 -3.76122108e-13\n",
       " -3.29814763e-13 -2.87615883e-13 -2.49275545e-13 -2.14557473e-13\n",
       " -1.83237630e-13 -1.55102892e-13 -1.29949808e-13 -1.07583453e-13\n",
       " -8.78163775e-14 -7.04676704e-14 -5.53621362e-14 -4.23296062e-14\n",
       " -3.12043925e-14 -2.18248951e-14 -1.40333743e-14 -7.67589565e-15\n",
       " -2.60245162e-15  1.33274078e-15  4.27077277e-15  6.34742190e-15\n",
       "  7.69239796e-15  8.42848892e-15  8.67064796e-15  8.52507094e-15\n",
       "  8.08832069e-15  7.44655832e-15  6.67494233e-15  5.83725208e-15\n",
       "  4.98578315e-15  4.16154774e-15  3.39479409e-15  2.70583646e-15\n",
       "  2.10616296e-15  1.59976499e-15  1.18461279e-15  8.54187739e-16\n",
       "  5.98977857e-16  4.07847906e-16  2.69210409e-16  1.71946744e-16\n",
       "  1.06055166e-16  6.30311837e-17  3.60109638e-17  1.97267674e-17\n",
       "  1.03327890e-17  5.15971756e-18  2.44844236e-18  1.10031050e-18\n",
       "  4.66550552e-19  1.85917852e-19  4.59212470e-20  0.00000000e+00]"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.gas.Sigma.jacobian() @ sim.gas.Sigma + sim.gas.S.ext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "Except for the values at the boundaries, the hydrodynamic and coagulation source terms are given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8.46614043e-07,  8.94324218e-07,  9.44723058e-07, ...,\n",
       "         2.50288884e-30,  3.47775276e-30,  4.83232179e-30],\n",
       "       [-7.89636355e-08, -8.00467456e-08, -4.40998597e-08, ...,\n",
       "        -5.22524344e-34, -8.24485964e-34, -1.24044831e-33],\n",
       "       [-6.18577567e-08, -6.27030779e-08, -3.45488454e-08, ...,\n",
       "        -8.02854010e-34, -1.31768376e-33, -2.10795758e-33],\n",
       "       ...,\n",
       "       [-1.01704177e-58, -1.26656663e-58, -1.57727488e-58, ...,\n",
       "        -1.41404653e-47, -1.76076662e-47, -2.19250146e-47],\n",
       "       [-3.09067346e-59, -3.84856568e-59, -4.79229156e-59, ...,\n",
       "        -4.29489341e-48, -5.34798875e-48, -6.65930000e-48],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, ...,\n",
       "         0.00000000e+00,  0.00000000e+00,  0.00000000e+00]])"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(sim.dust.Sigma.jacobian() @ sim.dust.Sigma.ravel()).reshape(sim.dust.Sigma.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the `ravel()` and `reshape()` operations, since the Jacobian works on a one-dimensional state vector.\n",
    "\n",
    "The total source terms except at the boundaries are given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[ 8.46614043e-07  8.94324218e-07  9.44723058e-07 ...  2.50288884e-30\n",
       "   3.47775276e-30  4.83232179e-30]\n",
       " [-7.89636355e-08 -8.00467456e-08 -4.40998597e-08 ... -5.22524344e-34\n",
       "  -8.24485964e-34 -1.24044831e-33]\n",
       " [-6.18577567e-08 -6.27030779e-08 -3.45488454e-08 ... -8.02854010e-34\n",
       "  -1.31768376e-33 -2.10795758e-33]\n",
       " ...\n",
       " [-1.01704177e-58 -1.26656663e-58 -1.57727488e-58 ... -1.41404653e-47\n",
       "  -1.76076662e-47 -2.19250146e-47]\n",
       " [-3.09067346e-59 -3.84856568e-59 -4.79229156e-59 ... -4.29489341e-48\n",
       "  -5.34798875e-48 -6.65930000e-48]\n",
       " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00 ...  0.00000000e+00\n",
       "   0.00000000e+00  0.00000000e+00]]"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(sim.dust.Sigma.jacobian() @ sim.dust.Sigma.ravel()).reshape(sim.dust.Sigma.shape) + sim.dust.S.ext"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of a successful integration step, the floor values and boundaries are enforced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writer\n",
    "\n",
    "`DustPy` uses by default the `hdf5writer` of `simframe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Writer (HDF5 file format using h5py)\n",
       "------------------------------------\n",
       "    Data directory : data\n",
       "    File names     : data/data0000.hdf5\n",
       "    Overwrite      : False\n",
       "    Dumping        : True\n",
       "    Options        : {'com': 'lzf', 'comopts': None}\n",
       "    Verbosity      : 1"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim.writer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is by default writing dump files and prevents overwriting of already existing data files."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}