4. Radiative environment

So far we have set up the grid and star objects, and one or more fluid objects. We can now set up a radiative environment that combines all of these ingredients to perform radiative transfer calculations. This is done with the fld.FLD.RadiativeEnvironment class, which is initialized with the following call signature:

fld.FLD.RadiativeEnvironment(
    fluids,
    verbose=True,
    constant_fluxlimiter=False,
    preconditioner='ilu',
    implicit_viscosity=False,
    debug=False
)

The only mandatory argument (and likely the only one you need to worry about) is the list of fluid objects that you want to include in the radiative transfer calculations. For example, with what we have set up so far, we can initialize the radiative environment with:

Disk = fld.FLD.RadiativeEnvironment([gas])

This is enough to declare a protoplanetary disk that accounts for viscous heating with α viscosity, ray-traced irradiation with gray opacities, and radiative transfer with Rosseland and Planck mean opacities.

The RadiativeEnvironment class is slightly different as most of the physics have already been set up in the fluid objects. However, we also need to prescribe boundary conditions for the radiative transfer calculations. This is done with the method:

Disk.declare_boundaries(bounds, vals)

The method expects two lists, 6 elements each, that specify the type of boundary condition and the corresponding value for each of the 6 boundaries (inner and outer in each of the 3 dimensions). The allowed boundary types are fld.grid.BOUNDARY_FIXEDVALUE for fixed-value boundaries, and fld.grid.BOUNDARY_REFLECTIVE for zero-gradient (reflective) boundaries. The expected values are the radiation energy density \(E_\mathrm{r}\) for fixed-value boundaries, and are not needed for reflective boundaries. For example, to set a fixed value boundary condition that corresponds to a radiation temperature of 10 K at the inner radial edge, and reflective boundaries everywhere else, we can do:

from astropy import units as u

bounds = [
    fld.grid.BOUNDARY_FIXEDVALUE,  # inner radial
    fld.grid.BOUNDARY_REFLECTIVE,  # outer radial
    fld.grid.BOUNDARY_REFLECTIVE,  # inner polar
    fld.grid.BOUNDARY_REFLECTIVE,  # outer polar
    fld.grid.BOUNDARY_REFLECTIVE,  # inner azimuthal
    fld.grid.BOUNDARY_REFLECTIVE,  # outer azimuthal
]

vals = [
    fld.tools.T_to_Er(10 * u.K), # Erad at inner radial edge for T=10 K
    0, 0, 0, 0, 0  # not needed for reflective boundaries
]

Disk.declare_boundaries(bounds, vals)

We are now ready to perform radiative transfer calculations with the Disk object. A typical call to the diffusion solver looks like:

Disk.advance(dt=<Quantity 1. yr>)

which advances the radiation field by a time step of 1 year.

The RadiativeEnvironment class also includes a method that calls compute_hydrostatic_equilibrium for all fluids and then applies the resulting density profiles, taking into account which fluids correspond to gas or dust. It is simply called with:

Disk.enforce_hydrostatic_equilibrium()

Additional control over the solver can be achieved with the Disk.configure_solver method, which accepts keyword arguments that are passed to the underlying solver:

  • solver: the name of the sparse matrix solver used. Options are bicgstab (default), spsolve, gmres, cg, cgs, lgmres, minres, and qmr, but note that the module has only really been tested with bicgstab and sometimes spsolve.

  • rtol: the relative tolerance for the iterative solver (default is \(10^{-8}\)).

  • maxiter: the maximum number of iterations for the iterative solver (default is 10000).

  • guess: if True, passes the current radiation energy density as an initial guess, which can be useful for speeding up convergence when performing multiple calls to advance in a time-dependent simulation. Otherwise, the solver starts with a zero-valued initial guess (default is True).

Regarding the preconditioner argument in the class initializer: the default ilu option applies an incomplete LU factorization to the sparse matrix, which can significantly speed up convergence for large problems. However, for very large problems, the factorization itself can be computationally expensive and may even fail due to memory constraints. In such cases, it may be better to set the preconditioner to jacobi (i.e., diagonal scaling).

In addition to the above preconditioners, the class supports “rescaled” variants of the default preconditioners, where the sparse matrix is rescaled by the inverse of its diagonal before applying the preconditioner. This is a typically cheap operation that can speed up convergence for some problems. To use the rescaled variants, simply set preconditioner to rescaled ilu or rescaled jacobi.

We can now follow this workflow to set up more meaningful problems and run radiative transfer calculations on them. We will go through some examples in the examples section.