%----------------------------------------
% -begin arches-component.tex
%----------------------------------------
% 
The ARCHES component was designed for the simulation of turbulent reacting flows
with participating media radiation.  It is a three-dimensional, Large Eddy
Simulation (LES) code described in \cite{Smith}. ARCHES uses a low-Mach number
(Ma$<\!0.3$), variable density formulation to simulate heat, mass, and
momentum transport in reacting flows.

The LES algorithm solves the filtered, density-weighted, time-dependent
coupled conservation equations for mass, momentum, energy, and particle moment
equations in a Cartesian coordinate system \cite{Arches}.
%
This set of filtered equations is discretized in space and time and solved on
a staggered, finite volume mesh.  The staggering scheme consists of four
offset grids, one for storing scalar quantities and three for each component
of the velocity vector.  Stability preserving, second and third order explicit
time-stepping schemes are used to advance the simulation in time.  
%
For the spatial discretization of the LES scalar equations, flux limiting
schemes for the convection operator are used to ensure that scalar values
remain bounded.  For the momentum equation, a central differencing scheme for
the convection operator is used for energy conservation.  All diffusion terms
are computed with a second-order approximation of the gradient. Overall,
ARCHES is second-order accurate in space and time.  
%
The ARCHES code is massively parallel and highly scalable through 
its integration in the Uintah framework \cite{Rawat}, and also through use of
parallel solvers like Hypre \cite{Falgout} and PETSc \cite{petsc-user-ref}.
As part of the ARCHES development, substantial research has been done on
radiative heat transfer using the parallel Discrete Ordinates Method and the
P1 approximation to the radiative transport equation \cite{Krishnamoorthy3}.\\*

In reacting flow simulations, the main computational cost is the solution of the large number of
systems of linear equations required by the Discrete Ordinates Method. While the solution
of these systems can be made to scale\cite{John-UUSCI-2012-002}, it is important to reduce this cost.
With this cost reduction in mind, more recent work has
been based upon the use of more efficient Reverse Monte Carlo Ray Tracing (RMCRT)
methods, e.g.~\cite{Modest,Sun,hunsaker_2011}. RMCRT lends itself to scalable parallelism
because the intensities of each ray are mutually exclusive. Therefore, multiple
rays can be traced simultaneously at any given cell and time step. Creating an efficient,
GPU-accelerated software component based on RMCRT methods is the focus of 
the Computational Experiments section (\ref{sec:computational-experiments}).

\begin{figure}
    \includegraphics[width=0.45\textwidth,draft=false]{3p3diagram.png}
    \caption{Outline of Reverse Monte Carlo Ray Tracing}
    \label{fig:RMCRT}
\end{figure}

\subsection{Developing a Uintah Radiation Model}
\label{sec:radiation_integration}
In this study, we propose to extend the Uintah
framework so that problems involving radiation can \emph{also} be directly supported within
Uintah. Some kinds of radiation transport problems already use CFD codes
and AMR techniques~\cite{AMR-radiation,pernice}; However, other problems require the concept
of tracing rays or particles, such as the simulation of light transport, heat, radiation, or
electromagnetic waves.

The approach adopted in Uintah is on using the RMCRT methods, as described by \cite{Modest}.
This approach has the important advantage that
by using the principle of reciprocity in radiative transfer, rays are traced backwards 
from the computational cell thus eliminating the need to track ray bundles that never reach 
that cell \cite{Modest}. In RMCRT, rather than following a ray forward and calculating the energy it has
lost, the amount of incoming intensity from its path absorbed by the origin where the ray was emitted
is calculated. As Sun \cite{Sun} points out, RMCRT is more amenable to domain decomposition 
and thus parallel implementation due to the backward nature of the process.
Figure \ref{fig:RMCRT} shows the back path of a ray
from S to the emitter E, on a nine cell structured mesh patch. Each $i^\mathrm{th}$ cell has its own temperature
$T_i$, absorption coefficient $\kappa_i$,  scattering coefficient $\sigma_i$ and appropriate pathlengths
$l_{i,j}$.  
In each case the incoming intensity is calculated, say in cell 4, and then
traced back through the other cells. The intensity is integrated along the ray path to compute 
a divergence of the heat flux or a surface flux.
When a ray hits a boundary (as on surface $17$ in the figure), it can be either reflected or
absorbed depending on the surface properties. Rays are terminated when their intensity is sufficiently small.

Despite the improved efficiency over forward MCRT, there are considerable
challenges in the efficient implementation of RMCRT as it is an all-to-all method, where all of the
geometry information and property model information for the entire computational domain must
be present on each processor \cite{Sun}. This nature severely limits the
size of the problem that can be computed due to memory constraints, especially with large highly resolved physical domains.
This challenge is being addressed by using the multi-level mechanisms within Uintah
to represent a portion of the domain at a coarser resolution, thus lowering the memory usage \cite{hunsaker_2011}.
The hybrid memory approach of Uintah also helps as only one copy of geometry is needed per multi-core node. 
In general, the data required by the RMCRT algorithm is projected
to all of the coarser levels, with each level spanning the entire domain.
For each fine level patch, data from the coarser levels is retrieved from the Uintah \emph{data warehouse} so
it encompasses the patch in a stair step fashion. 

CPU-only scalability studies of the RMCRT for the benchmark problem as described by Burns and Christon ~\cite{benchmark},
were run on on a single level\cite{hunsaker_2011} with $256^3$ cells, using 25 \& 100 rays per cell.  
Each scaling run was run for 10 timesteps, 1 patch per processor, and the mean time per timestep 
was computed. These preliminary results show reasonable scaling up to 768 cores, above this
the loss of scalability is perhaps due to increased communication costs and/or a load imbalance.
Nevertheless these results provide a good proof of concept and an excellent starting point for this work.
%
%----------------------------------------
% -end arches-component.tex
%----------------------------------------