Radiative transfer

RtSolver

Radiative transfer solvers are supported in the form of

Diffusion solvers, i.e. first order Eddington solvers, which take the form of a Helmholtz equation.
Particle solvers, which track photons are particles (e.g., Monte Carlo sampled solvers).

The solvers share a parent class RtSolver, and code that uses only the RtSolver interface will should be able to switch between the two implementations. Note, however, that the radiative transfer equation is inherently deterministic while Monte Carlo photon transport is inherently stochastic. The diffusion approximation relies on solving an elliptic equation in the stationary case and a parabolic equation in the time-dependent case, while the Monte-Carlo approach solves solves for fully transient or instantaneous transport.

Tip

The source code for the solver is located in $DISCHARGE_HOME/Source/RadiativeTransfer and it is a fairly lightweight abstract class. As with other solvers, RtSolver can use a specified Realm.

To use the RtSolver interface the user must cast from one of the inherited classes (see Diffusion approximation or Monte Carlo solver). Since most of the RtSolver is an interface which is implemented by other radiative transfer solvers, documentation of boundary conditions, kernels and so on are found in the implementation classes.

RtSpecies

The class RtSpecies is an abstract base class for parsing necessary information into radiative transfer solvers. When creating a radiative transfer solver one will need to pass in a reference to an RtSpecies instantiation such that the solvers can look up the required information. Currently, RtSpecies is a lightweight class where the user needs to implement the function

/*!
  @brief Get kappa (i.e. the inverse absorption length) at physical coordinates. 
  @param[in] a_pos Physical coordinates. 
*/
virtual Real
getAbsorptionCoefficient(const RealVect a_pos) const = 0;

This absorption coefficient is used in both the diffusion (see Diffusion approximation) and Monte Carlo (see Monte Carlo solver) solvers.

Important

Upon construction, one must set the class member m_name, which is the name passed to the actual solver.

Setting the source term

RtSolver stores a source term $\eta$ on the mesh, which describes the number of photons that are generated produced per unit volume and time. This variable can be set through the following functions:

/*!
  @brief Set source term
  @param[in] a_source Source term
*/
virtual void
setSource(const EBAMRCellData& a_source);

/*!
  @brief Set source
  @param[in] a_source Source term
*/
virtual void
setSource(const Real a_source);

/*!
  @brief Set source
  @param[in] a_source Source term (varies in space)
*/
virtual void
setSource(const std::function<Real(const RealVect a_pos)> a_source);

The usage of $\eta$ varies between the different solvers. It is possible, for example, to generate computational photons (particles) using $\eta$ when using Monte Carlo sampling, but this is not a requirement.

Diffusion approximation

EddingtonSP1

The first-order diffusion approximation to the radiative transfer equation is encapsulated by the EddingtonSP1 class which implements a first order Eddington approximation of the radiative transfer equation. EddingtonSP1 implements RtSolver using both stationary and transient advance methods (e.g. for stationary or time-dependent radiative transport). The source code is located in $DISCHARGE_HOME/RadiativeTransfer.

Equation(s) of motion

In the diffusion approximation, the radiative transport equation is

(10)\[\partial_t\Psi + \kappa\Psi - \nabla\cdot\left(\frac{1}{3\kappa}\nabla\Psi\right) = \frac{\eta}{c},\]

