Version 9
advanced radiative transfer in dusty systems
What's new in SKIRT 9

Introduction

During the last decade, SKIRT has been actively applied to various science cases and the code has grown into an advanced, state-of-the-art 3D dust continuum radiative transfer code. The recent version, SKIRT 8, has proven to be very robust and has been stable for over a year with just a few minor updates. While this is very satisfying, the core SKIRT team has set new long-term goals for the project.

Over time, we envision SKIRT to become a code that, while maintaining its current capabilities, also:

These ambitious goals have some fundamental consequences for the structure of the SKIRT code. Foremost, hydrogen line transfer requires support for kinematics because absorption and scattering cross sections may vary significantly over the wavelength Doppler-shifts caused by the relative motion of sources and media. This in turn requires tracking precise, variable wavelengths for photon packets moving through the medium, as opposed to using a fixed global wavelength grid as in the SKIRT 8 approach.

A second important challenge is to handle the growing complexity of the input model. Supporting kinematics requires specifying velocities across the spatial domain. Similarly, in addition to changes in the photon cycle, supporting spheroidal dust grains requires specifying the direction and degree of grain alignment across the spatial domain. The SKIRT 8 geometry and decorator system has been designed to model a density field normalized to unit mass. Generalizing these classes to also model vector fields (for velocities and directions) and unnormalized scalar fields (for alignment degrees) seems to be somewhere between daunting and impossible.

During the last year, we developed a new major version, SKIRT 9, restructuring the code to lay the foundation for the ambitious future outlined above. Because supporting the long-term goals required some far-reaching changes anyway, we took the opportunity to streamline or otherwise improve many other areas of the code along the way. The following section summarizes the most important changes from SKIRT 8 to SKIRT 9, and the remainder of this page provides more detail on some of these issues. The discussion is targeted towards readers who are familiar with SKIRT 8 at least as a user. Topics are covered in somewhat arbitrary order.

See also:

Summary of changes

SKIRT 9 includes all SKIRT 8 features (possibly with some adjustements) except those that were dropped on purpose for one of two reasons:

  1. The feature is very hard to implement in the new structure or dropping it allowed significantly streamlining the code. In most cases, there is an acceptable workaround for accomplishing the same goal.
  2. The feature seemed to be no longer used, but it could be added later if needed.

Major incompatibilities:

Other changes:

Major dropped features:

Minor dropped features (because no longer used):

Major new features:

Minor new features:

Performance

Photon packet properties

The fundamental properties of a photon packet include its wavelength and its weight. The wavelength property specifies the wavelength (or equivalently, the frequency) of all photons in the packet. The weight property specifies the number of photons carried by the packet, or more precisely the number of photons per unit of time.

At launch (from a primary or secondary source), a photon packet receives a wavelength (sampled from the source spectrum) and a luminosity (its share of the total luminosity of the source). The wavelength is stored as given. The luminosity is converted to a weight (number of photons) for storage in the photon packet. During a photon packet’s life cycle, updates can occur to both its weight (e.g. because of biasing) and its wavelength (e.g. after scattering by a moving medium). Because these updates can be fractional, both wavelength and weight are stored as floating-point values.

Adjusting a photon packet's wavelength indirectly affects the luminosity represented by the packet, because the latter is directly proportional to the frequency and thus inversely proportional to the wavelength.

Kinematics

The wavelength of a photon packet is defined relative to the model coordinate system (i.e. a medium at rest relative to the model coordinate system sees this wavelength). Velocities of sources and media are also defined relative to the model coordinate system. The bulk velocity of the aggregate medium in a spatial grid cell is determined by averaging the bulk velocities of all media components weighed by density. Instruments are considered to be at rest relative to the model coordinate system.

When a photon packet is launched, its wavelength is Doppler shifted according to the component of the source velocity in the photon packet’s direction. When a photon packet interacts with a medium, the perceived wavelength is derived by Doppler shifting the packet’s wavelength according to the component of the aggregate medium bulk velocity in the photon packet’s incoming direction. Specifically, registration of a photon packet’s contribution to the radiation field uses this perceived wavelength.

