Best Practices

Recommendations for getting the most out of pyFLD.

Choosing the right grid resolution

Rule of thumb: Start coarse, refine until results converge.

  • Radial resolution: 64–256 cells for protoplanetary disks. More for steep features (e.g., snow lines, opacity gaps).

  • Vertical resolution: 32–128 cells. More if you care about disk structure near the midplane or the photosphere.

  • Azimuthal resolution: 1 cell for axisymmetric; 64–512 for typical non-axisymmetric structures.

Always use logarithmic spacing in radius (and when appropriate, in other directions):

ri = fld.grid.get_domain_with_ghosts(0.1, 1000, Nx=128, log=True) * u.au

Gray vs. frequency-dependent opacities

Use gray opacities when:

  • You want quick estimates

  • Dust composition is uniform or mean properties suffice

  • Computational cost is critical

  • Wavelength-dependent properties don’t strongly affect temperature

Use frequency-dependent opacities when:

  • You want accurate temperature structure, especially for dust

  • Dust opacity has strong features (e.g., silicate, ice, graphite absorption)

  • Irradiation wavelength dependence matters for your problem

Wavelength grid for frequency-dependent opacities

Resolution matters:

  • Too few bins (~10): Mean opacities poorly estimated; misses absorption features

  • Sweet spot (~50–100): Good accuracy and reasonable speed

  • Too many (~1000): Slow, diminishing returns

Solver configuration

Default settings (usually fine):

Disk = fld.FLD.RadiativeEnvironment([gas], preconditioner='ilu')
Disk.configure_solver(solver='bicgstab', rtol=1e-8, maxiter=10000)

For faster convergence (at cost of potentially lower accuracy):

Disk.configure_solver(rtol=1e-4)

For very large problems (\(> 10^8\) cells):

Disk = fld.FLD.RadiativeEnvironment([gas], preconditioner='rescaled jacobi') # or 'rescaled ilu' if memory allows
Disk.configure_solver(rtol=1e-6)

To speed up multiple ``advance()`` calls (e.g., time evolution):

Disk.configure_solver(guess=True)  # Use previous solution as initial guess
for _ in range(100):
    Disk.advance(dt=1e5*u.yr)

Boundary conditions

For disk models (spherical geometry):

# Typical: fixed T at all boundaries (radiative equilibrium)
T_bc = 10 * u.K  # Outer boundary temperature
bounds = [fld.grid.BOUNDARY_FIXEDVALUE] * 6
vals = [fld.tools.T_to_Er(T_bc)] * 6

# Keep midplane reflective if computing its vertical structure from irradiation
bounds[3] = fld.grid.BOUNDARY_REFLECTIVE # the value is ignored when reflective

Disk.declare_boundaries(bounds, vals)

Time stepping

Choosing time step dt:

FLD is implicit, so large time steps are stable (unlike explicit methods). However:

  • Too small: Wastes computation time

  • Too large: May not resolve fast transients; convergence slower

Start with:

dt = (L**2 / D)  # Characteristic diffusion time
# where L is domain size, D = lambda*c/(kappa*rho) is diffusion coefficient

For a disk at 1 AU with typical parameters, try dt\(\sim 10^6-10^7\,\mathrm{yr}\) and adjust.

Opacity functions

Constant opacity (simplest for testing):

Gas.enroll_kappaR(lambda f: 0.1 * u.cm**2/u.g)
Gas.enroll_kappaP(lambda f: 0.1 * u.cm**2/u.g)

Temperature-dependent opacity (more realistic):

def opacity_P(f):
    T = f.tmp.to_value('K')
    # Simplified power law: kappa propto T^2 (example)
    return 1e-4 * u.cm**2/u.g * (T / 100)**2

Gas.enroll_kappaP(opacity_P)

Automatic from frequency-dependent data (recommended):

When using frequency-dependent opacities with build_mean_opacities=True (default), Rosseland and Planck means are computed automatically. No need to enroll them.

Viscosity

When to include:

  • Standard α-viscosity in protoplanetary disks

  • Studying viscous heating effects

  • Modeling accretion

Simple constant α:

Gas.enroll_alpha(lambda f: 1e-4)  # Shakura-Sunyaev alpha

variable α (e.g., α(T) from Cecil & Flock 2024):

Gas.enroll_alpha(fld.tools.get_alpha_Cecil_Flock_2024)

Dust without viscous heating (but with turbulent diffusion):

Dust = fld.fluids.Dust(..., viscosity=False)
Dust.enroll_alpha(lambda f: 1e-3)  # For settling, not heating

Irradiation

Initial setup:

# Gray irradiation: must enroll tau0 function
Gas.enroll_tau0(fld.tools.get_tau0_Flock_2019)

# Frequency-dependent: must enroll col0 function
Dust.enroll_col0(fld.tools.get_col0_Flock_2019)

Check irradiation is enabled:

print(Gas.irradiation)  # Should print 'gray' or 'freqdep', not False

Avoid common mistakes:

  • Do not forget to enroll tau0 (gray) or col0 (freqdep)

  • Ensure wavelength arrays match between Star and Dust

  • For frequency-dependent, provide absorption opacity (scattering can be zeros if not needed)

Hydrostatic equilibrium

When to use:

For self-consistent disk structure (gas pressure support, dust settling), especially with multiple species.

Best practice workflow:

# Initial setup
Disk = fld.FLD.RadiativeEnvironment([Gas, Dust])
Disk.declare_boundaries(bounds, vals)

# Solve to target temperature profile (here, to equilibrium)
for _ in range(50):
    Disk.advance(dt=1e8*u.yr)

# Update density to satisfy hydrostatic equilibrium
Disk.enforce_hydrostatic_equilibrium()

# Continue with refined calculation (T will adjust)
for _ in range(100):
    Disk.advance(dt=1e8*u.yr)

# Update HSE periodically if needed
for cycle in range(5):
    for _ in range(20):
        Disk.advance(dt=1e8*u.yr)
    Disk.enforce_hydrostatic_equilibrium()

Validation & testing

Analytic benchmarks:

Test against known solutions before tackling complex models:

Example: 1D diffusion — Linear temperature profile (see example-basic)

Convergence studies:

Run with different resolutions and check results converge:

for Nx in [32, 64, 128, 256]:
    grid = setup_grid(Nx=Nx)
    Disk = setup_model(grid)
    Disk.advance(dt=1e8*u.year)
    T_mid = Disk.tmp[...]
    print(f"Nx={Nx}: T_mid(r=10au) = {T_mid[...]} K")

Performance monitoring

Check solver performance:

Disk.verbose = True  # Print solver iterations
for i in range(10):
    Disk.advance(dt=1e6*u.yr)
    print(f"Step {i}: iters={Disk.iters}")

Profile the code using the debug flag to identify bottlenecks:

Gas = fld.fluids.Gas(..., debug=True)
Disk = fld.FLD.RadiativeEnvironment([Gas], debug=True)
Disk.advance(dt=1e6*u.yr) # will print progress info