Substantial internal update, fragmented dust mixes, dynamic medium state, and more

[petercamps]

This pull request incorporates a number of important updates that have been developed over the last few months in a separate branch. The changes include a substantial rewrite of some portions of the code to support new features, some of which are included in this update. The following list offers a summary of the changes.

Internal changes

• The absorption and scattering behavior of material mixes can now depend on any information in the medium state maintained by SKIRT for the corresponding medium component (as opposed to just the mass density) and on any information in the incoming photon packet (as opposed to just the wavelength). This allows implementing more material and interaction types without having to resort to ad hoc code.

• The number and type of variables maintained in the medium state for each spatial cell is now determined at runtime depending on the material mixes in the simulation. This results in less memory overhead for simple configurations; for example no memory is allocated for velocities or magnetic fields if these quantities are not used. Perhaps more importantly, material mixes can now request custom state variables for specific purposes. For example, an advanced dust mixture might want to import a separate density value for each of a number of grain size bins.

• As before, all medium state variables are initialized during setup, for example through importing information from snapshot files. However, it is now possible to self-consistently determine final values of the medium state variables based on the primary radiation field through an iterative process in the simulation. This can be used, for example, to self-consistently calculate dust destruction in a harsh radiation field.

• When calculating photon packet paths through the spatial grid, the segments of the path crossing each spatial cell are now generated one by one as opposed to being calculated and stored in one go. This allows efficient implementation of a photon cycle without forced scattering, where there is no need to calculate the path beyond the interaction location. This is mostly relevant for simulations with very short free path lengths, such as those including Lyman-alpha resonant scattering (see the seperate section on Lyman-alpha below).

New features related to these changes

• The new fragment dust mix decorator class splits another specified dust mix into separate fragments for each of its grain material types or even for each of a number of grain size bins per material type. The dust mix being fragmented can be a turn-key multi-grain mix or a fully user-configured dust mix, offering nearly limitless possibilities.

• A fragmented dust mix (i.e. a dust mix decorated as described in the previous item) can import a distinct density weight factor for each fragment from a snapshot input file. This allows, for example, to import dust densities from a hydro simulation that follows dust formation in a number of separate grain size bins using a single SKIRT medium component.

• There is a new set of options to enable and configure the capability of dynamically adjusting the medium state by iterating over the primary radiation field. This includes basic settings such as the minimum and maximum number of iterations. More importantly, it allows specifying one or more recipes for updating the medium state based on the radiation field at the end of each iteration.

• The new linear dust destruction recipe class implements a basic dust destruction recipe based on user-configured cutoff temperatures for silicate and graphite grains. It cooperates with a fragmented dust mix to determine the destruction fraction for each of the dust mix fragments separately based on the radiation field at the end of each iteration. Other recipes can be implemented by adding new subclasses.

• New probes allow outputting the values of custom state variables, such as the dust densities or destruction fractions for each fragment in a fragmented dust mix.

• A new advanced photon packet configuration option allows selecting a photon cycle without forced scattering, which is the default for the Lyman-alpha simulation mode (see the seperate section on Lyman-alpha below). The current implementation does not support storing the radiation field with a non-forced scattering photon cycle, so it cannot
  be used in simulations with secondary emission.

Other new features

• There are new classes representing dust grain size distributions loaded from a user-provided text column file or fully specified within the configuration file.

• The grain compositions representing spheroidal dust grains now offer built-in property tables and allow specifying an alignment fraction in addition to supporting user-provided tables. Note that support for spheroidal dust grains is limited to polarized thermal emission.

Performance

The effect of these changes on the speed of regular dust extinction and emission simulations is limited. For typical 3D simulations with a complex geometry using an octree or Voronoi spatial grid, the difference is insignificant and the new code might actually be slightly faster. For the simplest configurations with 1D or 2D symmetries the effect can be larger; worst case the new code is up to 25 per cent slower. Because these configurations are mostly used for benchmarks, this seems to be an acceptable trade-off.

Incompatibilities

• The classes representing spheroidal dust grains have been renamed and offer options with new names. Any ski files using spheroidal dust grains must be adjusted or recreated.

• Use of the builtin propery tables for spheroidal dust grains requires installing the Core resource pack version 4. No other changes were made to the resource pack so version 3 will work just fine for all other purposes.

• The class representing Lyman-alpha hydrogen has been renamed. Any Lyman-alpha ski files must be adjusted or recreated.

Testing

The new version of the code properly runs all existing functional tests, all new functional tests created to test the new capabilities, and all dust-related benchmarks published on the SKIRT web site. As a result, we believe this version of the code is ready for production use, with the exception of the Lyman-alpha features (see the section on Lyman-alpha below).

Lyman-alpha disclaimer

This update represents another progression in the support for Lyman-alpha radiative transfer through the implementation of an optimized non-forced scattering photon cycle. This speeds up the calculation of random walk photon packets in high optical depth regions. Unfortunately, the paths of peel-off photon packets sent for each scattering event to each of the instruments in the simulation still need to be calculated in full. Because a Lyman-alpha photon packet may scatter hundreds of thousands of times, this path calculation tends to dominate the processing time.

However, the Lyman-alpha capability is currently not ready for production use:

• While SKIRT seems to properly produce the results for some very basic spherically symmetric Lyman-alpha benchmark setups, we have not yet verified Lyman-alpha results for more complex geometries.

• There are indications that the results generated for complex, imported 3D geometries strongly depend on the type and resolution of the spatial grid employed in the simulation. It is at this point unclear whether this issue is inherent to the Lyman-alpha radiation transfer problem or whether it is caused by one or more errors in the code.

• In the current implementation it is quite tricky to configure a working Lyman-alpha simulation. One must properly confine the wavelength range of the sources and instruments to a narrow range around the Lyman-alpha line center and provide the appropriate wavelength biasing options.