After a scattering interaction, the photon packet’s wavelength is replaced by the perceived wavelength, Doppler shifted according to the component of the medium velocity in the photon packet’s outgoing direction.

Distributing photon packets among sources

The source system distributes the emitted photon packets among the primary sources in the simulation based on three types of information: the luminosity weight of each source, a user-configured weight for each source (defaulted to equal weight for all sources), and a user-configured bias factor (defaulted to 0.5) that shifts between distributing according to the above weights and simply allocating an equal number of packets to each source.

In turn, each source determines how to further distribute the photon packets allocated to it. For example, a geometric source distributes photon packets according to its (built-in) spatial density distribution. An imported source has a scheme similar to the source system, using a user-configured composite bias factor to shift between a distribution according to the luminosity weights of each sub-source (particle or cell) and a linear equal-weight allocation.

The secondary source system also uses the latter scheme to distribute photon packets among the cells of the spatial grid configured for the simulation.

Sampling wavelengths

Each radiation source provides a mechanism to sample a random wavelength from its spectral energy distribution (SED) and assign a corresponding luminosity such that the bolometric luminosity of the source is properly distributed across the photon packets.

All currently implemented SEDs employ the inversion method using the tabulated cumulative distribution to sample wavelengths. It would be possible, however, to employ a specialized sampling routine for analytically defined SEDs.

It is important that all aspects of the SED are properly sampled, for example wavelength ranges with low luminosities or narrow line features. To this end, the sampling procedure employs composite biasing to combine sampling from two distributions: the SED itself (favoring wavelengths that carry a lot of energy) and a custom bias distribution (favoring specific wavelength areas). Both the composite bias factor and the bias distribution are user-configurable. By default, the procedure has a bias factor of 0.5 and a bias distribution that is logarithmic in wavelength.

Wavelength grids

A simulation uses multiple wavelength grids, each specialized for a particular purpose:

The current implementation requires the radiation field wavelength grid and the dust emission grid to be configured by the user. Future implementations may allow determining this grid automatically using some heuristic.

Bands

SKIRT includes built-in bands, i.e. definitions of the transmission curves for many standard instruments used in astronomy such as GALEX, SDSS, 2MASS, Spitzer, Herschel and ALMA. Bands can be used for luminosity normalization and for detecting photon packets in instruments. To enable the use in instruments, SKIRT offers a “BandWavelengthGrid” consisting of a list of (possibly overlapping) bands. When attached to an instrument, an arriving photon packet is registered for each band in the wavelength grid after multiplying its luminosity by the band’s transmission factor corresponding to the packet’s wavelength.

Radiation field

When relevant for the simulation, information on the radiation field is accumulated for each spatial cell while photon packets are traced through the spatial domain. In SKIRT 9, the stored quantities do not depend on the properties of the medium residing in the cell other than through the indirect effects of attenuation by the medium in the spatial domain at large.

This has two important benefits. Firstly, it is easier to support multiple media types in the same simulation, because the stored quantity does not directly depend on any one of the media types specifically. Secondly, the radiation field is properly accumulated even if a cell is empty, allowing output of the radiation field (e.g., for each cell or along a planar cut) to be more complete.

Indicative dust temperature

A dust grain population consisting of identical grains (material type and size) which happen to be in local thermal equilibrium (LTE) with the embedding radiation field has a well-defined equilibrium temperature. Except for this idealized case, however, there is no such thing as "the temperature" of a dust mixture. Grains of different sizes are heated to different temperatures, and when grains are not in LTE, they even show a range of temperatures rather than a single equilibrium temperature.

Despite of these considerations, it is sometimes useful to represent the energetic state of a dust mixture with a single number. For example, spatial temperature cuts may help evaluating or debugging the configuration for a new simulated model, even if the displayed temperatures do not have a straightforward physical interpretation.

