Advanced example: protoplanetary disk

This example shows how to set up a protoplanetary disk with viscous heating, frequency-dependent ray-traced irradiation, and radiative transfer with Rosseland and Planck mean opacities. Several similar examples are included in the examples folder of the repository, which you can run and modify as needed.

To simplify the computation of frequency-dependent opacities, we will use the optool package. After running:

optool -radmc 01um -a 0.1

we have a file dustkappa_01um.inp containing absorption and scattering opacities and the scattering asymmetry parameter as a function of wavelength for 0.1 μm grains. We will use this to set up frequency-dependent irradiation in our example.

from astropy import units as u
import pyFLD as fld
import numpy as np

wl, kabs, ksca, g_asy = fld.tools.parse_optool_file('dustkappa_01um.inp')

# star object
star = fld.star.Star(M=1*u.M_sun, L=1*u.L_sun, T=5772*u.K, wl=wl)

# grid object
ri   = fld.grid.get_domain_with_ghosts(0.1, 10, Nx=128, log=True) * u.au
thi  = fld.grid.get_domain_with_ghosts(np.pi/2-0.35, np.pi/2, Nx=256) # half disk in polar direction
phii = fld.grid.get_domain_with_ghosts(0, 2*np.pi, Nx=1) # axisymmetric grid
grid = fld.grid.Grid(ri, thi, phii, geometry='spherical')

# fluid objects
Gas = fld.fluids.Gas(grid=grid, star=star, irradiation='gray', viscosity=True)
Dust = fld.fluids.Dust(grid=grid, star=star, irradiation='freqdep',
                    kabs_wl=kabs, ksca_wl=ksca, g_wl=g_asy,
                grain_rho=1*u.g/u.cm**3, grain_size=0.1*u.um)

# set up fluids: density has a radial power law and Gaussian vertical profile, temperature is a radial power law, and alpha viscosity is constant
R_ = grid.R.to_value('au')
z_ = grid.z.to_value('au')
H_ = 0.03 * R_ ** 1.25
rho = 3e-10 * u.g/u.cm**3 * R_ ** -1.5 * np.exp(-0.5*(z_/H_)**2)
T = 150 * u.K * R_ ** -0.5

# gas properties
Gas.set_density(rho * np.ones(grid.shape))
Gas.set_temperature(T * np.ones(grid.shape))
Gas.enroll_alpha(lambda f: 1e-4) # w/ viscous heating
Gas.enroll_tau0(fld.tools.get_tau0_Flock_2019)
Gas.enroll_kappaR(lambda f: 1.0 * u.cm**2/u.g)
Gas.enroll_kappaP(lambda f: 1.0 * u.cm**2/u.g)
Gas.enroll_kappaPstar(lambda f: 1.0 * u.cm**2/u.g)

# dust properties
Dust.set_density(rho * 1e-3 * np.ones(grid.shape))
Dust.set_temperature(T * np.ones(grid.shape))
Dust.enroll_alpha(lambda f: 1e-2) # w/o viscous heating
Dust.enroll_col0(fld.tools.get_col0_Flock_2019)
# no need to enroll kappaR/kappaP/kappaPstar for the dust as they will be automatically computed from the provided frequency-dependent opacity model

Gas.setup()
Dust.setup()

Disk = fld.FLD.RadiativeEnvironment([Gas, Dust], name='disk')

# 10 K in all boundaries
bounds = [fld.grid.BOUNDARY_FIXEDVALUE] * 6
vals   = [fld.tools.T_to_Er(10*u.K)] * 6
bounds[3] = fld.grid.BOUNDARY_REFLECTIVE # reflective midplane
Disk.declare_boundaries(bounds, vals)

# now run model for a very large time step to get close to radiative equilibrium. This will take several steps.

for _ in range(20): Disk.advance(dt=1e10*u.yr)

# you can now plot the results, for example the midplane temperature profile:

import matplotlib.pyplot as plt

r = grid.x1[grid.act1].to_value('au')
T_mid = Disk.tmp[grid.active][0,-1,:].to_value('K') # [phi=0, theta=-1 (midplane), r=all]
plt.plot(r, T_mid, label='pyFLD')
plt.plot(r, 145 * r ** -0.49, label='Power law fit')
plt.xlabel('Radius (au)')
plt.ylabel('Midplane Temperature (K)')
plt.title('Midplane Temperature Profile')
plt.loglog()
plt.legend()
plt.show()
plt.close()


# you can also opt to enforce hydrostatic equilibrium every few steps for a more self-consistent solution
Disk.configure_solver(rtol=1e-4) # for this example
Disk.verbose = False
for i in range(5):
    for _ in range(20): Disk.advance(dt=1e10*u.yr)
    Disk.enforce_hydrostatic_equilibrium()
    print(f'finished supercycle {i+1}')

# and then plot the results again
r = grid.x1[grid.act1].to_value('au')
T_mid = Disk.tmp[grid.active][0,-1,:].to_value('K') # [phi=0, theta=-1 (midplane), r=all]
plt.plot(r, T_mid, label='pyFLD')
plt.plot(r, 145 * r ** -0.49, label='Power law fit')
plt.xlabel('Radius (au)')
plt.ylabel('Midplane Temperature (K)')
plt.title('Midplane Temperature Profile')
plt.loglog()
plt.legend()
plt.show()