The Monte Carlo photoionization and moving-mesh radiation hydrodynamics code CMacIonize
Bert Vandenbroucke Kenneth Wood





adres]B. Vandenbrouckecor1
bv7@st-andrews.ac.uk
adres]K. Wood
kw25@st-andrews.ac.uk

[cor1]Corresponding author
[adres]SUPA, School of Physics & Astronomy, University of St Andrews,
North Haugh, St Andrews, KY16 9SS, United Kingdom


We present the public Monte Carlo photoionization and moving-mesh radiation
hydrodynamics code , which can be used to simulate the
self-consistent evolution of HII regions surrounding young O and B stars, or
other sources of ionizing radiation. The code combines a Monte Carlo
photoionization algorithm that uses a complex mix of hydrogen, helium and
several coolants in order to self-consistently solve for the ionization and
temperature balance at any given type, with a standard first order hydrodynamics
scheme. The code can be run as a post-processing tool to get the line emission
from an existing simulation snapshot, but can also be used to run full radiation
hydrodynamical simulations. Both the radiation transfer and the hydrodynamics
are implemented in a general way that is independent of the grid structure that
is used to discretize the system, allowing it to be run both as a standard fixed
grid code, but also as a moving-mesh code.



numerical hydrodynamics radiative transfer ISM: evolution






§ INTRODUCTION



Photoionization of hydrogen and helium in the interstellar medium (ISM) by
luminous UV sources has an important effect on the evolution and properties of
the ISM. Absorption of ionizing radiation through ionization is an important
source of energy that feeds the expansion of bubbles surrounding young O and B
stars, and hence shapes the structure of the ISM on small scales
<cit.>. In cases where the dynamical effect of ionizing radiation
is less pronounced, the presence of ionizing radiation will still alter the
overall ionization balance, not only of hydrogen and helium, but also of other
elements. This, together with an increase in temperature in ionized regions,
will have a visible impact on their emission spectrum, making them stand out as
HII regions. Detailed observations of HII emission spectra contain a wealth of
information about the local ISM and the incident radiation field, and modeling
them is important in understanding observational signatures of star formation
<cit.>, and diffuse emission from galactic discs
<cit.>.

On larger scales, photoionization also has important dynamical effects. The
combined UV emission of quasars and young stars in the early Universe generates
a UV background radiation field that is responsible for the reionization of the
Universe by redshift 6 <cit.>. This UV background field
affects the formation of galaxies by altering the abundances of ISM coolants
<cit.>, and is responsible for suppressing galaxy formation in low
mass haloes <cit.>. Futhermore,
radiative feedback might be an important mechanism to regulate star formation in
galactic discs <cit.>.

Modeling photoionization and in particular HII regions requires solving a very
complex system of ionization balance equations for the various elements present
in these regions, which is only possible if strong assumptions are made. The
widely used code Cloudy ascl:9910.001 <cit.> for example
assumes a simple 1D geometry, but keeps track of a large number of elements and
ionization stages. When a real 3D geometry is necessary, it is no longer
feasible to keep track of so many elements, and a selection has to be made,
depending on the problem at hand.

When the effect of ionizing radiation on the dynamics of the ISM is studied,
further assumptions need to be made about how to deal with the coupling between
radiation transfer and hydrodynamics, to get a radiation hydrodynamics (RHD)
scheme. Some methods treat the radiation field as a fluid governed by diffusion
equations <cit.>. These methods have the advantage that
they do not require extensive algorithmic changes and are relatively efficient.
However, they have undesired side effects, like e.g. the absence of shadows in
optically thin regions, and the fact that extra assumptions need to be made
about the propagation speed of the radiation field to prevent the integration
time step from getting very small. Alternative methods use an approximate ray
tracing scheme <cit.> which is more
complex to implement and preserves some directional information. However, these
schemes require careful fine tuning to make sure the radiation field accurately
covers the density structure, especially if the density structure is asymmetric
or clumpy.

A more accurate, albeit less efficient way to treat the radiation field is
provided by using Monte Carlo based RHD codes like torus ascl:1404.006
<cit.> and mocassin ascl:1110.010 <cit.>.
These codes have the advantage that they are also much more flexible and easier
to extend with extra physics, e.g. extra chemistry <cit.>.
Furthermore, Monte Carlo techniques are also widely used to model dust
scattering and absorption <cit.>, making it straightforward to
include dust scattering in Monte Carlo based RHD modeling.

In this work, we present our own Monte Carlo RHD code called , that
couples a basic finite volume hydrodynamics scheme with a Monte Carlo
photoionization code. Our code can use a variety of different grid types to
discretize the ISM, and can run with both a fixed grid and a fully adaptive
moving mesh. Apart from running as a RHD code,  can also be run as
a pure Monte Carlo photoionization code, and can be used to post-process density
fields from other simulations.