where $\Psi$ is the radiative intensity (i.e., photons absorbed per unit volume`. Here, $\kappa$ is the absorption coefficient (i.e., inverse absorption length). This value can be spatially dependent, and is passed in through the RtSpecies function getAbsorptionCoefficient that was discussed above. Note that in the context below, $\kappa$ is not the volume fraction of a grid cell but the absorption coefficient. The above equation is called the Eddington approximation, with the closure relation being that the radiative flux is given by $F = -\frac{c}{3\kappa}\nabla \Psi$.

In the stationary case this reduces to a Helmholtz equation

(11)\[\kappa\Psi - \nabla\cdot\left(\frac{1}{3\kappa}\nabla\Psi\right) = \frac{\eta}{c},\]

Implementation

EddingtonSP1 uses multigrid methods for solving Eq. 10 and Eq. 11, see Linear solvers. To advance the solution, one will call the member function

/*!
  @brief Advance RTE onto state a_phi
  @param[in]    a_dt      Time step
  @param[inout] a_phi     RTE solution
  @param[in]    a_source  Source term
  @param[in]    a_zeroPhi Set phi to zero in initial guess for multigrid solve
  @note If you're not doing a stationary solve, this does a backward Euler solve. 
*/
virtual bool
advance(const Real a_dt, EBAMRCellData& a_phi, const EBAMRCellData& a_source, const bool a_zeroPhi = false) override;

Internally, this version will perform one of the following:

Solve Eq. 10 if using a transient solver. This is done using a backward Euler solve:

\[\left(1+ \kappa \Delta t\right)\Psi^{k+1} - \Delta t \nabla\cdot\left(\frac{1}{3\kappa}\nabla\Psi^{k+1}\right) = \Psi^{k} + \frac{\Delta t\eta^{k+1}}{c},\]

This equation is a Helmholtz equation for $\Psi^{k+1}$ which is solved using geometric multigrid, see Linear solvers.
Solve Eq. 11 if using instantaneous photon transport. This is done directly with a geometric multigrid solver, see Linear solvers.

Boundary conditions

Simplified domain boundary conditions

It is possible to set the following simplified boundary conditions on domain faces and embedded boundaries: All of these boundary condition specifications take the form <type> <value>.

Dirichlet, with a fixed value of $\Phi$. E.g., dirichlet 0.0.
Neumann, using a fixed value of $\partial_n\Phi$. E.g., neumann 0.0.
A Larsen-type boundary condition, which is an absorbing boundary condition in the form

\[\kappa\partial_n\Psi + \frac{3\kappa^2}{2}\frac{1-3r_2}{1-2r_1}\Psi = g,\]

where $r_1$ and $r_2$ are reflection coefficients and $g$ is a surface source, see [Larsen et al., 2002]. Note that when the user specifies the boundary condition value (e.g. by setting the BC function), he is setting the surface sourge $g$. In the majority of cases, however, we will have $r_1 = r_2 = g = 0$ and the BC becomes

\[\partial_n\Psi + \frac{3\kappa}{2}\Psi = g.\]

The user must then pass a value larsen <value>, where the value corresponds to the souce term $g$. Typically, this term is zero.

Tip

For radiative transfer, the Larsen boundary condition is usually the correct one as it approximately describes outflow of photons on the boundary. In this case the correct boundary condition is larsen 0.0.

Custom domain boundary conditions

It is possible to use more complex boundary conditions by passing in dirichlet_custom, neumann_custom, or larsen_custom options through the solver configuration options (see Solver configuration). In this case the EddingtonSP1 solver will use a specified function at the domain edge/face which can vary spatially (and with time). To specify that function, EddingtonSP1 has a member function

/*!
  @brief Set the boundary condition function on a domain side
  @param[in] a_dir      Coordinate direction.
  @param[in] a_side     Side (low/high)
  @param[in] a_function Boundary condition function.
  @details This sets a boundary condition for a particular domain side. The user must also specify how to use this BC in the input script. 
*/
virtual void
setDomainSideBcFunction(const int                               a_dir,
                        const Side::LoHiSide                    a_side,
                        const EddingtonSP1DomainBc::BcFunction& a_function);

Here, the a_function argument is simply an alias:

/*!
  @brief Function which maps f(R^3,t) : R. Used for setting the associated value and boundary condition type. 
*/
using BcFunction = std::function<Real(const RealVect a_position, const Real a_time)>;

Note that the boundary condition type is still Dirichlet, Neumann, or Larsen (depending on whether or not dirichlet_custom, neumann_custom, or larsen_custom was passed in). For example, to set the boundary condition on the left $x$ face in the domain, one can create a EddingtonSP1DomainBc::BcFunction object as follows:

// Assume this has been instantiated.
RefCountedPtr<EddingtonSP1> eddingtonSolver;

// Make a lambda which we can bind to std::function.
auto myValue = [](const RealVect a_pos, const Real a_time) -> Real {
   return a_pos[0] * a_time;
}

// Set the domain bc function in the solver.
eddingtonSolver.setDomainSideBcFunction(0, Side::Lo, myValue);

Warning

A run-time error will occur if the user specifies one of the custom boundary conditions but does not actually set the function.

Embedded boundaries

On the EB, we currently only support constant-value boundary conditions. In the input script, the user can specify

dirichlet <value> For setting a constant Dirichlet boundary condition everywhere.
neumann <value> For setting a constant Neumann boundary condition everywhere.
larsen <value> For setting a constant Larsen boundary condition everywhere.

The specification of these boundary conditions occurs in precise analogy with the domain boundary conditions, and are therefore not discussed further here.

Solver configuration

The EddingtonSP1 implementation has a number of configurable options for running the solver, and these are given below:

Listing 29 EddingtonSP1 solver configuration options. Run-time configurable options are highlighted.

# ====================================================================================================
# EddingtonSP1 class options
# ====================================================================================================
EddingtonSP1.verbosity             = -1         ## Solver verbosity
EddingtonSP1.stationary            = true       ## Stationary solver
EddingtonSP1.reflectivity          = 0.         ## Reflectivity
EddingtonSP1.kappa_scale           = true       ## Kappa scale source or not (depends on algorithm)
EddingtonSP1.plt_vars              = phi src    ## Plot variables. Available are 'phi' and 'src'
EddingtonSP1.use_regrid_slopes     = true       ## Slopes on/off when regridding

EddingtonSP1.ebbc        = larsen 0.0    ## Bc on embedded boundaries
EddingtonSP1.bc.x.lo     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'
EddingtonSP1.bc.x.hi     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'
EddingtonSP1.bc.y.lo     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'
EddingtonSP1.bc.y.hi     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'
EddingtonSP1.bc.z.lo     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'
EddingtonSP1.bc.z.hi     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'
EddingtonSP1.bc.z.hi     = larsen 0.0    ## Bc on domain side. 'dirichlet', 'neuman', or 'larsen'

EddingtonSP1.gmg_verbosity         = -1           ## GMG verbosity
EddingtonSP1.gmg_pre_smooth        = 8            ## Number of relaxations in downsweep
EddingtonSP1.gmg_post_smooth       = 8            ## Number of relaxations in upsweep
EddingtonSP1.gmg_bott_smooth       = 8            ## NUmber of relaxations before dropping to bottom solver
EddingtonSP1.gmg_min_iter          = 5            ## Minimum number of iterations
EddingtonSP1.gmg_max_iter          = 32           ## Maximum number of iterations
EddingtonSP1.gmg_exit_tol          = 1.E-6        ## Residue tolerance
EddingtonSP1.gmg_exit_hang         = 0.2          ## Solver hang
EddingtonSP1.gmg_min_cells         = 16           ## Bottom drop
EddingtonSP1.gmg_bottom_solver     = bicgstab     ## Bottom solver type. Either 'simple <number>' and 'bicgstab'
EddingtonSP1.gmg_cycle             = vcycle       ## Cycle type. Only 'vcycle' supported for now
EddingtonSP1.gmg_ebbc_weight       = 1            ## EBBC weight (only for Dirichlet)
EddingtonSP1.gmg_ebbc_order        = 2            ## EBBC order (only for Dirichlet)
EddingtonSP1.gmg_smoother          = red_black    ## Relaxation type. 'jacobi', 'red_black', or 'multi_color'

The multigrid options are analogous to the multigrid options for FieldSolverGMG, see Tuning multigrid performance.

Monte Carlo solver

McPhoto defines a class which can solve radiative transfer problems using discrete photons. The class derives from RtSolver and can thus be used also be used by applications that only require the RtSolver interface. McPhoto can provide a rather complex interaction with boundaries, such as computing the intersection between a photon path and a geometry, and thus capture shadows (which EddingtonSP1 can not).

The Monte Carlo sampling is a particle-based radiative transfer solver, and particle-mesh operations (see Particle-mesh) are thus required in order to deposit the photons on a mesh if one wants to compute mesh-based absorption profiles.

Tip

The McPhoto class is defined in $DISCHARGE_HOME/Source/RadiativeTransfer/CD_McPhoto.H. See the McPhoto C++ API for further details.

The solver has multiple data holders for systemizing photons, which is especially useful during transport kernels where some of the photons might strike a boundary:

In-flight photons.
Bulk-absorbed photons, i.e., photons that were absorbed on the mesh.
EB-absorbed photons, i.e., photons that struck the EB during a transport step.
Domain-absorbed photons, i.e., photons that struck the domain edge/face during a transport step.
Source photons, for letting the user pass in externally generated photons into the solver.

Various functions are in place for obtaining these particles:

/*!
  @brief Get m_photons
  @return m_photons
*/
virtual ParticleContainer<Photon>&
getPhotons();

/*!
  @brief Get bulk photons, i.e. photons absorbed on the mesh
  @return m_bulkPhotons
*/
virtual ParticleContainer<Photon>&
getBulkPhotons();

/*!
  @brief Get eb Photons, i.e. photons absorbed on the EB
  @return m_ebPhotons
*/
virtual ParticleContainer<Photon>&
getEbPhotons();

/*!
  @brief Get domain photons, i.e. photons absorbed on domain edges/faces
  @return m_domainPhotonsl
*/
virtual ParticleContainer<Photon>&
getDomainPhotons();

/*!
  @brief Get source photons
  @return m_sourcePhotons
*/
virtual ParticleContainer<Photon>&
getSourcePhotons();

Photon particle

The Photon particle is a simple encapsulation of a computational photon which is used by McPhoto. It derives from GenericParticle<2,1> and stores (in addition to the particle position):

The particle weight.
The particle mean absorption coefficient.
The particle velocity/direction.

Tip

The Photon class is defined in $DISCHARGE_HOME/Source/RadiativeTransfer/CD_Photon.H

When defining the McPhoto class, the particle’s absorption coefficient can be computed from the implementation of the absorption function method in RtSpecies.

Generating photons

There are several ways users can generate computational photons that are to be transported by the solver.

Fetch the source photons by calling

/*!
  @brief Get source photons
  @return m_sourcePhotons
*/
virtual ParticleContainer<Photon>&
getSourcePhotons();

The source photons can then be filled and added to the other photons.

Add photons directly, by first obtaining the in-flight photons through

/*!
  @brief Get m_photons
  @return m_photons
*/
virtual ParticleContainer<Photon>&
getPhotons();

Photons can then be added directly.

If the source term $\eta$ has been filled, the user can call McPhoto::advance to have the solver generate the computational photons and then advance them. This is the correct approach for, e.g., applications that always use mesh-based photon source terms and want to have the computational photons be generated on the fly.

Warning

The advance function is only meant to be used together with a mesh-based source term that the user has filled prior to calling the method.

When using the advance, the number of photons that are generated are limit to a user-specified number (see Solver configuration for further details).

Transport modes

McPhoto can be run as a fully transient, in which photons are tracked in time, or as an instantaneous solver. For the instantaneous mode, photon absorption positions are stochastically sampled with Monte Carlo procedure and the photons are immediately absorbed on the mesh. For the transient mode the photon advancement occurs over $\Delta t$, so there is a limited distance ($c \Delta t$) that the photons can propagate. In this case, only some of the photons will be absorbed on the mesh whereas the rest may continue their propagation.

Instantaneous transport

When using instantaneous transport, any photon generated in a time step is immediately absorbed on the boundary through the following steps:

Optionally, have the solver generate photons to be transport (or add them externally).
Draw a propagation distance $r$ by drawing random numbers from an exponential distribution $p(r) = \kappa \exp\left(-\kappa r\right)$. Here, $\kappa$ is computed by calling the underlying RtSpecies absorption function. The absorbed position of the photon is set to $\mathbf{x} = \mathbf{x}_0 + r\mathbf{n}$.

Warning

In instantaneous mode photons might travel infinitely long, i.e. there is no guarantee that $c\Delta t \leq r$.
Deposit the photons on the mesh.

Transient transport

The transient Monte Carlo method is almost identical to the stationary method, except that it does not deposit all generated photons on the mesh but tracks them through time. For each photon, do the following:

Compute an absorption length $r$ by sampling the absorption function at the current photon position.
Each photon is advanced over the time step $\Delta t$ such that the position is

\[\mathbf{x} = \mathbf{x}_0 + \mathbf{c}\Delta t.\]
Check if $\left|\mathbf{x}-\mathbf{x}_0\right| < r$ and if it is, absorb the photon on the mesh.

Other transport kernels

In addition to the above two methods, the solver interface permits users to add e.g. source photons externally and add them to the solvers’ transport kernel.

Solver configuration

McPhoto can be configured through its input options, see below:

# ====================================================================================================
# McPhoto class options
# ====================================================================================================
McPhoto.verbosity                = -1               ## Solver verbosity
McPhoto.instantaneous            = true             ## Instantaneous transport or not
McPhoto.max_photons_per_cell     = 32               ## Maximum no. generated in a cell (<= 0 yields physical photons)
McPhoto.num_sampling_packets     = 1                ## Number of sub-sampling packets for max_photons_per_cell. Only for instantaneous=true
McPhoto.blend_conservation       = false            ## Switch for blending with the nonconservative divergence
McPhoto.transparent_eb           = false            ## Turn on/off transparent boundaries. Only for instantaneous=true
McPhoto.plt_vars                 = phi src phot     ## Available are 'phi' and 'src', 'phot', 'eb_phot', 'dom_phot', 'bulk_phot', 'src_phot'
McPhoto.intersection_alg         = bisection        ## EB intersection algorithm. Supported are: 'raycast' 'bisection'
McPhoto.bisect_step              = 1.E-4            ## Bisection step length for intersection tests
McPhoto.bc_x_low                 = outflow          ## Boundary condition. 'outflow', 'symmetry', or 'wall'
McPhoto.bc_x_high                = outflow          ## Boundary condition
McPhoto.bc_y_low                 = outflow          ## Boundary condition
McPhoto.bc_y_high                = outflow          ## Boundary condition
McPhoto.bc_z_low                 = outflow          ## Boundary condition
McPhoto.bc_z_high                = outflow          ## Boundary condition
McPhoto.photon_generation        = deterministic    ## Volumetric source term. 'deterministic' or 'stochastic'
McPhoto.source_type              = number           ## 'number'      -> Source term contains the number of photons produced
                                                    # 'volume'      -> Source terms contains the number of photons produced per unit volume
                                                    # 'volume_rate' -> Source terms contains the volumetric rate
                                                    # 'rate'        -> Source terms contains the rate
McPhoto.deposition               = cic              ## 'ngp'  -> nearest grid point, 'cic' -> cloud-in-cell
McPhoto.deposition_cf            = transition       ## 'interp', 'halo', 'halo_ngp', 'transition'.

Tip

The McPhoto class includes a hidden input parameter McPhoto.dirty_sampling = true/false which enables a cheaper sampling method for discrete photons when calling the advance method. The caveat is that the method does not incorporate boundary intersect, only works for instantaneous propagation, and avoids filling the data holders that are necessary for load balancing.

Clarifications

When computational photons are generated through the solver, users might have filled the source term differently depending on the application. For example, users might have filled the source term with the number of photons generated per unit volume and time, or the physical number of photons to be generated. The two input options McPhoto.photon_generation and McPhoto.source_type contain the necessary specifications for ensuring that the user-filled source term can be translated properly for ensuring that the correct number of physical photons are generated. Firstly, McPhoto.source_type contains the specification of what the source term contains:

number if the source term contains the physical number of photons.
volume if the source terms contains the physical number of photons generated per unit volume.
volume_rate if the source terms contains the physical number of photons generated per unit volume and time.
rate if the source terms contains the physical number of photons generated per unit time.

When McPhoto calculates the number of physical photons in a cell, it will automatically determine from McPhoto.source_type, $\Delta V$ and $\Delta t$ how many physical photons are to be generated in each grid cell.

McPhoto.photon_generation permits the user to turn on/off Poisson sampling when determining how many photons will be generated. If this is set to stochastic, the solver will first compute the number of physical photons $\overline{N}_\gamma^{\text{phys}}$ following the procedure above, and then run a Poisson sampling such that the final number of physical photons is

\[N_{\gamma}^{\text{phys}} = P\left(\overline{N}_{\gamma}^{\text{phys}}\right).\]

Otherwise, if McPhoto.photon_generation is set to deterministic then the solver will generate

\[N_{\gamma}^{\text{phys}} = \overline{N}_{\gamma}^{\text{phys}}\]

photons. Again, these elements are important because users might have chosen to perform the Poisson sampling outside of McPhoto.

Important

All of the above procedures are done per-cell.

Example application

Example applications that use RtSolver are found in: