Dust evolution in protoplanetary disks
Standard functions for gas disk physics including thermodynamics, viscosity, pressure profiles, and hydrodynamic evolution. These functions implement the gas disk physics with performance-critical calculations accelerated by Fortran extensions.
Functions for calculating gas temperature, pressure, and sound speed.
from dustpy.std import gas
def gas.cs_isothermal(sim):
"""
Calculate isothermal sound speed.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Isothermal sound speed [cm/s]
"""
def gas.P_midplane(sim):
"""
Calculate midplane gas pressure.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Midplane pressure [g/cm/s²]
"""
def gas.T_passive(sim):
"""
Calculate passive irradiation temperature.
Computes gas temperature from stellar irradiation
assuming passive heating without viscous dissipation.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Temperature [K]
"""
def gas.rho_midplane(sim):
"""
Calculate midplane gas mass density.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Midplane mass density [g/cm³]
"""
def gas.n_midplane(sim):
"""
Calculate midplane number density.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Midplane number density [1/cm³]
"""Functions for calculating characteristic length scales in the gas disk.
def gas.Hp(sim):
"""
Calculate pressure scale height of the gas disk.
The pressure scale height determines the vertical
extent of the gas disk and affects dust settling.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Pressure scale height [cm]
"""
def gas.eta_midplane(sim):
"""
Calculate midplane pressure gradient parameter.
The eta parameter quantifies the radial pressure gradient
and determines the strength of gas drag on dust particles.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Pressure gradient parameter
"""Functions for calculating viscosity and transport coefficients.
def gas.nu(sim):
"""
Calculate kinematic viscosity of the gas.
Uses the alpha-disk model to compute viscosity from
the turbulent alpha parameter and local disk properties.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Kinematic viscosity [cm²/s]
"""
def gas.mfp_midplane(sim):
"""
Calculate mean free path at the midplane.
The mean free path determines the gas drag regime
(Epstein vs. Stokes) for dust particles.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Mean free path [cm]
"""Functions for calculating gas velocity components from disk evolution.
def gas.vrad(sim):
"""
Calculate radial gas velocity from viscous evolution.
Includes contributions from viscous spreading and
any external torques or mass sources.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Radial velocity [cm/s]
"""
def gas.vvisc(sim):
"""
Calculate viscous radial velocity component.
The radial velocity driven purely by viscous torques
in the standard alpha-disk model.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Viscous radial velocity [cm/s]
"""
def gas.vtorque(sim):
"""
Calculate velocity contribution from external torques.
Additional radial velocity from gravitational torques
due to planets or other perturbing bodies.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Torque-driven velocity [cm/s]
"""Functions for calculating mass transport and evolution source terms.
def gas.Fi(sim):
"""
Calculate mass flux through radial interfaces.
The mass flux determines how gas surface density
evolves due to radial transport processes.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Mass flux through interfaces [g/cm/s]
"""
def gas.S_hyd(sim):
"""
Calculate hydrodynamic source terms.
Source terms from viscous evolution, including
effects of radial transport and viscous heating.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Hydrodynamic source terms [g/cm²/s]
"""
def gas.S_tot(sim):
"""
Calculate total source terms for gas evolution.
Combines hydrodynamic and external source terms
for the complete gas evolution equation.
Parameters:
- sim: Simulation object
Returns:
numpy.ndarray: Total source terms [g/cm²/s]
"""Functions for setting up initial gas surface density profiles.
def lyndenbellpringle1974(r, rc, p, Mdisk):
"""
Calculate self-similar disk profile from Lynden-Bell & Pringle (1974).
Creates an exponentially tapered power-law surface density
profile commonly used as initial conditions for disk evolution.
Parameters:
- r: Radial grid [cm]
- rc: Characteristic radius [cm]
- p: Power-law index
- Mdisk: Total disk mass [g]
Returns:
numpy.ndarray: Surface density profile [g/cm²]
"""Functions for numerical integration of gas evolution equations.
def gas.jacobian(sim, x, *args, **kwargs):
"""
Calculate Jacobian matrix for gas evolution.
Used in implicit integration schemes for stable
evolution of the gas disk equations.
Parameters:
- sim: Simulation object
- x: Current state vector
Returns:
numpy.ndarray: Jacobian matrix
"""
class gas.impl_1_direct:
"""
Direct implicit integration scheme for gas evolution.
Provides first-order implicit integration for stable
evolution of viscous disk equations.
"""
def __init__(self):
"""Initialize implicit integration scheme."""
...Helper functions for gas simulation setup and maintenance.
def gas.boundary(sim):
"""
Set boundary conditions for gas quantities.
Applies inner and outer boundary conditions for
gas surface density and other gas properties.
Parameters:
- sim: Simulation object
"""
def gas.enforce_floor_value(sim):
"""
Enforce minimum floor values for gas quantities.
Prevents numerical issues by ensuring gas properties
remain above physically reasonable minimum values.
Parameters:
- sim: Simulation object
"""
def gas.dt(sim):
"""
Calculate appropriate time step for gas evolution.
Determines stable time step based on viscous time scale
and numerical stability requirements.
Parameters:
- sim: Simulation object
Returns:
float: Time step [s]
"""
def gas.prepare(sim):
"""
Prepare gas integration step.
Stores the current value of surface density in a hidden
field for use in integration schemes.
Parameters:
- sim: Simulation object
"""
def gas.finalize(sim):
"""
Finalize gas integration step.
Applies boundary conditions, enforces floor values, and
updates dependent quantities after integration.
Parameters:
- sim: Simulation object
"""
def gas.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 gas
import dustpy.constants as c
# Create and initialize simulation
sim = Simulation()
sim.makegrids()
sim.initialize()
# Calculate basic gas properties
sound_speed = gas.cs_isothermal(sim) # [cm/s]
pressure = gas.P_midplane(sim) # [g/cm/s²]
scale_height = gas.Hp(sim) # [cm]
density = gas.rho_midplane(sim) # [g/cm³]
# Convert to useful units
print(f"Sound speed: {sound_speed/1e5:.1f} km/s")
print(f"Scale height: {scale_height/c.au:.2f} AU")# Calculate transport properties
viscosity = gas.nu(sim) # [cm²/s]
radial_velocity = gas.vrad(sim) # [cm/s]
viscous_velocity = gas.vvisc(sim) # [cm/s]
# Calculate mass flux
mass_flux = gas.Fi(sim) # [g/cm/s]
# Time scales
viscous_time = sim.grid.r**2 / viscosity # [s]
print(f"Viscous time at 1 AU: {viscous_time[50]/c.year:.1e} years")# Set up Lynden-Bell & Pringle profile
r = sim.grid.r # [cm]
rc = 10 * c.au # Characteristic radius
p = -1.5 # Power-law index
Mdisk = 0.01 * c.M_sun # Disk mass
# Calculate initial profile
initial_sigma = lyndenbellpringle1974(r, rc, p, Mdisk)
sim.gas.Sigma[:] = initial_sigma # Set as initial condition# Calculate pressure gradient parameter
eta = gas.eta_midplane(sim)
# This affects dust drift velocities
print(f"Pressure gradient parameter at 1 AU: {eta[50]:.3f}")
# Mean free path affects drag regime
mfp = gas.mfp_midplane(sim) # [cm]
particle_size = 0.01 # 100 micron particles
print(f"Knudsen number: {particle_size/mfp[50]:.2e}")# Calculate passive disk temperature
temperature = gas.T_passive(sim) # [K]
# Set up temperature-dependent properties
sim.gas.T[:] = temperature
# Recalculate sound speed with new temperature
updated_cs = gas.cs_isothermal(sim) # [cm/s]Install with Tessl CLI
npx tessl i tessl/pypi-dustpy