Dust evolution in protoplanetary disks
Standard functions for dust physics including particle dynamics, coagulation kernels, velocities, probabilities, fluxes, and source terms. These functions implement the core physics of dust evolution in protoplanetary disks with performance-critical calculations accelerated by Fortran extensions.
Functions that calculate fundamental dust properties from the simulation state.
from dustpy.std import dust
def dust.a(sim):
"""
Calculate particle size from density and filling factor.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Particle sizes [cm]
"""
def dust.rho_midplane(sim):
"""
Calculate midplane mass density of dust.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Midplane mass density [g/cm³]
"""
def dust.eps(sim):
"""
Calculate vertically integrated dust-to-gas ratio.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Dust-to-gas ratio
"""
def dust.H(sim):
"""
Calculate dust scale height using Dubrulle et al. (1995) prescription.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Dust scale heights [cm]
"""
def dust.St_Epstein_StokesI(sim):
"""
Calculate Stokes number in Epstein and Stokes I regimes.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Stokes numbers
"""Functions for calculating dust diffusivity and transport coefficients.
def dust.D(sim):
"""
Calculate dust diffusivity from turbulent mixing.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Dust diffusivity [cm²/s]
"""Functions for calculating particle collision rates and coagulation kernels.
def dust.kernel(sim):
"""
Calculate vertically integrated collision kernel.
The collision kernel determines the rate at which particles
of different masses collide and potentially stick or fragment.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Collision kernel [cm²/s]
"""Functions for calculating various components of dust particle velocities.
def dust.vrad(sim):
"""
Calculate radial dust velocity including drift and diffusion.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Radial velocity [cm/s]
"""
def dust.vdriftmax(sim):
"""
Calculate maximum drift velocity including back reaction effects.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Maximum drift velocity [cm/s]
"""
def dust.vrel_tot(sim):
"""
Calculate total relative velocity between particles (RMS).
Combines contributions from various relative velocity sources.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Total relative velocity [cm/s]
"""
def dust.vrel_azimuthal_drift(sim):
"""
Calculate relative velocity from azimuthal drift differences.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Relative velocity from azimuthal drift [cm/s]
"""
def dust.vrel_brownian_motion(sim):
"""
Calculate relative velocity from Brownian motion.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Relative velocity from Brownian motion [cm/s]
"""
def dust.vrel_radial_drift(sim):
"""
Calculate relative velocity from radial drift differences.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Relative velocity from radial drift [cm/s]
"""
def dust.vrel_turbulent_motion(sim):
"""
Calculate relative velocity from turbulent motion.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Relative velocity from turbulence [cm/s]
"""
def dust.vrel_vertical_settling(sim):
"""
Calculate relative velocity from vertical settling differences.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Relative velocity from vertical settling [cm/s]
"""Functions for calculating sticking and fragmentation probabilities.
def dust.p_stick(sim):
"""
Calculate sticking probability for particle collisions.
Determines the likelihood that colliding particles will
stick together rather than bounce or fragment.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Sticking probability (0-1)
"""
def dust.p_frag(sim):
"""
Calculate fragmentation probability for particle collisions.
Determines the likelihood that colliding particles will
fragment into smaller pieces.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Fragmentation probability (0-1)
"""Functions for calculating mass transport through advection and diffusion.
def dust.F_adv(sim, Sigma=None):
"""
Calculate advective mass flux.
Parameters:
- sim: Simulation object
- Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
Returns:
numpy.ndarray: Advective flux [g/cm/s]
"""
def dust.F_diff(sim, Sigma=None):
"""
Calculate diffusive mass flux.
Parameters:
- sim: Simulation object
- Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
Returns:
numpy.ndarray: Diffusive flux [g/cm/s]
"""
def dust.F_tot(sim, Sigma=None):
"""
Calculate total mass flux (advective + diffusive).
Parameters:
- sim: Simulation object
- Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
Returns:
numpy.ndarray: Total mass flux [g/cm/s]
"""Functions for calculating mass source and sink terms from various physical processes.
def dust.S_coag(sim, Sigma=None):
"""
Calculate coagulation source terms.
Computes mass exchange between size bins due to particle
collisions, sticking, and fragmentation processes.
Parameters:
- sim: Simulation object
- Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
Returns:
numpy.ndarray: Coagulation source terms [g/cm²/s]
"""
def dust.S_hyd(sim, Sigma=None):
"""
Calculate hydrodynamic source terms.
Computes source terms from hydrodynamic processes including
radial transport and diffusion.
Parameters:
- sim: Simulation object
- Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
Returns:
numpy.ndarray: Hydrodynamic source terms [g/cm²/s]
"""
def dust.S_tot(sim, Sigma=None):
"""
Calculate total source terms (coagulation + hydrodynamic).
Parameters:
- sim: Simulation object
- Sigma: Surface density array (optional, uses sim.dust.Sigma if None)
Returns:
numpy.ndarray: Total source terms [g/cm²/s]
"""Functions for setting up initial dust distributions and parameters.
def dust.MRN_distribution(sim):
"""
Set up MRN (Mathis, Rumpl, Nordsieck) initial particle mass distribution.
Creates a power-law size distribution commonly used as initial
conditions for dust evolution simulations.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Initial surface density distribution [g/cm²]
"""
def dust.SigmaFloor(sim):
"""
Calculate floor value for dust surface density.
Prevents numerical issues by enforcing minimum dust densities.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Floor values [g/cm²]
"""Functions for numerical integration of dust evolution equations.
def dust.jacobian(sim, x, dx=None):
"""
Calculate Jacobian matrix for implicit integration schemes.
Parameters:
- sim: Simulation object
- x: Current state vector
- dx: State perturbation (optional)
Returns:
numpy.ndarray: Jacobian matrix
"""
def dust.Sigma_deriv(sim, t, Sigma):
"""
Calculate time derivative of dust surface density.
Core function for explicit integration schemes.
Parameters:
- sim: Simulation object
- t: Current time [s]
- Sigma: Current surface density [g/cm²]
Returns:
numpy.ndarray: Time derivative [g/cm²/s]
"""
class dust.impl_1_direct:
"""
Direct implicit integration scheme for dust evolution.
Provides first-order implicit integration with direct
matrix solving for stable evolution of stiff systems.
"""
def __init__(self):
"""Initialize implicit integration scheme."""
...Helper functions for dust simulation setup and maintenance.
def dust.boundary(sim):
"""
Set boundary conditions for dust quantities.
Applies inner and outer boundary conditions based on
the configured boundary condition types and values.
Parameters:
- sim: Simulation object
"""
def dust.enforce_floor_value(sim):
"""
Enforce minimum floor values for dust quantities.
Prevents numerical issues by ensuring all dust quantities
remain above specified minimum values.
Parameters:
- sim: Simulation object
"""
def dust.coagulation_parameters(sim):
"""
Calculate and update coagulation-related parameters.
Updates collision kernels, probabilities, and other
coagulation parameters based on current simulation state.
Parameters:
- sim: Simulation object
"""
def dust.dt(sim):
"""
Calculate appropriate time step for dust evolution.
Determines stable time step based on current dust
evolution rates and numerical stability criteria.
Parameters:
- sim: Simulation object
Returns:
float: Time step [s]
"""
def dust.dt_adaptive(sim):
"""
Calculate adaptive time step for dust evolution.
Uses error estimation to adjust time step dynamically
for optimal balance of accuracy and efficiency.
Parameters:
- sim: Simulation object
Returns:
float: Adaptive time step [s]
"""
def dust.prepare(sim):
"""
Prepare implicit dust integration step.
Stores the current value of surface density in a hidden
field for use in implicit integration schemes.
Parameters:
- sim: Simulation object
"""
def dust.finalize_explicit(sim):
"""
Finalize explicit integration step.
Applies boundary conditions and enforces floor values
after an explicit integration step.
Parameters:
- sim: Simulation object
"""
def dust.finalize_implicit(sim):
"""
Finalize implicit integration step.
Applies boundary conditions, enforces floor values, and
updates dependent quantities after implicit integration.
Parameters:
- sim: Simulation object
"""
def dust.set_implicit_boundaries(sim):
"""
Calculate fluxes at boundaries after implicit integration.
Sets boundary source terms based on the change in surface
density during the implicit integration step.
Parameters:
- sim: Simulation object
"""from dustpy import Simulation
from dustpy.std import dust
# Create and initialize simulation
sim = Simulation()
sim.makegrids()
sim.initialize()
# Calculate dust properties
particle_sizes = dust.a(sim) # [cm]
stokes_numbers = dust.St_Epstein_StokesI(sim)
dust_scale_heights = dust.H(sim) # [cm]
collision_kernel = dust.kernel(sim) # [cm²/s]# Calculate various velocity components
radial_velocity = dust.vrad(sim) # [cm/s]
total_relative_v = dust.vrel_tot(sim) # [cm/s]
brownian_v = dust.vrel_brownian_motion(sim) # [cm/s]
drift_v = dust.vrel_radial_drift(sim) # [cm/s]
# Maximum drift velocity with back reaction
max_drift = dust.vdriftmax(sim) # [cm/s]# Calculate source terms
coag_sources = dust.S_coag(sim) # [g/cm²/s]
hydro_sources = dust.S_hyd(sim) # [g/cm²/s]
total_sources = dust.S_tot(sim) # [g/cm²/s]
# Analyze mass fluxes
advective_flux = dust.F_adv(sim) # [g/cm/s]
diffusive_flux = dust.F_diff(sim) # [g/cm/s]# Calculate appropriate time steps
dt_dust = dust.dt(sim) # [s]
dt_adaptive = dust.dt_adaptive(sim) # [s]
print(f"Dust time step: {dt_dust/dust.year:.2e} years")
print(f"Adaptive time step: {dt_adaptive/dust.year:.2e} years")Install with Tessl CLI
npx tessl i tessl/pypi-dustpy