For these purposes, SKIRT 9 allows calculating an "indicative dust temperature" for each spatial cell. This value is obtained by mass-averaging the LTE equilibrium temperatures for the various dust populations present in the cell in a well-defined way. Please keep in mind that the indicative dust temperature does not really correspond to a physical temperature. For more information, refer to the MediumSystem::indicativeDustTemperature() function.

Importing snapshots

Import snapshot files must be in a well-defined, SKIRT-specific column text format. There are three supported snapshot paradigmns, each with their own format for describing the spatial geometry. In addition, the input file must contain the properties required by the component being imported. The number and type of such properties depends on the component's type (source, medium, geometry) and on the component's configuration options.

The supported snapshot paradigmns and the corresponding input formats are:

In all cases, columns are in a fixed order prescribed by SKIRT. Values are either in the default units prescribed by SKIRT, or in the units defined through header lines in the data file being imported. In the latter case, any of the ski file units corresponding to the imported quantity can be used.

Simulation structure

The following simulation item hierachy illustrates the run-time structure of a typical simulation, which is reflected in the structure of the configuration file:

MonteCarloSimulation
   sourceSystem
      sources[]
   mediumSystem
      options
      media[]
      grid
   instrumentSystem
      defaultWavelengthGrid
      instruments[]
   probeSystem
      probes[]

The top-level MonteCarloSimulation item has a SimulationMode property. The value of this property determines the overall operative mode of the simulation. Here are the currently supported simulation modes; other modes may need to be added when implementing new physics:

The oligochromatic wavelength regime does not allow kinematics.

Sources

Each source defines the following information at each point in the spatial domain of the source (which is independent of the spatial grid in the simulation):

The following class hierachy illustrates the available sources and their inheritance relationships:

Source
   NormalizedSource
      GeometricSource
      PointSource
      CenteredSource
         CubicalBackgroundSource
         SphericalBackgroundSource
         StellarSurfaceSource
   ImportedSource
      ParticleSource
      MeshSource
         VoronoiMeshSource
         AdaptiveMeshSource

Note that the anisotropic geometries offered by SKIRT 8 (such as the background and surface sources) have been replaced by specialty sources. Other combinations of the generic source capabilities can be provided later as the need arises.

Media

Each medium defines the following information at each point in the spatial domain of the source (which is independent of the spatial grid in the simulation):

The following class hierachy illustrates the available media and their inheritance relationships:

Medium
   GeometricMedium
   ImportedMedium
      ParticleMedium
      MeshMedium
         VoronoiMeshMedium
         AdaptiveMeshMedium

There are three types of material: dust, electrons, and gas. Currently, no gas material types have been implemented. A material mix can represent any of these material types. Regardless of its material type, a material mix must provide the material properties required for tracing photon packets through a material of this type, in other words, for processing absorption and scattering. The currently implemented photon cycle loop can work with:

The set of additional material properties that may be needed to calculate secondary emission spectra (e.g. thermal emission from dust grains) and/or medium state updates (e.g. hydrogen ionization fraction) differs between fundamental material types.

Currently a single material mix family is implemented that allows selecting one of a set of user-configured dust mixes based on an index value imported from the medium snapshot. Other families can be implemented depending on future needs. Note however, that each dust mix (including each dust mix "spawned" from a family) pre-calculates and stores a potentially substantial amount of information, of the order of 500 MB for a typical dust mix supporting stochastic heating. Consequently, the number of dust mixes that can be spawned from a family is usually severely limited by the amount of available memory.

Geometries and decorators

Geometries focus on defining a normalized spatial density distribution and no longer support anisotropic emission. Sources requiring anisotropic emission, including point sources and radiation backgrounds, are implemented as specialty sources (see Sources). Because it is no longer possible to apply geometry decorators in these cases, the specialty sources offer an explicit position property replacing the offset decorator. Most other decorator transforms are irrelevant for specialty sources anyway.

To allow the application of geometry decorators to arbitrary density distributions, imported geometries are provided in addition to the built-in semi-analytical geometries. Imported geometries are of course subject to the limitations inherent in the geometric sources and media, which essentially means that properties other than mass density are constant across the spatial domain. Also, because geometries are normalized to unity, a luminosity or mass normalization must be supplied separately in the configuration even though the information can be deduced from the contents of the imported file.