Our code is written in modular C++ and is meant to be both user-friendly and
efficient by combining a well-structured and documented design with an
implementation that makes use of new features of modern C++11. The code has a
limited number of dependencies and can be run in parallel using a hybrid OpenMP
and MPI parallelization strategy. Some parts of the code are wrapped into a
Python library using Boost
Python[<http://www.boost.org/doc/libs/release/libs/python/doc/html/index.html>].
The photoionization part of the code can also be used as an external C or
Fortran library, facilitating coupling our code to other simulation codes.

This paper is structured as follows: in section <ref> we discuss
the physics that has been implemented in , and give a short
overview of the Monte Carlo photoionization technique and the finite volume
hydrodynamics scheme. In section <ref> we describe the design
considerations that were used during the development of the code, and detail
their implementation. We conclude in section <ref> with the
results of a number of benchmark tests that are part of the public code
repository and that show its accuracy and performance.




§ PHYSICS




The emission line spectrum of a star forming nebula is determined by its thermal
equilibrium, which is a steady-state equilibrium between heating through
photoionization by UV sources, and cooling by various atomic processes in the
nebula. <cit.> identify four important sources of cooling:

  
  *  Energy loss by recombination of hydrogen and helium, i.e. the reverse
  of the photoionization process,
  
  *  energy loss by bremsstrahlung emitted by free electrons,
  
  *  energy loss by collisionally excited line radiation from some abundant
  metals, and
  
  *  energy loss by collisionally excited line radiation from hydrogen.

In order to compute photoionization and recombination rates, we need to know the
ionization structure of the gas in the nebula, and the temperature of the
nebula. Since the temperature itself is the solution of the thermal equilibrium,
this can only be solved for iteratively.

The thermal and ionization equilibrium is also important for the dynamics of the
nebula: ionized regions have more free particles and hence a higher specific
energy than neutral regions, so that photoionization effectively acts as a
heating term in the hydrodynamics of the gas. In order to properly model this
effect, combined radiation hydrodynamics (RHD) simulations are necessary.

 can be run in two different modes: either as a pure Monte-Carlo
photoionization code that ray traces the radiation of an ionizing UV radiation
field through a density field and self-consistently solves for the ionization
and temperature structure, or as a radiation hydrodynamics code that uses the
output of the photoionization code as a heating source in a hydrodynamical
simulation. The former is essentially a completely rewritten version of the
photoionization code of <cit.>, while the latter combines this code
with a standard finite volume method which is a simplified version of the
algorithm implemented in  ascl:1605.003 <cit.>.
We will summarize the most important physical ingredients of both methods below.

Note that in the current version of the code, we do not include a treatment of
non-ionizing radiation, nor do we take into account dust scattering and the
dynamical effect of radiation pressure on dust. The treatment of these processes
uses algorithms that are very similar to the ones used for photoionisation, and
it is straightforward to extend the code with these processes in the future.




 §.§ Photoionization


As our initial research focusses on diffuse ionized gas in star forming nebulae,
we only model photoionization of hydrogen and helium self-consistently, for UV
radiation in the energy range [13.6, 54.4] eV, corresponding to the ionization
threshold for hydrogen and the second ionization threshold for helium. As in
<cit.>, we do not trace double ionized helium, and we only care about
photons that are energetic enough to ionize hydrogen.

To model the various cooling mechanisms correctly, we also need to know the
ionization structure of a number of coolants, i.e. C^+, C^++, N^0,
N^+, N^++, O^0, O^+, O^++, Ne^+, Ne^++, S^+, S^++, and
S^+++ (see <ref>). These are treated approximately,
where we make the assumption that the number of free electrons released by
photoionization of these elements is neglible compared to the total number of
free electrons, which allows us to use a simplified ionization balance equation.

This approximation only holds in regions that are sufficiently ionized, as the
total number of free electrons is mainly determined by the ionization of
hydrogen and helium for realistic elemental abundances.

Note that <cit.> did not include cooling due to S^+++, and does
not mention the use of carbon cooling rates (although they were used). However,
we found that not including these coolants leads to excessively high
temperatures in the Lexington benchmark tests (see
<ref>).



  §.§.§ Monte Carlo technique


The local photoionization rate depends on various factors, the most important of
which are the position, direction and energy of the incoming ionizing radiation,
and the local ionization state. Due to the strong non-linearity of the
photoionization process, it is impossible to exactly solve for the ionization
balance except for a very limited number of cases, so that approximate
techniques are required.

As a first step, we discretize the density field of interest on a geometrical
grid structure consisting of a (large) number of small cells. Each cell contains
a compact subset of the total physical region of interest and is bounded by a
discrete number of planar faces, which separate it from neighbouring cells. Our
grid can be a regular Cartesian grid consisting of cubical cells, but can also
be a hierarchical adaptive mesh refinement (AMR) grid <cit.>, or an
unstructured grid <cit.>.

The radiation field is also split up into a (very large) number of photon
packets with a certain weight, which are sampled using a Monte Carlo technique.
Each packet represents a fraction of the total radiation field, and is emitted
by randomly sampling its properties (origin, travel direction and
energy/wavelength) from underlying distribution functions. The photon packet is
then ray-traced through the density grid by computing geometric path lengths
until a randomly sampled optical depth is reached, or the photon packet leaves
the simulation box. For each discrete cell that lies in the path of the photon
packet, we keep track of the total path length covered by the photon packet
within that cell, to get an approximation to the photoionization integral
I_i,X^0 → X^+ for photoionization from ion X^0 to level X^+
in that cell <cit.>:

    I_i,X^0 → X^+   =   ∫_ν_X^0 → X^+^∞4πJ_i(ν')/hν'σ_X^0 → X^+(ν')
    dν' 
       ≈   Q/W_tV_i∑_j w_j l_i,jσ_X^0 → X^+(ν_j),

where J_i(ν) is the mean intensity of radiation in the cell as a function
of frequency, ν_X^0 → X^+ is the threshold ionization
frequency for ionization of ion X^0, σ_X^0 → X^+(ν)
is the ionization cross section as a function of photon frequency, Q is the
total luminosity of all UV sources, V_i the volume of the cell, w_j is the
weight of the individual photon packets j that pass through the cell (with
W_t = ∑_j w_j the total weight of all packets), each of which covers a path
length l_i,j through the cell, and
h=6.626×10^-34 J s is Planck's constant.

Note that the original code of <cit.> assumed equal weights for all
photon packets. We generalized this to improve the sampling of low luminosity
external radiation fields.

To account for the diffuse radiation field caused by recombination of ionized
hydrogen and helium, we perform an extra sampling step when a photon packet has
reached the desired optical depth and is still within the simulation box. We
first randomly decide if the photon is absorbed by hydrogen or helium, with the
probability for absorption by hydrogen given by <cit.>

    P_i,j(H^0) = n_i,H^0σ_H^0 →H^+(ν_j)/n_i,H^0σ_H^0 →H^+(ν_j) + n_i,He^0σ_He^0 →He^+(ν_j),

where n_i,H^0 and n_i,He^0 are the number densities of
neutral hydrogen and neutral helium in the cell respectively.

Depending on the element that absorbed the photon, there are various reemission
channels, some of which give rise to ionizing UV radiation:

  
  *  hydrogen Lyman continuum radiation,
  
  *  helium Lyman continuum radiation,
  
  *  19.8 eV radiation from the resonant 2^3S → 1^1S
  transition in neutral helium,
  
  *  ionizing radiation for one of the two photons in the helium two photon
  continuum,
  
  *  helium Lyman α radiation.

Each reemission channel has a specific probability associated with it, which
will depend on the local temperature in the grid cell <cit.>. We use
these probabilities to randomly pick a channel. For each channel there is an
associated probability of actually producing an ionizing photon, and an
associated spectrum for the resulting reemitted radiation. We randomly decide if
the photon is reemitted in the ionizing part of the spectrum. If it is, we
sample a new random direction and optical depth for the photon and repeat the
ray-tracing step. If it is not, the photon is assumed to be reemitted as
non-ionizing line continuum which escapes from the system, and it is terminated.

As in <cit.>, we assume a medium that is optically thick to Lyman
alpha radiation, so that we do not explicitly ray-trace helium Lyman α
photons, but assume that they are absorbed on the spot. We need to explicitly
take this into account as an extra term when solving the ionization balance
within the cell:

    n_i,H^0 I_i,H^0 →H^+ +
    P_i(H_OTS) n_i,e n_i,He^+α_He^0,2^1P
     = n_i,e n_i,H^+α_H^+ →H^0(T_i),

where α_He^0,2^1P is the recombination rate of the 2^1P
level of neutral helium, α_H^+ →H^0(T)
is the recombination rate from ionized to neutral hydrogen as a function of
temperature, and T_i is the temperature in the cell.

n_i,e is the electron density in the cell, which is approximately given by

    n_i,e = n_i,H^+ + n_i,He^+

if we neglect the free electrons released by ionized metals.

P_i(H_OTS) is the probability of on the spot absorption of
helium Lyman α radiation, which is approximately given by
<cit.>

    P_i(H_OTS) = ( 1 +
    0.77/√(T_i/10^4 K)f_i,He^0/f_i,H^0)^-1,

with f_i,H^0 and f_i,He^0 the neutral fractions of
hydrogen and helium in the cell, defined as

    f_i,X^0 = n_i,X^0/n_i,X^0 + n_i,X^+.


For helium the ionization balance is simply given by

    n_i,He^0 I_i,He^0 →He^+ = n_i,e
    n_i,He^+α_He^+ →He^0(T_i).


Equations (<ref>) and
(<ref>) are solved simultaneously for a given input
temperature T_i. Once the ionization state of hydrogen and helium is known, we
can compute the density of free electrons using equation
(<ref>). With these densities, we can solve for the
ionization state of the metals, for which the ionization balance is generally
given by

    n_i,X^0 I_tot,i,X^0 → X^+(T_i) = n_i,X^+α_tot,i,X^+ → X^0(T_i),

with I_tot,i,X^0 → X^+(T_i) and
α_tot,i,X^+ → X^0(T_i) the total ionization and
recombination rate in the cell.

The ionization rate is generally given by

    I_tot,i,X^0 → X^+(T_i)    = I_i,X^0 → X^+
        + n_i,H^+ I_C,H,X^0 → X^+(T_i) 
        + n_i,He^+ I_C,He,X^0 → X^+ (T_i),

where I_C,H,X^0 → X^+(T) and I_C,He,X^0
→ X^+(T) are the ionization rates for ion X^0 through a charge
transfer reaction with ionized hydrogen or helium respectively, as a function of
temperature.

Likewise, the total recombination rate is given by

    α_tot,i,X^+ → X^0(T_i)    = n_i,eα_X^+
    → X(T_) 
        + n_i,H^0α_C,H,X^+ → X(T_i) 
        + n_i,H^0α_C,He,X^+ → X(T_i),

with α_C,H,X^+ → X(T) and
α_C,He,X^+ → X(T) charge transfer recombination
rates.

We will generally not include charge transfer ionization rates for reactions
involving ionized helium, and only include charge transfer ionization rates for
hydrogen and charge transfer recombination rates for some of the metals.



  §.§.§ Data




  
Spectra

We support a number of input spectra for the ionizing radiation, ranging from
single frequency spectra that are used for benchmark tests, over black body
spectra, to realistic stellar atmosphere spectra from the models of
<cit.>. For complex spectra, we precompute the cumulative number
distribution function of ionizing photons for a discrete number of frequencies,
and then use linear interpolation on this table to sample random frequencies at
runtime.

We also support external radiation fields, like a redshift-dependent cosmic UV
background field which can be used to model the ISM of high redshift galaxies.
For this, we use the spectra of <cit.>, as downloaded from
their website.



  
Ionization cross sections

We use fits to the photoionization cross sections of hydrogen, helium, and the
various coolants from <cit.>. These fits smooth out over
resonances, but have the advantage that they are relatively cheap to compute at
runtime.



  
Recombination rates

We use radiative recombination rate fits of <cit.>. For the
coolants, these are supplemented with dielectronic recombination rate fits of
<cit.>, <cit.>, <cit.>, and
<cit.>.



  
Charge transfer rates

We use the charge transfer ionization and recombination rates for hydrogen from
<cit.>, and the helium charge transfer recombination rates of
<cit.>.



  
Alternative data

For some of the benchmark tests described in section <ref> we
need simplified values for the photoionization cross sections and radiative
recombination rates. To this end, we made sure that the rates can be easily
changed in our implementation, as part of the modular design of our code (see
<ref>).




 §.§ Heating and cooling


When a photon packet is absorbed by hydrogen or helium, an amount of energy
equal to the excess w.r.t. the ionization energy of that element is converted
into heating of the local gas. The total integrated heating H_i,X^0
→ X^+ for a grid cell is given by <cit.>

    H_i,X^0 → X^+ = ∫_ν_X^0 → X^+^∞4πJ_i(ν')/hν'σ_X^0 → X^+(ν') 
    
    h ( ν' - ν_X^0 → X^+) dν'.

and is hence very similar to the ionization rate estimator, and can be treated
in the same way during our Monte Carlo photon propagation scheme. As for the
ionizing luminosity, the optical depth for a cell depends on the temperature and
ionization state of the gas in that cell, so that a self-consistent solution can
only be obtained with an iterative scheme.

For cooling by recombination of hydrogen and helium, we use the recombination
cooling rates of <cit.>. The cooling rate due to bremsstrahlung from
free electrons is given by <cit.>, where we use the fits to the
mean Gaunt factor given by <cit.>.




 §.§ Line cooling



Despite the low abundances of metals such as C, N, O, Ne, and S in star forming
nebulae, line emission from these elements contributes signicantly to the
radiative cooling, as their low-lying energy levels can be easily excited
through collisions with free electrons <cit.>. To model this
process, we need to keep track of the ionization state of these coolants, and
model their line emission. The details of this treatment are explained below.



  §.§.§ Mechanism and data


In general, the level population x_X,i = n_X,i/n_X of the ith
energy level of an ion X with density n_X is the solution of
<cit.>

    ∑_j > ix_X,jn_e q_X,ji + ∑_j > i x_X,j A_X,ji
    
    = ∑_j < i x_X,i n_e q_X,ij + ∑_j < i x_X,i A_X,ij,

where q_X,ij is the collisional excitation or deexcitation rate from level
j to level i, A_X,ij is the radiative deexcitation rate from level i
to level j, and n_e is the electron density.

The collisional deexcitation rate q_X,ji is given by

    q_X,ji = h^2/√(k)(2πm_e)^3/2Υ_X(i, j)/ω_X,j,

with k=1.38 × 10^-23 J K^-1 Boltzmann's constant,
m_e = 9.109 × 10^-31 kg the mass of an electron,
Υ_X(i,j) the velocity-averaged collision strength, and
ω_X,j the statistical weight of the jth level of ion X. It is
linked to the collisional excitation rate q_X,ij through the relation

    q_X,ij = ω_X,j/ω_X,i q_X,jiexp(-χ_X,ij/kT),

with χ_X,ij the energy difference between level i and level j.

Velocity-averaged collision strengths, radiative recombination rates, energy
differences, and statistical weights can be measured experimentally or modelled
quantum mechanically. We use data from a large number of sources, as detailed in
table:line_cooling_data.



The specific ions we use can be classified into two categories: ions with two
low lying energy levels (N^++ and Ne^+), and ions with five low lying
levels (C^+, C^++, N, N^+, O, O^+, O^++, Ne^++, S^+ and
S^++). For the former, we solve the system of two equations
((<ref>) for i = 2, and x_X,1 + x_X,2 = 1)
analytically, while for the latter we need to solve the full set of five coupled
linear equations.



  §.§.§ New data fits


The velocity-averaged collision strengths used above vary with temperature. To
account for this fact, we fitted a general curve of the form

    Υ_X(i,j,T) = T^a_ij + 1( b_ij + c_ij/T + d_ijlog(T) + . 
    . e_ij T ( 1 + (f_ij - 1) T^g_ij) )

to all data, where T is the temperature (in K) and a_ij, b_ij,
c_ij, d_ij, e_ij, f_ij and g_ij are fitting parameters. The
general form of this curve was inspired by <cit.>. It is worth
noticing that there are no reliable data for some of the fine structure levels
of O^+ above 10,000 K, so that we needed to extrapolate the low temperature
data. The same is true for N^++ above 40,000 K. For C^++ and Ne^+
the data points are sparse, so the fits are less reliable. For S^+ and
S^++ the values of e, f and g were kept zero, as this provided a
better fit than when they were allowed to vary.

The values of the fitting parameters are listed in
table:collision_strength_fits_carbon
table:collision_strength_fits_sulphur, the corresponding fits and relative
differences between fit and data are shown in
figure:CII_fitfigure:SIV_fit.




 §.§ Hydrodynamics


Hydrodynamical integration is performed using a finite volume method
<cit.> on a generally unstructured, (co-)moving mesh
<cit.>. The Euler equations of hydrodynamics in conservative form

    ∂U/∂t + ∇⃗.F⃗( U ) = 0,

with

    U = [  ρ; v⃗;  P ],


    F⃗(U) = [                      ρv⃗;           ρv⃗v⃗ + P 1⃗⃗⃗; ρ( u + 1/2|v⃗|^2 ) + Pv⃗ ],

and an adiabatic equation of state P = (γ-1)ρu, where
ρ is the mass density, v⃗ the flow velocity, P the pressure, and
u the thermal energy per unit mass, are discretized on an unstructured Voronoi
mesh. Integration of the primitive variables over the volume V_i of each cell
leads to a set of conserved variables (mass m_i, momentum p⃗_i, and
total energy E_i) for each cell. Integrating the conservation law
(<ref>) over the volume of the cell allows us to
reduce the time integration of these conserved quantities as a flux exchange
between the cell and its neighbouring cells:

    dQ_i/dt = -∑_j F⃗_ij . A⃗_ij,

where A⃗_ij represents the oriented surface area of the geometrical
face between cell i and cell j, F⃗_ij is an appropriate estimate of
the fluxes (<ref>) at a representative location on the face, and
Q_i is the vector of conserved variables for that cell.

To obtain appropriate fluxes, we use the values on both sides of the face as
input for an exact Riemann solver, which gives the exact physical solution for
the two state problem defined by the variables at both sides up to a desired
precision. By sampling this analytic solution we obtain new values for the
primitive variables that can be used to compute fluxes to be used in
(<ref>).

Note that our current implementation is only first order in space and time. It
is however very straightforward to extend this to higher order by the
introduction of appropriate spatial gradients, see e.g.
<cit.>; this will be implemented in future versions of the
code. Our current implementation also does not yet include external forces or
self-gravity.



  §.§.§ Moving mesh scheme


For the specific case of an unstructured Voronoi grid, we can make the method
Lagrangian by allowing the generators of the Voronoi grid to move in between
time steps of the integration scheme. To account for this movement, we add
correction terms to the flux expressions (<ref>), as detailed in
<cit.>. Note that the
hydrodynamical integration is completely independent of the movement of the
generators. If the generators do not move, the method is Eulerian. If on the
other hand the generator movement is set to the local fluid velocity in the
corresponding cell, the method is fully Lagrangian. In the latter case, we can
solve the Riemann problem across the cell faces in the rest frame of the faces,
which leads to better accuracy in the presence of large bulk velocities, and
allows us to use a larger integration time step than would be used in an
equivalent Eulerian scheme, as the time step only depends on the relative
velocity of the fluid w.r.t. the grid.



  §.§.§ Radiation hydrodynamics


To couple the radiation to the hydrodynamics, an operator splitting method is
used, whereby the photoionization heating term is added after each step of the
hydrodynamics scheme, assuming photoionization equilibrium. The current values
of the density in each cell are converted into number densities that are then
used as input for the photoionization code. The photoionization code computes a
self-consistent ionization structure for each cell, which is then used to decide
what the temperature of the cell should be. The number of iterations and number
of photon packets used in this step is a simulation parameter. The resulting
temperature is compared with the actual hydrodynamical temperature of the cell,
and the energy difference is added to the cell as an energy source term.

In our current implementation, we do not subtract energy for cells that were
ionized and become neutral again, since this involves a careful treatment of the
time step to prevent negative energies. We use an almost isothermal equation of
state with an adiabatic index γ=1.0001, and assume a two-temperature
medium, whereby the temperature in the neutral phase is assumed to be T_n =
100 K, while the temperature in the ionized medium is assumed to be T_i
= 10,000 K. For a cell with hydrogen neutral fraction x_H, the
assumed average temperature is then

    T = x_H T_n + ( 1 - x_H) T_i.

This approach is sufficient to reproduce a basic benchmark expansion test (see
<ref>). We could of course also get the photoionization
heating directly from the photoionization step; this more general approach will
be subject of future work.

Note that it is possible to couple  as an external library to other
hydrodynamics codes using a very similar approach, see
<ref>.




§ CODE




The most important difference between  and the original
photoionization code of <cit.> is a complete redesign of the structure
of the code (accompanied by a migration from Fortran to C++11), including a full
in-line documentation of the code using
Doxygen[<http://www.doxygen.org>] (the documentation for the latest
stable version of the code can be found on a dedicated
webpage[<http://www-star.st-and.ac.uk/ bv7/CMacIonize_documentation/>]).
Below we outline the main design considerations and detail how they are
implemented in the code.




 §.§ Design considerations




  §.§.§ User friendliness


The photoionization part of  is primarily focused on
post-processing output from other simulation codes as part of a simulation
analysis work flow. This means that the code will be used by researchers that
are not necessarily very familiar with the details of the code, but that still
want to produce reproducible science products. To accommodate this, we aim to
minimize the learning curve for using . In cases where writing
additional code is inevitable (like for example when reading a new type of input
file), we want to limit the effort necessary to accomplish this: the new code
should be limited to a single function or C++ class, and the user should be able
to write this code without worrying about the details of the photoionization
algorithm.



  §.§.§ Reproducibility


A properly designed computer algorithm should be deterministic, so that running
the same simulation twice with the same input and the same version of a code
should produce exactly the same result, independent of the hardware
architecture. Parallelization and system specific optimizations might cause tiny
changes in round off that could cause minuscule changes in result between
different runs (especially in Monte Carlo algorithms), but even then a
simulation code should be very close to a unique mapping from input data to an
output solution. Reproducibility is hence an inherent feature of computer
simulations.

However, keeping track of input parameters and code versions can be tedious,
especially when simulations are combined with the design of improved algorithms
and code changes are made. For this reason, we aim to provide a robust system to
log parameters and code versions, so that all published  results
should be perfectly reproducible.



  §.§.§ Modularity


Complex algorithms combine a large number of components that each have their
specific complexities. However, most of these components are predominantly
independent of each other, and only interact with the other components through
narrowly defined interfaces. Isolating complex components into separate entities
or modules increases the readability of a code, and makes the code more robust
if combined with an appropriate unit testing strategy. We will therefore aim to
produce a modular code, whereby separate components are identified and isolated
into separate code entities.



  §.§.§ Scalability


Modern computing architecture is highly parallel, with the computing power of a
typical high performance computer spread out over a large number of separate
computing units or nodes that are interconnected through a high speed
network, and with each of these nodes in turn consisting of up to 128 separate
CPU cores that share a single memory space. In order to efficiently use
these machines, it is crucial that an algorithm is designed with a parallel
mindset. We cannot think of the algorithm as a serial list of instructions that
are executed one by one, but instead need to think of the tasks that need
to be performed by the algorithm, the data that is needed to perform
these tasks (and that might be shared with other tasks), and the
dependencies that govern which tasks can be executed in parallel and
which tasks are mutually exclusive due to conflicts.

In the current version of , we aim to provide a reasonable scaling
by distributing computations across multiple cores on a single node, and across
multiple nodes. Our current parallelization strategy does not address the need
to distribute data across multiple nodes in order to efficiently use the
available memory. This will be addressed in future versions of the code.




 §.§ Design implementation




  §.§.§ User interface



For standard users that do not plan to add additional code, the interaction with
 is limited to

  
  *  calling the command line program CMacIonize, and
  
  *  writing a parameter file that contains the parameters for the
  simulation.

The command line program has a very limited set of options that control the
number of shared memory parallel threads used to run, and the mode in which to
run (photoionization only or RHD). The only other parameter to the program is
the name of the parameter file. An example command line call to 
could like this:

This will run  in the standard photoionization mode using 8 shared
memory parallel threads, and using the parameters provided in the file
parameterfile.param.

The parameter file contains all information needed to set up and run the
simulation, and maps to the underlying modular structure of the code (all
parameters are linked to a specific C++ class). It is a simple text file in YAML
format[<http://yaml.org/>], and a very basic example could look
like this:

This parameter file sets up a Strömgren test in a box of 10 pc×10  pc×10  pc containing gas with a
density of 100 cm^-3 at a temperature of 8,000 K, with a star at the
centre with a total luminosity of 4.26×10^49 s^-1 and a
monochromatic spectrum with a frequency equivalent to a photon of
13.6 eV. The simulation uses a Cartesian grid of
64×64×64 cells, and uses 20 iterations with 10^6 photon
packets for each iteration (see <ref> for details of
this test).

The example above illustrates how easy it is to read and understand a parameter
file. It also illustrates another key feature of the code: the use of units.
Internally, we consistently use SI units throughout the code to avoid any
confusion about units. However, we also require the user to specify units for
all physical quantities that are used as an input, so that the user does not
need to worry about unit conversions at all. We support a variety of different
units, including complex unit conversions (e.g. photon energy in eV to photon
frequency in Hz), and adding new units is very straightforward.

Apart from supporting units, the parameter file also supports various number
formats and 3D vectors.

When the program is started, the parameter file is parsed and translated into a
corresponding simulation structure. Before the actual simulation starts, the
actually used parameters are written to a reference file. Most parameters have
default values and need not be specified in the parameter file; when written to
the reference file, the default values will be displayed, and the file will
clearly state that the default value was used. If a parameter is not used, the
reference file will mention this as well. For parameters that have units, the
reference file will contain the value in SI units, as well as the original
value.

If no parameter file is given, default values will be used for all parameters
(that correspond to a low resolution version of the Strömgren benchmark
test, see <ref>). In this case it is possible to use
the reference file as a first guess for sensible parameter values, and iterate
on it to construct an actual useful parameter file. All available parameters are
also extensively documented in the Doxygen documentation of the corresponding
classes.



  §.§.§ Reproducibility



In order to guarantee reproducibility, we use a strategy that consists of three
pillars:

  
  *  version control as a way to uniquely identify a specific code version,
  
  *  parameter logging in output files as a way to keep track of used
  parameters, and
  
  *  unit testing to guarantee the same results across code versions.




  
Version control

The code is stored in a public online Git
repository[<https://github.com/bwvdnbro/CMacIonize>] to make it
easier to keep track of the code history, and to facilitate collaboration on the
code. Git keeps track of the changes that are made to the code in between so
called commits, i.e. logged checkpoints of the code status. Each commit
has an associated key that uniquely identifies it, and it is always possible to
return to a specific version of the code using the appropriate commit key.

Moreover, Git also provides a command line tool called git describe that
can be used to check the current version of the local copy of the repository
that a user is using, and that checks if the repository is dirty, i.e.
contains uncommitted code changes.

We have incorporated git describe into our code configuration chain, so
that the compiled code knows (a) what the unique commit key of the current
version is, and (b) if the current code version is exactly equal to that code
version, or contains uncommitted changes. If the code contains uncommitted
changes, it will refuse to run any simulation, so that the user is forced to
commit changes and make sure the code is identifiable before running scientific
simulations.



  
Logging

The output of any simulation consists of snapshots, i.e. dumps of specific
quantities at some time during the simulation.  supports various
types of snapshot files, the default being the same HDF5 format that is also
used by Gadget2 ascl:0003.001,
SWIFT[<https://gitlab.cosma.dur.ac.uk/swift/swiftsim>],
and .

In order to exactly reproduce a snapshot file, we need to know

  
  *  what version of the code was used to generate it, and
  
  *  which parameters were used to run the code.

The former can be easily realised by storing the unique commit key for the
current code version in the snapshot files. To guarantee the latter, we also
store all parameters in the snapshot files, as they were used. This corresponds
to the values that are part of the reference parameter file (see
<ref>). By also storing the parameters for which
default values were used, we guarantee reproducibility across different code
versions, if at some point the default value for a parameter were to change.
The public version of the code contains a Python script that can extract the
parameter file that was used from a snapshot file in the default HDF5 format.

For simulations that use input from external files, as e.g. simulations that
post-process the density field from another simulation, we also need these
external files in order to reproduce the results. Since these files can be quite
large, storing them as part of the snapshot files is not an option. In this
case, we rely on the user using a convenient method of keeping track of those
files to guarantee reproducibility.

Apart from the code version and the parameter values we also log configuration
flags and system specific information. This is not strictly necessary in order
to reproduce simulation results, but might nonetheless be helpful e.g. during
debugging.



  
Unit testing

One of the key issues when developing a simulation code is making sure that the
results are scientifically accurate, and that they stay accurate throughout the
further development of the code. Unit testing is a very powerful tool to achieve
this, especially when combined with a good modular design (see
<ref>).

A unit test is a small independent program that calls a small part of the code
with known input values and checks its output against a known solution. If the
output matches the expected result, the test passes and we know that part of the
code behaves as expected.

When properly designed, a unit test covers all possible paths through the code
that is being tested, e.g. if the code contains conditions that check for
strange input values, the test will call the code with strange input values and
check that these are handled correctly.

Unit tests were an integral part of the early development process of
, as all new code was tested against reference values of the old
code of <cit.>. During the addition of new features that were not part
of the old code, we still tried to start from the unit test as much as possible,
which meant first thinking about what the expected behaviour of a code component
(function or class) should be, before actually implementing it. The overhead
this implies is quickly recovered by the ease with which we can locate bugs in
new code.

The current version of the code contains almost 70 separate unit tests, which
are managed as part of our code configuration and run using CTest, the
CMake[<https://cmake.org/>] unit testing framework. Depending on
the hardware and system configuration, the tests take less than a minute to a
few minutes to complete, and can be run as part of the standard compilation
process.

When new code is added to the stable version on the Git repository, the new code
is automatically compiled with a number of different compilers on different
systems using the continuous integration environment Travis
CI[<https://travis-ci.org>], and the unit tests are run. Code is
only allowed to be merged into the stable repository if it passes all the tests.
This way we ensure that new code never breaks or alters old functionality,
unless this is done on purpose (in which case the corresponding unit test needs
to be modified).



  §.§.§ Object oriented design



To implement modularity in our code design, we use C++ objects as the building
blocks of the code. An object has a limited number of responsibilities, and is
as unaware as possible of the rest of the code (unless the interaction of
various objects is the responsibility of the class). Most objects are covered by
a corresponding unit test (see <ref>), although
the unit tests for some classes are grouped together if this makes more sense.

We use a number of different design patterns <cit.>. Basic simulation
components like the density grid and the source distribution use
inheritance combined with a factory class to provide different
interchangeable implementations (e.g. the density grid can be a regular
Cartesian grid, an AMR mesh, or an unstructured Voronoi grid). The density grid
itself makes extensive use of iterators to provide grid type unaware
access to cells, while most grid computations are performed using
visitors that perform a single task for each cell of the grid.

The object oriented design is tightly interwoven with the parameter file used as
the user interface (see <ref>), with objects
mapping to specific blocks in the parameter file (and the  keyword
always referring to a factory class that provides multiple implementations of a
general interface). Most objects have a so-called parameter file constructor,
which can create an object instance based on the parameter values given in the
parameter file, with parameters mapping directly to object properties.



  §.§.§ Task based design





To provide inherent parallelism in our code, we think of the actions that need
to be performed by the algorithm in terms of small tasks, that perform a
limited set of actions on a small part of the computational domain. By limiting
the amount of work done by a single task, we can improve the load balancing
between different parallel threads significantly, while limiting the fraction of
the computational domain that is affected minimizes conflicts.

The photon traversal algorithm e.g. can be done independently for small batches
of photons. Each part of the path of a photon is only a single cell, so that we
only need to worry about two threads accessing the same cell at the same time to
prevent conflicts. Similarly, the ionization balance computation for different
cells can be done completely independently, and so can the temperature balance
computation or even the cell initialization.

In practice, our task based design is implemented using a number of interacting
classes called Worker, Job, and JobMarket. Job and
JobMarket make use of compile time polymorphism and are template
interfaces, meaning that these classes do not actually exist, but are abstract
interfaces that define common functionality for classes that can be used as C++
template arguments for other classes. This offers the same flexibility as run
time polymorphism, but without the computational overhead. A Worker is
our abstract representation of a thread, while a Job is the abstract
representation of a task that needs to be performed. The JobMarket is
responsible for spawning Job instances. The Worker instances in
turn are spawned by a WorkDistributor, which is the only class that needs
to know about the underlying parallel environment that is used. Our current
implementation only supports OpenMP, but it would be straightforward to replace
this by e.g. a POSIX threads or Intel Threading Building Blocks implementation.

The general workflow for a shared memory parallel run is illustrated in
figure:workdistributor. When a parallel part of the execution is
started, a corresponding JobMarket implementation is created and passed
on as a template argument to the WorkDistributor. The
WorkDistributor then generates a number of Worker instances that
are run in parallel, and that call the get_job() function of the
JobMarket instance to get actual Job instances that need to be
executed. The worker then calls the execute() method of the Job
instance to perform the task at hand. When the task is finished, the
Worker goes back to the JobMarket to get the next Job,
until no more tasks are available.





This paradigm nicely divides the responsibilities of the parallelization
process: the Job provides the actual task at hand, the JobMarket
regulates how tasks are divided and hence controls the load balancing, while the
WorkDistributor is responsible for handling the underlying parallel
environment. figure:photon_shoot_jobmarket and
figure:ionization_photon_shoot_job show how this works for a
photoionization simulation.

We need a locking mechanism to protect common variables in the shared memory
domain, like e.g. cells in the photon traversal algorithm or counter values
within the JobMarket. To this end, we provide our own Lock class
that is unaware of the underlying parallel environment. We also experimented
with using atomic data operations, but found this to be slower than a locking
mechanism in most cases, predominantly because of the lack of hardware support
for floating point atomic operations, and because the large number of variables
updated per cell access in the photon traversal algorithm reduces the impact of
using an expensive locking mechanism.




 §.§ Library exposure



To improve the usability of , we also provide a library interface
to the code, which can be used in both C, C++ and Fortran2003. This library
interface is in all ways equivalent to the standard  program (the
same IonizationSimulation class is used to run an actual photoionization
simulation), but uses special input and output classes to directly obtain the
density field from another code and return the resulting ionization structure to
that code without the need to output anything to disk.

Our current library implementation has already successfully been used to couple
the code to the SPH code PHANTOM ascl:1709.002.

The library needs to be initialized using a parameter file that is identical to
the one used for the actual  program, and provides a single
function that takes an array of positions, smoothing lengths and masses as an
input, and outputs an array of neutral hydrogen fractions. Under the hood, it
converts the positions, smoothing lengths and masses into a density field that
is mapped onto a density grid, and then uses sources from the parameter file to
illuminate this density field and compute self-consistent neutral fractions.
These neutral fractions are then mapped back to the original SPH positions using
the provided smoothing lengths.

It should be fairly straightforward to use the same approach to couple the code
to other types of hydrodynamical codes, like AMR codes.

Some parts of the code are also wrapped into a Python library. This library is
not meant to run full photoionization simulations like the C/C++/Fortran
counterpart, but instead can be used to facilitate the analysis of simulation
snapshots, by e.g. providing access to the line cooling data used by the code.




 §.§ Unstructured grid generation


An important feature of  is the option to use an unstructured
Voronoi grid as the main grid structure for both the radiation transfer and the
hydrodynamics. Due to the poor scaling properties of the incremental
construction algorithm used in <cit.>, we decided to
implement two new algorithms. The first algorithm (which for historical reasons
is called the “old” Voronoi algorithm) is our own rewritten version of the no
longer actively supported voro++
library[<http://math.lbl.gov/voro++/>]. This algorithm works in
most cases. However, we found that in some very specific degenerate cases, the
voro++ algorithm can produce the wrong Voronoi grid without crashing (as can be
confirmed graphically, or by computing the total volume of all cells). This is
due to the way the algorithm handles degeneracies.

Since these degenerate cases do in fact happen when starting a moving-mesh
hydrodynamical simulation from a regular initial grid, we also implemented a
new, more scalable version of the incremental construction algorithm used in
 (called the “new” algorithm). The current version of this grid
construction algorithm is about a factor 3 slower than the voro++ algorithm, but
is completely robust to any degeneracies due to the usage of arbitrary precision
arithmetics. We plan to release this algorithm as a standalone library that can
be used as a replacement for voro++ (Vandenbroucke et al., in
prep.).

In this work, we will use the old Voronoi algorithm whenever an unstructured
grid is used. This does not affect the results we show in any way, as the
Voronoi grid for a set of points is a unique geometrical structure, and is
independent of the way it is computed.




§ BENCHMARKS




In order to verify that our code produces physically accurate results in an
efficient way, we run a number of benchmark tests. Tests are grouped together
into four categories:

  
  *  tests that verify the ionization algorithm,
  
  *  tests that verify the combined ionization and temperature computation
  algorithm,
  
  *  tests that verify the radiation hydrodynamics algorithms, and
  
  *  tests that check the parallel efficiency of the code.


For the fourth category, we just reuse tests from the three other categories.
The initial conditions and analysis scripts for all benchmark tests are part of
the public version of the code.




 §.§ Strömgren sphere



To test the ionization algorithm, we run a simple test, inspired by the work of
<cit.>, in which a single ionizing source is at the centre of a
box containing a homogeneous density field consisting of only hydrogen. If we
assume that the source completely ionizes all hydrogen within a radius R_s
(the Strömgren radius), while the hydrogen outside this region stays
neutral, then the ionization balance equation for the ionized region is
<cit.>:

    Q = 4π/3 R_s^3 n_H^2 α_H^+ →H^0 (T),

where we made the assumption that n_H = n_H^+ within the
ionized region. If we also assume a fixed temperature T_0, then we get an
analytic expression for the constant Strömgren radius:

    R_s = ( 3Q/4πn_H^2α_H^+ →H^0(T_0))^1/3.


If the ionized region itself emits ionizing radiation (through the diffuse
field), the ionization balance equation changes to

    Q ( 1 + P_r,H(T) + P_r,H(T)^2 + ... ) 
    
    = 4π/3 R_s^3 n_H^2 α_H^+ →H^0 (T),

where P_r,H(T) is the reemission probability for ionizing radiation.

Using 1 + x + x^2 + ... = 1/1-x, (<ref>) now
becomes:

    R_s' = ( ( Q/1 - P_r,H(T_0))
    3/4πn_H^2α_H^+ →H^0(T_0))^1/3.


We will hence run two different versions of the test, that test different parts
of the algorithm:

  
  *  a version that does not include the diffuse reemission field and
  should reproduce (<ref>), and
  
  *  a version that includes diffuse radiation and should produce a larger
  ionization region, as given by (<ref>).


For both tests, we will use a cubic box of 10 × 10 × 10 pc
containing gas with a hydrogen number density of n_H =
100 cm^-3. At the centre of the box, we put a single source with a
luminosity of Q = 4.26 × 10^49 s^-1 with a monochromatic
spectrum that emits photons at the ionization threshold energy for hydrogen,
ν = 13.6 eV. We assume a constant photoionization cross section for
neutral hydrogen of σ_H^0 →H^+ = 6.3
× 10^-18 cm^-2, and a constant radiative recombination rate
α_H^+ →H^0 (T_0) = 4 ×
10^-13 cm^3 s^-1. The abundances, photoionization cross
sections and recombination rates for all other elements and ions are set to
zero.

We use a Cartesian density grid of 64 × 64 × 64 cells, and use
10^6 photon packets for 20 iterations to get a converged result.



  §.§.§ No diffuse field






This test corresponds to the benchmark test stromgren, and the setup is
as described above. The resulting hydrogen neutral fraction as a function of
radius is shown in the top panel of figure:stromgren_results. The
code accurately reproduces the expected Strömgren radius given by
(<ref>).

figure:stromgren_convergence shows the evolution of the hydrogen
neutral fraction profile with the number of iterations used. Initially, the
neutral fraction is set to a very low value everywhere in the box, so that the
ionizing radiation can illuminate a large region efficiently. After the first
iterations of the ionization state computation, the neutral fraction in the
outer regions quickly goes up until a converged result is reached. The result is
already well converged after 6 iterations.



  §.§.§ Diffuse field


This test corresponds to the benchmark test stromgren_diffuse, and
includes diffuse reemission with a reemission probability P_r,H(T_0) =
0.364 (corresponding to T_0 = 8000 K). The resulting hydrogen neutral
fraction profile is shown in the bottom panel of
figure:stromgren_results. As expected, the ionized region is larger
in this case. Our code still accurately reproduces the expected Strömgren
radius given by (<ref>).




 §.§ Lexington benchmarks



To test the combined temperature and ionization calculation for the full system
including metals, we run two of the benchmark tests that were the result of a
1995 workshop in Lexington and that are known as the Lexington benchmarks
<cit.>. The initial conditions for these tests can be found in
<cit.>, and correspond to the HII40 and HII20 model in that work.

The tests consist of a uniform density box with a hydrogen number density
n_H = 100 cm^-3, in which a central spherical region with
radius R_in = 3 × 10^16 m is evacuated. In the centre
of the evacuated region we place a single stellar source with a luminosity Q =
10^49 s^-1 for the low temperature benchmark, and Q = 4.26 ×
10^49 s^-1 for the high temperature benchmark, with a black body
spectrum. The temperature T_BB of the black body spectrum is also
different for the two tests, with T_BB = 20,000 K for the low
temperature benchmark, and T_BB = 40,000 K for the high
temperature version. The abundances (relative to hydrogen) of helium and the
various coolants are set to the values listed in
table:lexington_abundances.



To compare the test results, a number of quantities are computed:

  
  *  The total Hβ luminosity, which is computed from a power law
  fit to the data of <cit.>, using the electron density and
  temperature derived from the photoionization simulation.
  
  *  The height of the Balmer Jump Δ (BC 3645), defined as the
  jump in the Balmer continuum flux in a synthetic spectrum between the flux at
  3643 Å and 3681 Å. To obtain synthetic spectra, we use
  the continuum emission coefficients for hydrogen and helium from
  <cit.>. Note that <cit.> and other authors wrongly
  quote this value in units of Å, while the actual value is in
  Å^-1.
  
  *  The inner temperature T_inner at the boundary of the
  evacuated region.
  
  *  The average density weighted temperature <cit.>
  
    ⟨ T [ n_e n_H^+ ] ⟩ = ∫ n_e n_H^+
      T dV/∫ n_e n_H^+dV.

  
  *  The outer radius of the ionization region, defined as the average
  radius of cells with hydrogen neutral fractions in the range f_H^0∈ [0.1, 0.2].
  
  *  The ratio of the density weighted ionized fractions of hydrogen and
  helium, defined as
  
    ⟨ f_He^+⟩/⟨ f_H^+⟩
      = ∫ n_e f_He^+dV/∫ n_e n_H^+dV.



Apart from those, we also compute the line strengths of a number of emission
lines, relative to the total Hβ luminosity. These lines are a subset of
the emission lines that are used for the metal line cooling (see
<ref>), and their strength is computed in the same way
(summed over all cells).

For both tests, we set up a box of size 6 × 6 × 6 pc,
using a Cartesian density grid of 64 × 64 × 64 cells. We use
10^8 photons and 20 iterations to get a converged result for all coolants.

To set up the initial condition with a vacuum region in the centre, we use a
special implementation of the DensityFunction used to set up the density
field, called BlockSyntaxDensityFunction. This implementation uses a very
simple geometrical block description of the initial condition, which is the same
as used by the initial condition generator of 
<cit.>.



  §.§.§ Low temperature benchmark






This corresponds to the lexingtonHII20 benchmark test, and uses a black
body spectrum with T_BB = 20,000 K. The resulting ionic
fraction profiles for hydrogen, helium and several coolants are shown in the top
panel of figure:lexington_profiles, together with reference values
from the 1D code Cloudy. They follow the same trends as observed in
e.g. <cit.>, and overall agreement is pretty good. The resulting
temperature profile is shown in the top panel of
figure:lexington_temperatures. Our result follows the same overall
trend as the reference curve, but we systematically overestimate the temperature
in the ionized region. We do reproduce the correct peak temperature at the
ionization radius. This difference can be attributed to the different cooling
rates we use compared to Cloudy.

table:lexingthonHII20 lists the line strengths and comparison
quantities used by <cit.>, and compares them with the median
values given in that paper. All values have the expected order of magnitude,
although some values deviate significantly. This can be attributed to a
combination of our overall higher temperature, and different line emission
rates.





  §.§.§ High temperature benchmark


This corresponds to the lexingtonHII40 benchmark test, and uses a black
body spectrum with T_BB = 40,000 K. The resulting ionic
fraction profiles for hydrogen, helium and several coolants are shown in the
bottom panel of figure:lexington_profiles, and they again follow the
same trends as observed in <cit.>. The resulting temperature profile
is shown in the bottom panel of figure:lexington_temperatures. We
again notice an overall higher temperature in most of the ionized region
(although we actually underestimate the central temperature), which again is due
to the different data values used by our code.

table:lexingthonHII40 lists the line strengths and comparison
quantities. Note that we cannot compute the [C II] 1335 Å
line strength given in <cit.>, as this corresponds to a
transition outside the five lowest lying levels of C^+. Neither can we compute
the [O IV] 25.9 μm line strength, since O^++ has an
ionization potential that is (just) higher than the 54.4 eV upper limit
of our energy interval, and we hence cannot track O^+++. The values
generally agree with the <cit.> values, although there are again
significant differences, especially for the sulphur lines.






 §.§ STARBENCH benchmark



To test the coupling between the radiation transfer algorithm and our
hydrodynamical integration scheme, we use the benchmark D-type expansion of an
HII region which is part of the STARBENCH project <cit.>. The setup
is similar to the Strömgren test introduced above, but instead of assuming a
static solution, we allow the gas to react to the increased pressure due to the
higher temperature of the ionised region, and study the expansion of the
resulting ionized region over time.

We put a central source with a luminosity of 10^49 s^-1 and a
monochromatic spectrum that emits at 13.6 eV in a box of 2.512
× 2.512 × 2.512 pc containing only hydrogen with a density
of 3113 cm^-3. The hydrogen photoionization cross section is set to
6.3×10^-18 cm^2, while the hydrogen recombination rate is set
to 2.7×10^-13 cm^3 s^-1. We assume there is no
diffuse radiation field.

We assume a very simple isothermal equation of state (in practice this is
realized by setting the adiabatic index to γ = 1.0001), and assume that
the ionized region has a constant temperature T_i = 10,000 K, while the
neutral region has a constant temperature T_n = 100 K.

The system is evolved in time until t = 0.141 Myr, and we keep track of
the evolution of the ionization front.

There is no strict analytic solution for this problem, but there are two
reference solutions for the evolution of the ionization front as a function of
time. The first is the so called Spitzer solution <cit.>

    R_Sp(t) = R_s ( 1 + 7/4c_s,i t/R_s)^4/7,

where c_s,i is the (constant) sound speed in the ionized region, and R_s
is the Strömgren radius, as defined in (<ref>).

The second solution is due to <cit.>:

    R_HI(t) = R_s ( 1 + 7/4√(4/3)c_s,i
    t/R_s)^4/7,

and evolves at a somewhat faster rate. It is worth pointing out that we do not
require our simulation to reproduce any one of these solutions, but we do
require it to be close to them.

As a measure of the ionization front radius, we will use the average radius of
cells with neutral fractions in the range [0.8, 0.9]. Due to the sharp
transition from ionized to neutral (as can be seen from
figure:starbench_profile), using different boundaries for this
interval does not change the ionization radius much, as long as we make sure we
exclude noisy cells with x_H < 0.1 or x_H > 0.9.

We will run two versions of this test: a version that uses a static Cartesian
grid of 64 × 64 × 64 cells (the Eulerian solution), and a
version that uses a co-moving Voronoi mesh with 10,000 grid generator positions
sampled from a uniform distribution and regularized using Lloyd's algorithm for
10 iterations (the Lagrangian solution). For both, we apply the photoionisation
algorithm after every hydrodynamics step, using 10 iterations. The Eulerian
version uses 10^6 photon packets, while the Lagrangian version (with a lower
effective grid resolution) uses 10^5. These values were found to give a good
trade-off between accuracy and computational efficiency.



  §.§.§ Eulerian solution






This version corresponds to benchmark test starbench. It evolves the
system forward in time using 2,048 fixed size time steps. The evolution of the
ionization front is shown in figure:starbench_radius, and very
closely follows the Spitzer solution (<ref>). The left panels
in figure:starbench_profile show the density and neutral fraction as
a function of radius for t=0.0987 Myr.



  §.§.§ Lagrangian solution


This version corresponds to benchmark test star­bench_­vo­ro­noi. It
evolves the system forward in time using only 256 fixed size time steps (since
the average cell size is larger in this case). The evolution of the ionization
front is again shown in figure:starbench_radius. This time, the
ionization front first follows the Hosokawa-Inutsuka solution
(<ref>), and then slows down to the Spitzer solution
(<ref>). The density and neutral fraction profiles are shown in
the right panels of figure:starbench_profile.




 §.§ Parallel efficiency


As mentioned in <ref>, we have made most of our
algorithm inherently parallel by designing it in terms of small tasks. In this
subsection, we will illustrate how this affects the parallel scaling of the
code. Since different modes of the code use different parts of the algorithm, we
will rerun all of the benchmarks tests introduced above with different numbers
of shared memory threads and distributed memory processes.

For all runs, we use a single node of our local high performance computing
cluster Kennedy. This node has a 2.10 GHz Intel Xeon E5-2683 v4 processor with
16 physical cores, hyper-threaded to run 32 threads in parallel. Note that
because of the hyper-threading, we do not expect perfect scaling if more than 16
cores are used, since then different threads will be competing for resources.

We will quantify the scaling through the speed up S(n) of the algorithm
as a function of the number of computing units n. It is given by

    S(n) = t(1)/t(n),

where t(n) is the total runtime of the algorithm when using n computing
units. For a hypothetical code that scales perfectly, the speed up is simply
given by S_p(n)=n.

In practice, there are a number of important factors that can affect the
parallel scaling, and cause S(n) to lie below S_p(n):

  
  *  The existence of serial parts of the code that are only executed by
  one computing unit while the other computing units are idle, or that are
  executed by all computing units (and hence duplicate work). These are located
  in parts of the algorithm that are not parallel. If t_s(n) and t_p(n) are
  respectively the runtime of the serial and the parallel part of the code when
  using n computing units, then the theoretical maximum speed up is always
  lower than S_p(n) (assuming t_p(n) = t_p(1)/n):
  
    S_p,s(n) = ( t_s(1) + t_p(1)/nt_s(1) + t_p(1)) n.

  
  *  The overhead caused by running the algorithm in parallel. This
  overhead can be caused by extra code that needs to be executed as part of the
  parallelization strategy, or by delays caused by different computing units
  fighting over hardware access. Overhead will cause t_p,o(n) >
  t_p(1)/n, i.e. t_p,o(n) = t_p(1)/n + t_o(n) (with t_o(1)
  = 0).
  If the overhead is parallel, it is shared among the various computing
  units (t_o(n) = t_o,c/n). We can determine the constant from the 2
  computing units measurement, and the speed up including overhead is
  
    S_p,s+op(n) = ( t_s(1) + t_p(1)/n t_s(1) + 2 t_p,o(2))
      n.

  If on the other hand the overhead is serial, it is constant per computing unit
  (t_o(n) = t_o,u), and the speed up is given by
  
    S_p,s+os(n) = 
    ( t_s(1) + t_p(1)/n t_s(1) + n t_p,o(2) + ( 1 - n/2) t_p(1)) n.

  
  *  The occurrence of load imbalances between different computing units
  that cause some computing units to sit idle while they are waiting for other
  computing units to finish a task. While load imbalances affect the scaling in
  the same way as serial parts of the code, we cannot directly measure them, and
  treat them in the same way as the overhead.


For this version of , we made the decision to only parallelize
those parts of the algorithm that would lead to the most significant speed up,
i.e. the photon traversal algorithm and the cell based computations. This means
that there is still a significant fraction of serial code left in the algorithm,
which we will address in future versions of the code. This will limit the
expected scaling behaviour of the code, especially in simulations with small
grid sizes and small photon packet numbers.

Since the serial version of the code uses the same task based strategy as the
shared memory parallel version, the latter has no code based overhead (unless
the code is configured without OpenMP, in which case the OpenMP calls will be
absent). The only overhead will be caused by hardware concurrency issues: since
different threads might attempt to write to the same grid cell during the photon
traversal step, we need to lock each cell upon access. If another thread already
obtained access to the cell, the thread trying to acquire the lock will sit idle
until that cell becomes available, causing a small overhead.

The distributed memory parallel version of the code has a much larger code based
overhead, as it requires communication between the different parallel processes
which is completely absent from the serial version. The overhead is particularly
large in our code, as we made the decision to store variables per cell rather
than per variable type, which means we need to repack variables into separate
communication buffers before we can send them to another process. We plan to
address this issue as part of a future distributed memory parallel version that
includes a domain decomposition to distribute the grid across multiple
processes. We hence do not expect very good distributed memory scaling for our
current version.

The task based shared memory parallelization strategy we use automatically takes
care of the load balancing on a single node, as different threads get tasks from
a shared task pool. The time a single thread will need to wait cannot be longer
than the longest time it takes to finish a single task, so we can control the
former by adapting the size of the latter. When a JobMarket (see
<ref>) starts spawning tasks of a specific type,
it usually starts with reasonably large task sizes, and then gradually makes the
tasks smaller, until some lower limit is reached. This is illustrated in
figure:jobtimes. Using large initial task sizes limits the overhead
caused by calls to the JobMarket, while the lower size limit sets the
maximum load imbalance between different threads.

We do not perform any specific load balancing between distributed memory
parallel processes, and simply try to assign equal photon numbers and cell
numbers to each process when executing a step in parallel. A better load
balancing scheme will be part of the future domain decomposed distributed memory
algorithm.

Below, we will give the scaling results for the various benchmark tests. We will
focus on shared memory parallelization for all tests, and only show distributed
memory parallelization scaling for the Strömgren benchmark, as we still plan
to change our distributed memory parallelization strategy in future versions of
the code.



  §.§.§ Photoionization only




  
Shared memory scaling





table:stromgren_scaling_openmp and
figure:stromgren_scaling_openmp show the shared memory scaling
measurements for the default version of the Strömgren benchmark test
(without diffuse field). It is immediately obvious that there is a large
overhead in the parallel runs. Comparing the speed up with the two speed up
curves including overhead, we conclude that this overhead is parallel, and hence
shared among the threads.

The most likely cause of the constant overhead is the locking mechanism, since
the number of times a lock is used depends on the number of visited cells and is
hence constant for a given simulation. In the serial run, no locks need to be
set, and there is almost no overhead when locking and unlocking a cell. As soon
as more than 1 thread is used, locking is necessary and the constant overhead
will be added to the simulation time.

If we compare the scaling to S_p,s+op(n), the parallel overhead scaling
curve, we see good scaling up to 20 threads. After that, scaling decreases
significantly. This is expected, as the 32 threads have to compete for the 16
available physical cores.



  
Distributed memory scaling





table:stromgren_scaling_mpi and
figure:stromgren_scaling_mpi show the distributed memory scaling
measurements for the default version of the Strömgren benchmark test
(without diffuse field). Again there is a large overhead, but this time it is a
serial overhead that is almost constant per process. This overhead is due to the
repacking and communication, and due to load imbalances. Although there is a
speed up and hence some gain from using multiple processes, the parallel scaling
is poor.



  
Time line



figure:jobtimes shows a time line of the Strömgren benchmark
test, run on 4 shared memory threads on the same node. The different coloured
bars represent different tasks being executed, while the whitespace represents
serial parts of the code, or parts where the threads are effectively waiting
until all threads finished a specific type of jobs. Overall, the task based
parallelism works well to reduce load imbalances between different threads.
However, there are still some serial parts of the code that limit the
scalability. This graph also does not show the time threads spend waiting to
acquire locked resources, which can affect the summed total runtime of the
photon traversal task.



  §.§.§ Cooling and heating






table:lexington_scaling and figure:lexington_scaling show
the shared memory scaling measurements for the default high temperature
Lexington benchmark test. As in the Strömgren run, there is a considerable
overhead due to the locking mechanism, but since we do a lot more work, this
overhead is less noticeable. Overall, the code scales well up to 20 threads.



  §.§.§ Radiation hydrodynamics






table:starbench_scaling and figure:starbench_scaling show
the shared memory scaling measurements for a short (10 time step) version of the
default STARBENCH benchmark test. In this case, there is a significant serial
fraction, and a considerable amount of overhead as well. The overhead is
parallel at first, but switches to serial for high thread number. This is likely
due to hardware issues, as the hydrodynamical integration scheme is more
computation bound, and hence does not benefit from core hyper threading.
Overall, there is still a lot of room for improvement of the scalability.




§ CONCLUSION



We presented the public Monte Carlo photoionization and moving-mesh RHD code
, highlighted the implemented physics, and the most important
design considerations during code development. We illustrated the usage and
performance of the code with a number of relevant benchmark tests, and showed
that it produces accurate and reproducable scientific results.

This work accompanies the first official release of the code (
1.0), and describes the code as it is for this version. A number of important
improvements to the code are already on the way:

  
  *  implementation of a second order hydrodynamics solver,
  
  *  implementation of a distributed memory domain decomposition, using a
  scheme similar to <cit.>, which would enable us to overcome
  current memory limitations that make it impossible to run large grids,
  
  *  optimization and full implementation of an alternative Voronoi grid
  construction algorithm (Vandenbroucke et al., in prep.), and
  
  *  implementation of more extensive atomic data using the Stout
  database <cit.>.

These will be part of future code releases.




§ ACKNOWLEDGEMENTS



We want to thank the anonymous referee for a positive and constructive report
that helped clarify some key points in the paper.
We acknowledge support from STFC grant ST/M001296/1. We would like to thank Maya
Petkova for trying out the code and coupling it to PHANTOM, Daniel
Hanaway for suffering through the complications of trying to get the code to
work under Windows, and Diego Gonçalves and Nina Sartorio for useful
discussions about the coupling between radiation and hydrodynamics. Many thanks
as well to Pedro Gonnet and Matthieu Schaller for their constructive input about
task based parallelism; the latter also for advocating the use of a linear scale
in scaling plots.




§ REFERENCES




elsarticle-harv







§ FITTING PARAMETERS FOR METAL LINE COOLING VELOCITY-AVERAGED COLLISION
STRENGTHS


Below we list the fit parameters to the collision strength data for all elements
with low lying excited levels that are used in the code, grouped together per
element and per ion of that element. The figures show the data points that were
used in the fit, the actual fitting curve (with the same parameter precision as
displayed in the tables), and the relative difference between the data values
and the fitting curve, given by

    ξ_i → j (T) = | d(i → j, T) - f(i
    → j, T)|/d(i → j, T) + f(i → j, T),

where d(i → j, T) and f(i → j, T) are the data value
and fit at temperature T for the transition from level i to level j
respectively.





