Instruments

The instrument hierarchy has been restructured to make things more flexible and orthogonal. All instruments use the FluxRecorder helper class (not exposed to the user configuration) to actually record SEDs and/or frames. This allows offering optional functionality such as tracking output statistics in all instruments with minimal effort. Each instrument can have its own wavelength grid.

When recording of statistics is turned on for an instrument, it also outputs information intended for calculating statistical properties of the results in each bin (a wavelength bin for an SED, or the surface brightness pixel for a particular wavelength in a frame). For example, the recorded information allows calculating second order statistical properties such as the relative error R and fourth order statistical properties such as the variance of the variance VOV. Analysis of these statistical properties may help answering the difficult question of how many (or how few) photon packets are necessary to provide "acceptable" results.

The following class hierachy illustrates the available instruments and their inheritance relationships:

Instrument
   DistantInstrument            // parallel projection
      SEDInstrument
      FrameInstrument
      FullInstrument
   PerspectiveInstrument        // perspective projection
   AllSkyInstrument             // all-sky projection (configurable)

Probes

SKIRT essentially writes two kinds of output:

SKIRT 9 moves responsibility for producing the second type of output to a separate set of classes/objects, together called the probe system, as opposed to embedding it in the simulation code itself. From a user perspective, the SKIRT 8 writeXXX attributes scattered around the various simulation items are replaced by a list of probe objects living in a dedicated area of the run-time hierarchy (see Simulation structure).

Other than resulting in cleaner code, this design has the following benefits:

Currently, probes can get invoked at two points in the simulation: at the end of the setup phase, and at the very end of the simulation run. Additional probe points may be added as they become relevant in the future, for example at the end of each iteration towards self-consistency.

The following tables list the currently implemented probes.

Probes invoked after setup Description of output
InstrumentWavelengthGridProbe Text column files with the representative points and bin widths for the instrument wavelength grids in the simulation
LuminosityProbe Text column file with the primary source luminosities as a function of wavelength
LaunchedPacketsProbe Text column file with the number of photon packets launched from primary sources as a function of wavelength
SpatialGridPlotProbe Text data file(s) for plotting the grid structure
SpatialGridConvergenceProbe Text file for human consumption with theoretical versus gridded medium masses and optical depths along selected sightlines
TreeSpatialGridTopologyProbe Text data file representing the topology of the simulation’s spatial tree grid, which can be loaded by a future simulation
DefaultMediaDensityCutsProbe FITS files with theoretical and gridded medium density cuts along the default coordinate planes (depending on symmetries)
OpticalDepthMapProbe FITS file for each material type with an optical depth map seen from a given location
SpatialCellPropertiesProbe Text column file with cell center, volume, optical depth, and medium densities for each cell in the spatial grid
SpatialGridSourceDensityProbe Text column file with the normalized density of geometric source components sampled at the center of each cell in the spatial grid
OpticalMaterialPropertiesProbe Text column files listing the optical properties (cross sections, mass coefficients, asymmetry parameter) for the media in the simulation as a function of wavelength
DustGrainPopulationsProbe Text column files listing information (grain type, grain size range, mass contribution) on the dust grain populations of multi-grain dust mixes in the simulation
DustGrainSizeDistributionProbe Text column files tabulating the size distributions of the dust grain populations of multi-grain dust mixes in the simulation
DustEmissivityProbe Text column files with the emissivity of the dust mixtures in the simulation when embedded in a set of predefined radiation fields
Probes invoked after run Description of output
DefaultRadiationFieldCutsProbe FITS file with mean radiation field intensity cuts along the default coordinate planes (depending on symmetries) as a function of wavelength
RadiationFieldPerCellProbe Text column file with the mean radiation field intensity for each cell in the spatial grid as a function of wavelength
DefaultDustTemperatureCutsProbe FITS files with indicative dust temperature cuts along the default coordinate planes (depending on symmetries)
DustTemperaturePerCellProbe Text column file with the indicative dust temperature for each cell in the spatial grid
LinearDustTemperatureCutProbe Text column file with an indicative dust temperature cut along a given line segment
MeridionalDustTemperatureCutProbe Text column file with an indicative dust temperature cut along a given meridian line (i.e., one half of a great circle)
DustAbsorptionPerCellProbe Text column file with the spectral luminosity absorbed by dust for each cell in the spatial grid as a function of wavelength

Configuration user interface

SKIRT 9 offers even more options than previous versions. This is inevitable because it supports more physics (e.g. media other than dust) and has more flexibility (e.g. biasing options for launching photon packets). This extra complexity might be confusing, especially to beginning or occasional users.

To help alleviate these concerns, the configuration mechanism has been equipped with more powerful capabilities to dynamically adjust the displayed options and the default values depending on choices made earlier in the configuration process. These earlier choices may represent various aspects of the configuration; for example:

All the metadata needed to control these capabilities are defined in the SKIRT class headers using macro’s as usual.

Resources

SKIRT uses various data tables representing SEDs, SED families, optical dust properties, enthalpies, and so forth. Implementing the planned new physics will require an ever growing set of potentially very large resource files (many Gigabytes). Reading and interpreting the many different third-party formats in the SKIRT C++ code is cumbersome and slow. Therefore, in SKIRT 9, all tabular resource data are converted to a binary, SKIRT-specific format (called "stored table format") using Python procedures included with PTS. The SKIRT code then needs to support just a single format, which can be optimized, for example, to allow memory-mapped access.

Furthermore, to avoid burdening the source code repository with large data files, resources are made available on the SKIRT web site and need to be downloaded separately as part of the installation procedure.

Parallelization

SKIRT 9 implements hybrid parallelization similar to the SKIRT 8 task parallelization mode. The total number of tasks at hand (e.g. the number of photon packets to be emitted) is divided into relatively small chunks. These chunks are handed out dynamically to the various processes and/or threads in the simulation until all chunks are completed. This “on-call” scheme should provide close to perfect load balancing.

Unfortunately, the SKIRT 8 data parallelization mode has not been implemented in SKIRT 9. In this mode, SKIRT 8 distributes the largest memory data tables of the simulation across multiple processes, lowering the overall memory requirements for complex models. To achieve this, SKIRT 8 assigns a specific set of wavelengths to each process. However, in SKIRT 9, photon packets can change wavelength during their life cycle, so implementing this scheme is impossible, or at least far from straightforward.

Performance

SKIRT 9 simulations using a single optical wavelength in oligochromatic regime proceed at essentially the same speed (within 10%) as those in SKIRT 8. In this special case, any overhead caused by the more complex wavelength framework (see below) is limited and often fully compensated by the gain obtained from the improved parallel scaling.

On the other hand, simulations in panchromatic mode using the full range from UV to submm wavelengths usually run slower in SKIRT 9 by up to a factor of two compared to SKIRT 8. There are various reasons for this performance impact, and most of them are related to the new treatment of wavelengths. In SKIRT 8, the user configures a global wavelength grid for a simulation, and all photon packets carry an integer wavelength index into this grid. Instead, in SKIRT 9, for each photon packet a new wavelength is sampled from the source distribution (usually tabulated on 1000 or more wavelength points) in combination with a corresponding biasing mechanism to ensure that all spectral regions are properly sampled.

The most important performance impact, however, is caused by the need to convert the photon packet wavelength (a floating point number) to the corresponding bin index (an integer number) in the wavelength grids encountered for each step along the way in the photon packet lifecycle. This includes retrieving optical media properties, storing the radiation field in each spatial cell, and detecting photon packets in the instruments. It is essentially impossible to accelerate these binning operations through, for example, caching, because of the large amount of generality and flexibility allowed in SKIRT models. The wavelength of a packet can change throughout its lifetime (for example, caused by kinematic effects), and each medium component and instrument may have its own distinct wavelength grid.