# Benchmarks for infinite medium, time dependent transport problems with
isotropic scattering

William Bennett and Ryan G. McClarren  Department of Aerospace and Mechanical
Engineering, University of Notre Dame, Notre Dame, Indiana, United States

(August 27, 2024)

###### Abstract

The widely used AZURV1 transport benchmarks package provides a suite of
solutions to isotropic scattering transport problems with a variety of initial
conditions (Ganapol et al. 2001). Most of these solutions have an initial
condition that is a Dirac delta function in space; as a result these
benchmarks are challenging problems to use for verification tests in computer
codes. Nevertheless, approximating a delta function in simulation often leads
to low orders of convergence and the inability to test the convergence of
high-order numerical methods. While there are examples in the literature of
integration of these solutions as Green’s functions for the transport operator
to produce results for more easily simulated sources, they are limited in
scope and briefly explained. For a sampling of initial conditions and sources,
we present solutions for the uncollided and collided scalar flux to facilitate
accurate testing of source treatment in numerical solvers. The solution for
the uncollided scalar flux is found in analytic form for some sources. Since
integrating the Green’s functions is often nontrivial, discussion of
integration difficulty and workarounds to find convergent integrals is
included. Additionally, our uncollided solutions can be used as source terms
in verification studies, in a similar way to the method of manufactured
solutions.

##  1 Introduction

The AZURV1 benchmark suite, developed by Ganapol et al. (2001), is an
indispensable verification tool in the transport community. Some of the works
that rely on these benchmarks include (Variansyah and McClarren 2022; Harel,
Burov, and Heizler 2021; Peng and McClarren 2021; Heizler 2010; Garrett and
Hauck 2013; Seibold and Frank 2014; Schlachter and Schneider 2018; Hauck and
Heningburg 2019; Heningburg and Hauck 2020). The AZURV1 benchmark extends the
work of Monin (1956) and contains solutions for delta-function initial
conditions in planar, line, point and spherical shell shapes. These solutions
can be considered Green’s functions111There is contention whether one should
use “Green’s function” or “Green function”. We follow the convention
recommended by Wright (2006) in retaining the possessive. for the respective
geometries.

Since running one of the AZURV1 problems with a numerical code requires the
difficult approximation of a delta function, Garrett and Hauck (2013) present
a method for integrating Ganapol’s line source problem to find the exact
solution for a Gaussian initial condition that resembles the line source but
is more manageable for a numerical solver. This nearby problem is cited by
(Seibold and Frank 2014; Schlachter and Schneider 2018; Hauck and Heningburg
2019; Heningburg and Hauck 2020). The popularity of this solution confirms
that there is interest for such solutions in the transport community. In this
work we present transport solutions for a variety of initial conditions and
sources to address this need. These solutions are considerably easier to
handle than Green’s functions for numerical codes and are, therefore, more
useful for convergence studies. Additionally, we give the solution for the
uncollided scalar flux for these problems. The uncollided flux can be used as
the source in a computer code that solves that transport equation. This method
is similar to the Method of Manufactured Solutions (MMS) (Salari and Knupp
2000; McClarren and Lowrie 2008), where a known solution is used to solve for
a source term (that is typically a complicated function of space, angle, and
time), which is then given to the numerical algorithm. Other researchers can
use the uncollided solutions we present as the prescribed source in a
verification test.

The remainder of this paper is organized as follows. In Section 2 we present
the multiple collision approach and the original plane pulse solution from
(Ganapol et al. 2001). We then integrate this solution over square and
Gaussian spatial distributions for an initial pulse and for a source that is
on for a fixed time, $t_{0}$ in Section 3. The pulsed line source solution is
presented in Section 5, as well as the integral over this pulse in a Gaussian
configuration.

##  2 Uncollided-Collided Split Transport Model

We begin with the neutral particle transport problem in an infinite medium
with isotropic scattering (Ganapol et al. 2001)

$\left(\frac{\partial}{\partial t}+\mu\frac{\partial}{\partial
x}+1\right)\psi(x,t,\mu)=\frac{c}{2}\,\phi(x,t)+\frac{1}{2}\,S(x,t),$ (1)

where the angular flux is represented by $\psi(x,t,\mu)$, the scalar flux by
$\phi(x,t)=\int^{1}_{-1}\\!d\mu^{\prime}\,\psi(x,t,\mu^{\prime})$, $S$ is a
source, and $c$ is the scattering ratio. The spatial coordinate $x$ is
measured in units of mean-free path and with a particle speed of unity, $t$
measures time in units of mean-free time; $\mu\in[-1,1]$ is the cosine of the
angle between a direction of travel and the $x$-axis. Because we are in an
infinite medium, we do not specify boundary conditions. We do, however, assert
that the initial conditions are zero-flux, unless otherwise specified.

We can split Eq. (1) into an equation for the angular flux of particles that
have not undergone a collision and an equation for the particles that have
undergone a collision. To do this we write
$\psi(x,t,\mu)=\psi_{\mathrm{u}}(x,t,\mu)+\psi_{\mathrm{c}}(x,t,\mu)$.

The equation for the uncollided angular flux, $\psi_{\mathrm{u}}$ is

$\left(\frac{\partial}{\partial t}+\mu\frac{\partial}{\partial
x}+1\right)\psi_{\mathrm{u}}(x,t,\mu)=\frac{1}{2}\,S(x,t).$ (2)

Notice that this equation is a purely absorbing transport equation. This fact
will allow us to write closed-form solutions for the uncollided scalar flux in
many instances.

The equation for the collided angular flux, $\psi_{\mathrm{c}}$, looks like
the original transport equation, Eq. (1), with the addition of a source term
from the uncollided solution:

$\displaystyle\left(\frac{\partial}{\partial t}+\mu\frac{\partial}{\partial
x}+1\right)\psi_{\mathrm{c}}(x,t,\mu)$
$\displaystyle=\frac{c}{2}\int^{1}_{-1}\\!d\mu^{\prime}\left(\psi_{\mathrm{c}}(x,t,\mu^{\prime})+\psi_{\mathrm{u}}(x,t,\mu^{\prime})\right)$
$\displaystyle=\frac{c}{2}\int^{1}_{-1}\\!d\mu^{\prime}\psi_{\mathrm{c}}(x,t,\mu^{\prime})+S_{\mathrm{u}}(x,t).$
(3)

Here we have defined a source term for the collided equation as

$S_{\mathrm{u}}(x,t)=\frac{c}{2}\int^{1}_{-1}\\!d\mu^{\prime}\psi_{\mathrm{u}}(x,t,\mu^{\prime})\\\
=\frac{c}{2}\phi_{\mathrm{u}}(x,t).$ (4)

A few points are in order regarding the uncollided and collided transport
equations, Eqs. (2) and (2). Firstly, it is clear that adding these two
equations together yields the original transport equation, Eq. (1).
Furthermore, we note that if the uncollided scalar flux, $\phi_{\mathrm{u}}$,
is known then one can solve a standard transport equation with source given by
$c\phi_{\mathrm{u}}/2$ to get the collided solution $\psi_{\mathrm{c}}$. This
means that one could use the uncollided solutions we give below as a source
term in a verification test because the solution will be an approximation to
$\phi_{\mathrm{c}}$.

Numerical methods that employ decomposition based on collisions are often more
efficient as they are able to allocate computational resources separately to
solve for the typically less smooth uncollided fluxes, as in (Alcouffe,
O’Dell, and Brinkley 1990; Hauck and McClarren 2013; Walters and Haghighat
2017). Our proposed source treatment is similar to these methods, except
instead of solving for the uncollided flux with a more refined mesh or more
angular degrees of freedom, an analytic or semi-analytic solution for the
uncollided scalar flux is inserted into Eq. (2), acting like a source term.

##  3 Planar pulse Green’s functions

For a one dimensional infinite plane pulsed source with isotropic scattering
the corresponding source is $S(x,t)=\delta(x)\delta(t)$. The solution for the
uncollided scalar flux is (Ganapol et al. 2001)

$\phi_{\mathrm{u}}^{\mathrm{pl}}(x,t)=\frac{e^{-t}}{2t}\Theta(1-|\eta|),$ (5)

with

$\eta\equiv\frac{x}{t},$ (6)

and $\Theta$ is a step function that returns unity for positive arguments and
zero otherwise. The solution for the collided flux is,

$\phi_{\mathrm{c}}^{\mathrm{pl}}(x,t)=c\left(\frac{e^{-t}}{8\pi}\left(1-\eta^{2}\right)\int^{\pi}_{0}\\!du\,\mathrm{sec}^{2}\left(\frac{u}{2}\right)\mathrm{Re}\left[\xi^{2}e^{\frac{ct}{2}(1-\eta^{2})\xi}\right]\right)\Theta(1-|\eta|),$
(7)

where the complex valued function $\xi$ is ,

$\xi(u,\eta)=\frac{\log q+i\,u}{\eta+i\,\mathrm{tan}(\frac{u}{2})}$ (8)

and,

$q=\frac{1+\eta}{1-\eta}.$ (9)

Equations (5) and (7) will be treated as Green’s function kernels and
integrated to find solutions for a variety of sources and initial conditions.
The general form for this integration is:

$\phi_{\mathrm{j}}(x,t)=\int^{t}_{0}\\!d\tau\,\int_{-\infty}^{\infty}\\!ds\,S(s,\tau)\,\phi_{j}^{\mathrm{pl}}(x-s,t-\tau),$
(10)

here $S$ is an arbitrary source, and the subscript $j=$ u or c denotes the
uncollided or collided scalar flux.

The total scalar flux ($\phi_{\mathrm{u}}+\phi_{\mathrm{c}}$) and the
uncollided scalar flux for this problem are shown in Figure 1. Because the
uncollided solution decays exponentially in time and is effectively zero on
the scale of the plots, we do not show the uncollided solution in the later
time panels. We also point out that at $t=1$ the presence of the wavefront is
noticeable at $x=1$; this feature also decays exponentially and is hardly
noticeable in the $t=5$ panel at $x=5$. This wavefront is one of the features
in the AZURV1 solution that can make it challenging to use this benchmark in
convergence studies for high-order numerical methods.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 1: Plane pulse scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{pl}}+\phi_{\mathrm{c}}^{\mathrm{pl}}$, for $c=1$
at several times; panel (a) also contains the uncollided scalar flux,
$\phi_{\mathrm{u}}^{\mathrm{pl}}$, denoted by a dashed line.

###  3.1 Square pulse

The square pulse integrates the solution from the plane pulse over a finite
spatial range. We consider a square pulse of width $x_{0}$ and magnitude one
centered on the origin,

$S(x,t)=\Theta(x_{0}-|x|)\delta(t).$ (11)

This source can also be written as initial condition
$\phi(x,t=0,\mu)=\Theta(x_{0}-|x|)$. With Eq. (11) as the source, Eq. (10)
gives the uncollided solution:

$\phi_{\mathrm{u}}^{\mathrm{sp}}(x,t)=\int^{x_{0}}_{-x_{0}}\\!ds\,\frac{e^{-t}}{2t}\Theta\left(1-|\eta^{\prime}|\right),$
(12)

where we have defined $\eta^{\prime}$ to contain the integration variable,

$\eta^{\prime}=\frac{x-s}{t}.$ (13)

The Dirac delta function in the source removes the integration over $\tau$.
With the integration limits changed to $[-x_{0},x_{0}]$, the step function in
the source is always unity. Evaluating the integral in Eq. (12) requires
considering all of the possible relationships of the parameters $x,t$ and
$x_{0}$. With these cases considered and the integral evaluated, the
uncollided solution can be written as a piecewise function:

$\phi_{\mathrm{u}}^{\mathrm{sp}}(x,t)=\begin{cases}0&|x|-t>x_{0}\\\
x_{0}\frac{e^{-t}}{t}&t>x_{0}\>\&\>x_{0}-t\leq x\leq t-x_{0}\\\ e^{-t}&t\leq
x_{0}\>\&\>t-x_{0}\leq x\leq x_{0}-t\\\
\frac{e^{-t}(t+x+x_{0})}{2t}&-t-x_{0}<x<t+x_{0}\>\&\>x_{0}+x\leq t\leq
x_{0}-x\\\ \frac{e^{-t}(t-x+x_{0})}{2t}&-t-x_{0}<x<t+x_{0}\>\&\>x_{0}-x\leq
t\leq x_{0}+x\end{cases}$ (14)

Eq. (10) with a square pulse as the source and the plane pulse collided
solution as the kernel gives,

$\phi_{\mathrm{c}}^{\mathrm{sp}}(x,t)=\int^{x_{0}}_{-x_{0}}\\!ds\,\int^{\pi}_{0}\\!du\,F_{1}(x,s,t,u)\,\Theta\left(1-|\eta^{\prime}|\right),$
(15)

with identical simplifications to the integrals as with the uncollided case.
The integrand function is defined as

$F_{1}(x,s,t,u)=\frac{c\,e^{-t}}{8\pi}\left(1-\eta^{\prime
2}\right)\mathrm{sec}^{2}\left(\frac{u}{2}\right)\mathrm{Re}\left[\xi^{2}e^{\frac{ct}{2}(1-\eta^{\prime
2})\xi}\right],$ (16)

with $\xi$ given by Eq. (8). We have not determined a simple, closed form for
the integrals in Eq. (15), though we note that the integrand appears to be a
well-behaved function, and we have had no trouble performing the integration.

The solutions for this problem are shown in Figure 2 for times $t=1,5,$ and
$10$ when the source width parameter is $x_{0}=0.5$. As in the plane pulse
case we also show the uncollided solution at $t=1$. Though the uncollided
solution is still exponentially decaying, we notice that the uncollided
solution approaches zero linearly, and is continuous, but nonsmooth at
$|x|=x_{0}.$

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 2: Square pulse scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{sp}}+\phi_{\mathrm{c}}^{\mathrm{sp}}$, for $c=1$
and $x_{0}=0.5$ at several times; panel (a) also contains the uncollided
scalar flux, $\phi_{\mathrm{u}}^{\mathrm{sp}}$, denoted by a dashed line.

###  3.2 Gaussian pulse

Next, we consider an initial pulse with a Gaussian spatial profile with
standard deviation $\sigma$ centered on the origin,

$S(x,t)=\exp\left({\frac{-x^{2}}{\sigma^{2}}}\right)\,\delta(t).$ (17)

As in Section 3.1, the integration over time simplifies and the uncollided
solution from Eq. (10) becomes,

$\phi_{\mathrm{u}}^{\mathrm{gp}}(x,t)=\int^{\infty}_{-\infty}\\!ds\,\frac{e^{-t}}{2t}\,\exp\left({\frac{-s^{2}}{\sigma^{2}}}\right)\Theta\left(1-|\eta^{\prime}|\right).$
(18)

Equation (18) may be solved analytically. The step function defined in the
Green’s kernel changes the integration limits to $x-t$ and $x+t$ and the
solution is

$\phi_{\mathrm{u}}^{\mathrm{gs}}(x,t)=\sigma\,\sqrt{\pi}\,e^{-t}\,\frac{\text{erf}\left(\frac{t-x}{\sigma}\right)+\text{erf}\left(\frac{t+x}{\sigma}\right)}{4t}.$
(19)

Unlike the solution for a plane pulse or a square pulse, the uncollided
angular flux induced by Eq. (17) is smooth, which has implications for the
convergence of numerical solvers.

The expression for the collided flux in this configuration is very similar to
Eq. (15), but with a different source term and integration limits,

$\phi_{\mathrm{c}}^{\mathrm{gp}}(x,t)=\int^{\infty}_{-\infty}\\!ds\,\int^{\pi}_{0}\\!du\,\exp\left({\frac{-s^{2}}{\sigma^{2}}}\right)\,F_{1}(x,s,t,u)\,\Theta\left(1-|\eta^{\prime}|\right),$
(20)

where $F_{1}$ is given by Eq. (16) and $\eta^{\prime}$ by Eq. (13). In this
case, finding the effective integration limits is simple. Solving
$|\eta^{\prime}|=1$ for $s$ gives the effective integration limits,

$\phi_{\mathrm{c}}^{\mathrm{gp}}(x,t)=\int^{x+t}_{x-t}\\!ds\,\int^{\pi}_{0}\\!du\,\exp\left({\frac{-s^{2}}{\sigma^{2}}}\right)\,F_{1}(x,s,t,u).$
(21)

Figure 3 shows the solution to this problem at times $t=1,5$ and $10$. In this
case all of the solutions, including the uncollided solutions, are smooth.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 3: Gaussian pulse scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{gp}}+\phi_{\mathrm{c}}^{\mathrm{gp}}$, for $c=1$
and $\sigma=0.5$ at several times; panel (a) also contains the uncollided
scalar flux, $\phi_{\mathrm{u}}^{\mathrm{gp}}$, denoted by a dashed line.

###  3.3 Square source

We now consider sources that are nonzero for a finite length of time. The
source can be understood as superpositions of pulses, like those presented in
Sections 3.1 and 3.2 from $t=0$ to $t=t_{0}$. This method uncovers
difficulties in the Green’s kernel that were overlooked in the construction of
solutions for pulses. Namely, the expression $\frac{e^{-t}}{2t^{n}}$, where
$n$ is any integer, that appears in the uncollided and collided kernels, is
singular as $t$ approaches zero and the step function
$\Theta(1-|\eta^{\prime}|)$ behaves more erratically. This erratic behavior
causes integrands to be poorly behaved in numerical integration since the
effective domain changes wildly which can cause quadrature points to be wasted
on regions where the integrand is zero and miss important features.

We first consider a source of width $x_{0}$ and magnitude one turned off at
$t=t_{0}$ centered on the origin,

$S(x,t)=\Theta(x_{0}-|x|)\,\Theta(t_{0}-t).$ (22)

Using this source in Eq. (10) we arrive at the following integral to define
the uncollided square source solution

$\phi_{\mathrm{u}}^{\mathrm{ss}}(x,t)=\int^{\mathrm{min}(t,t_{0})}_{0}\\!d\tau\,\int^{x_{0}}_{-x_{0}}\\!ds\,\frac{e^{-(t-\tau)}}{2(t-\tau)}\,\Theta\left(1-|\eta^{\prime\prime}|\right),$
(23)

with $\eta^{\prime\prime}$ that contains the time integration variable,

$\eta^{\prime\prime}=\frac{x-s}{t-\tau}.$ (24)

As in Section 3.2, the spatial step function in the source has been absorbed
into the integration limit over $s$. The time dependent step function in the
source is always one for the integration limits of Eq. (23). We note that the
solution for the inner integral over $s$ of Eq. (23) is Eq. (14) where $t$ is
replaced with $t-\tau$. We have,

$\phi_{\mathrm{u}}^{\mathrm{ss}}(x,t)=\int^{\mathrm{min}(t,t_{0})}_{0}\\!d\tau\,\phi_{\mathrm{u}}^{\mathrm{sp}}(x,t-\tau).$
(25)

Solving Eq. (25) is an exercise of accounting for all of the possible values
for $\tau$. For example, the first case in the square pulse solution, the
scalar flux is zero if $|x|-t>x_{0}$. Solving $|x|-(t-\tau)>x_{0}$ for $\tau$
tells us that the square source solution is zero if $\tau>t+x_{0}-|x|$.
Considering all of the cases in Eq. (14) allows us to come to an analytic
solution for the uncollided flux from a square source,

$\phi_{\mathrm{u}}^{\mathrm{sp}}(x,t)=\left[-x_{0}\mathrm{E_{i}}(\tau-t)\right]\bigg{|}_{0}^{b}+\frac{1}{2}\left[(|x|-x_{0})\text{Ei}(\tau-t)+e^{\tau-t}\right]\bigg{|}_{b}^{c}+\left[e^{-(t-\tau)}\right]\bigg{|}_{c}^{d}.$
(26)

Where the evaluation intervals of $\tau$ are defined by

$b=\left[\mathrm{min}\left(d,t-|x|-x_{0}\right)\right]_{+},$ (27)
$c=\left[\mathrm{min}\left(d,t+|x|-x_{0}\right)\right]_{+},$ (28)
$d=\left[\mathrm{min}\left(t_{0},t,t-|x|+x_{0}\right)\right]_{+},$ (29)

where $[\cdot]_{+}$ returns the positive part of its argument, and
$\mathrm{E_{i}}$ is the exponential integral.

The expression for the collided solution takes the form of an integral over
Eq. (15):

$\phi_{\mathrm{c}}^{\mathrm{ss}}(x,t)=\int^{\mathrm{min}(t,t_{0})}_{0}\\!d\tau\,\int^{x_{0}}_{-x_{0}}\\!ds\,F_{2}(x,s,t,\tau,u)\,\Theta\left(1-|\eta^{\prime\prime}|\right),$
(30)

where the integrand $F_{2}$ is slightly different from $F_{1}$ from the square
pulse case,

$F_{2}(x,s,t,\tau,u)=\frac{ce^{-(t-\tau)}}{8\pi}\left(1-\eta^{\prime\prime
2}\right)\int^{\pi}_{0}\\!du\,\mathrm{sec}^{2}\left(\frac{u}{2}\right)\mathrm{Re}\left[\xi^{2}e^{\frac{c(t-\tau)}{2}(1-\eta^{\prime\prime
2})\xi}\right],$ (31)

with $\eta^{\prime\prime}$ given by Eq. (24).

While we were able to avoid singularity difficulties by integrating
analytically in finding the solution for the uncollided case in this
configuration, the collided case must be integrated numerically. Due to the
behavior of the step function, the effective integration domain of Eq. (30)
varies drastically with $\tau$. Therefore, the integral is not well-suited to
numerical integration. Switching the order of integration and merging the step
function with the integration limits over $\tau$ gives

$\phi_{\mathrm{c}}^{\mathrm{ss}}(x,t)=\int^{x_{0}}_{-x_{0}}\\!ds\,\int^{\mathrm{min}(t,t-|x-s|)_{+}}_{0}\\!d\tau\,F_{2}(x,s,t,\tau,u),$
(32)

and allows us to cast the integral in a form that we have found to converge at
an reasonable rate.

For this problem the solutions, as shown in Figure 4 for $x_{0}=0.5$ and
$t_{0}=5$, are smoother than the square pulse solutions shown previously.
Also, because the source is one until $t=5$, there is a noticeable uncollided
solution at that later time.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 4: Square source scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{ss}}+\phi_{\mathrm{c}}^{\mathrm{ss}}$, for $c=1$,
$t_{0}=5$ and $x_{0}=0.5$ at several times; panels (a) and (b) also contain
the uncollided scalar flux, $\phi_{\mathrm{u}}^{\mathrm{ss}}$, denoted by a
dashed line.

##  4 Gaussian source

We next consider a Gaussian source with standard deviation $\sigma$ that is
turned off at time $t_{0}$ where $S$ is given as

$S(x,t)=\exp\left(\frac{-x^{2}}{\sigma^{2}}\right)\,\Theta(t_{0}-t).$ (33)

Like the square pulse and the square source, the Gaussian source can be
considered a superposition of Gaussian pulses, as Eq. (33) is a superposition
of pulses defined in Eq. (17). Therefore, the uncollided scalar flux solution
can be found with a variable change from $t$ to $t-\tau$ in Eq. (19) and
integration over the time that the source is on,

$\phi_{\mathrm{u}}^{\mathrm{gs}}(x,t)=\int_{0}^{\mathrm{min}(t,t_{0})}\\!d\tau\,\sigma\,\sqrt{\pi}\,e^{-(t-\tau)}\,\frac{\text{erf}\left(\frac{t-\tau-x}{\sigma}\right)+\text{erf}\left(\frac{t-\tau+x}{\sigma}\right)}{4(t-\tau)}.$
(34)

Solving this integral requires numerical integration. While this integral
involves evaluating the integrand when $\tau=t$, the behavior of the error
function allows the integrand to be well behaved. Like the pulse source, this
uncollided scalar flux is associated with a smooth angular flux solution.

To find the collided flux, we integrate Eq. (20) over time,

$\phi_{\mathrm{c}}^{\mathrm{gs}}(x,t)=\int_{0}^{\mathrm{min}(t,t_{0})}\\!d\tau\,\int^{\infty}_{-\infty}\\!ds\,\exp\left({\frac{-s^{2}}{\sigma^{2}}}\right)\,F_{2}(x,s,t,\tau,u)\,\Theta\left(1-|\eta^{\prime\prime}|\right).$
(35)

Just as with the square source case (Section 3.3), the integration orders for
$\tau$ and $s$ are switched to find a better behaved integrand,

$\phi_{\mathrm{c}}^{\mathrm{gs}}(x,t)=\int^{\infty}_{-\infty}\\!ds\,\int^{\mathrm{min}(t,t_{0})}_{0}\\!d\tau\,\exp\left({\frac{-s^{2}}{\sigma^{2}}}\right)\,F_{2}(x,s,t,\tau,u)\Theta(1-|\eta^{\prime\prime}|).$
(36)

To find a more efficient integration interval over $s$, we unravel
$\eta^{\prime\prime}$ to get $t-\tau-x\leq s\leq t-\tau+x$. With the floor of
$\tau$ defined so that $\tau\geq 0$, new integration limits for $s$ are found
and the most efficient form of the integral for the collided solution is,

$\phi_{\mathrm{c}}^{\mathrm{gs}}(x,t)=\int^{x+t}_{x-t}\\!ds\,\int^{\mathrm{min}(t,t_{0},t-|x-s|)_{+}}_{0}\\!d\tau\,\exp\left({\frac{-s^{2}}{\sigma^{2}}}\right)\,F_{2}(x,s,t,\tau,u),$
(37)

where $\mathrm{min}(\cdot)_{+}$ returns the minimum of its arguments or zero
if the minimum is negative. Figure 5 shows the solutions for this problem with
$\sigma=0.5$ and $t_{0}=5$. These solutions resemble the square source
solutions, especially at the later times. At $t=1$ the Gaussian source
solution appears to be narrower than that for the square source solution.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 5: Gaussian source scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{gs}}+\phi_{\mathrm{c}}^{\mathrm{gs}}$, for $c=1$,
$t_{0}=5$ and $\sigma=0.5$ at several times; panels (a) and (b) also contain
the uncollided scalar flux, $\phi_{\mathrm{u}}^{\mathrm{gs}}$, denoted by a
dashed line.

##  5 Cylindrical geometry

Ganapol et al. (2001) provides the solution for cylindrical geometry where the
source is an infinite line pulse aligned with the $z$ axis. This solution is
obtained using the plane-to-point transform and integration of the point
source over an infinite line. The resulting uncollided expression is

$\phi_{\mathrm{u}}^{\mathrm{l}}(r,t)=\frac{e^{-t}}{2\pi
t^{2}}\frac{1}{\sqrt{1-\eta_{p}^{2}}}\Theta\left(1-\eta_{p}\right).$ (38)

Here the superscript $\mathrm{l}$ denotes a line source, $r$ is the radial
coordinate, and

$\eta_{p}\equiv\frac{r}{t},$ (39)

where the $p$ subscript stands for “polar”. The absolute value in the step
function becomes irrelevant in this geometry and is discarded. Unlike the
solutions presented thus far where there were no singularities apart from
$t=0$, Eq. (38) is singular as $\eta_{p}$ approaches one. The collided flux is
found by integrating the collided flux for a point source, which Ganapol also
provides,

$\phi_{\mathrm{c}}^{\mathrm{l}}(r,t)=\left[2t\int_{0}^{\sqrt{1-\eta_{p}^{2}}}\\!d\omega\,\phi_{\mathrm{c}}^{\mathrm{pt}}\left(t\sqrt{\eta_{p}^{2}+\omega^{2}},t\right)\right]\Theta\left(1-\eta_{p}\right),$
(40)

where the point source collided flux is,

$\phi_{\mathrm{c}}^{\mathrm{pt}}(r,t)=\Theta(1-\eta_{p})\times\\\
\left(\frac{e^{-t}}{4\pi
rt^{2}}(ct)\log\left[\frac{1+\eta_{p}}{1-\eta_{p}}\right]+\frac{1}{2\pi}\frac{e^{-t}}{4\pi
rt^{2}}\left(\frac{ct}{2}\right)^{2}\left(1-\eta_{p}^{2}\right)\int_{0}^{\pi}\\!du\,\sec^{2}\left(\frac{u}{2}\right)\mathrm{Re}\left[\left(\eta_{p}+i\tan\left(\frac{u}{2}\right)\right)\xi^{3}e^{\frac{ct}{2}\left(1-\eta_{p}^{2}\right)\xi}\right]\right),$
(41)

where the superscript $\mathrm{pt}$ is short for point, $\eta_{p}$ is given by
Eq. (39), and $\xi$ by Eq. (8). The step function is redundant since it has
been absorbed into the integration limits of Eq. (40).

The line source solution is useful for verification of 2-D transport codes.
However, as we can see from the solutions in Figure 6, there is a singularity
at the wavefront that is still present at $t=5$ on the scale of the figure.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 6: Line pulse scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{l}}+\phi_{\mathrm{c}}^{\mathrm{l}}$, for $c=1$, at
several times; panels (a) and (b) also contain the uncollided scalar flux,
$\phi_{\mathrm{u}}^{\mathrm{l}}$, denoted by a dashed line.

###  5.1 Gaussian pulse

We consider an infinite cylindrical Gaussian pulse of standard deviation
$\sigma$,

$S(r,t)=\exp\left(-\frac{r^{2}}{\sigma^{2}}\right)\delta(t).$ (42)

For the uncollided flux, it is actually necessary to transform back into
Cartesian coordinates since the integration in polar coordinates causes the
integrand to be badly behaved. Using the relationship $r^{2}=x^{2}+y^{2}$, we
introduce $r^{\prime}$,

$r^{\prime 2}=(x-s)^{2}+(y-v)^{2},$ (43)

Where $s$ and $v$ are dummy variables that are integrated over Cartesian
space. Defining a new $\eta$ for the Green’s kernel,

$\eta_{p}^{\prime}=\frac{r^{\prime}}{t}.$ (44)

Now the uncollided flux for a Gaussian pulse may be written in integral form,

$\phi_{\mathrm{u}}^{\mathrm{gp}}(x,y,t)=\frac{e^{-t}}{2\pi
t^{2}}\int^{\infty}_{-\infty}\\!dv\,\int^{\infty}_{-\infty}\\!ds\,\frac{1}{\sqrt{1-\eta_{p}^{\prime
2}}}\exp\left(-\frac{s^{2}+v^{2}}{\sigma^{2}}\right)\Theta\left(1-\eta^{\prime}_{p}\right).$
(45)

The integrand of Eq. (45) is still poorly behaved, but assimilating the step
function into the integration limits of the integral over $s$ will cast it
into a well behaved form. This is done by finding the roots of $s$ for the
equation, $\eta^{\prime}_{p}=1$. The expression for the uncollided flux
becomes,

$\phi_{\mathrm{u}}^{\mathrm{gp}}(x,y,t)=\frac{e^{-t}}{2\pi
t^{2}}\int^{\infty}_{-\infty}\\!dv\,\int^{s_{b}}_{s_{a}}\\!ds\,\frac{1}{\sqrt{1-\eta_{p}^{\prime
2}}}\exp\left(-\frac{s^{2}+v^{2}}{\sigma^{2}}\right),$ (46)

where the integration limits over $s$ are,

$s_{a}=x-\sqrt{[t^{2}-v^{2}+2vy-y^{2}]_{+}}$ (47)
$s_{b}=x+\sqrt{[t^{2}-v^{2}+2vy-y^{2}]_{+}}$ (48)

where $[\cdot]_{+}$ returns the positive part of its argument. Since the
solution is symmetric about the pole, it is not necessary to integrate Eq.
(46) over a two dimensional domain. We can choose $y=0$ and find the
uncollided solution as a function of $r$.

$\phi_{\mathrm{u}}^{\mathrm{gp}}(r,t)=\frac{e^{-t}}{2\pi
t^{2}}\int^{\infty}_{-\infty}\\!dv\,\int^{s_{b}}_{s_{a}}\\!ds\,\frac{1}{\sqrt{1-\eta_{p}^{\prime\prime
2}}}\exp\left(-\frac{s^{2}+v^{2}}{\sigma^{2}}\right),$ (49)

where

$\eta_{p}^{\prime\prime}=\frac{\sqrt{(r-s)^{2}+v^{2}}}{t}.$ (50)

The expression for the collided flux is better behaved in polar coordinates,
where the variables over which the Green’s kernel is integrated become $\rho$
and $\theta^{\prime}$. With this transformation, we define a new radius in
polar coordinates,

$r^{\prime}=\sqrt{(r\cos(\theta)-\rho\cos(\theta^{\prime}))^{2}+(r\sin(\theta)-\rho\sin(\theta^{\prime}))^{2}}.$
(51)

$\eta_{p}^{\prime}$ is given by Eq. (44) with Eq. (51) as $r^{\prime}$. Since
the solution will be symmetric about $r$, the value of angular coordinate
$\theta$ is arbitrary. However, to properly integrate the solution kernel, the
angular coordinate, $\theta^{\prime}$, must be independent from $\theta$.
Therefore, the collided flux requires integration over angle and radius,

$\phi_{\mathrm{c}}^{\mathrm{gp}}(r,\theta,t)=2t\int^{2\pi}_{0}\\!d\theta^{\prime}\,\int^{\infty}_{0}\\!d\rho\,\rho\exp\left(-\frac{\rho^{2}}{\sigma^{2}}\right)\left[\int_{0}^{\sqrt{1-\eta_{p}^{\prime
2}}}\\!d\omega\,\phi_{\mathrm{c}}^{\mathrm{pt}}\left(t\sqrt{\eta_{p}^{\prime
2}+\omega^{2}},t\right)\right]\Theta\left(1-\eta_{p}^{\prime}\right).$ (52)

Evaluating Eq. (52) requires four integrals over a difficult integrand. To
simplify the integral, we first recast it in a more explicit form where the
step function from the point source is brought out of the function,

$\phi_{\mathrm{c}}^{\mathrm{gp}}(r,\theta,t)=2t\int^{2\pi}_{0}\\!d\theta^{\prime}\,\int^{\infty}_{0}\\!d\rho\,\int^{\sqrt{1-\eta_{p}^{\prime
2}}}_{0}\\!d\omega\,\int^{\pi}_{0}\\!du\,Q(\rho)\,\rho\,F^{\mathrm{pt}}_{2}\left(t\sqrt{\eta_{p}^{\prime
2}+\omega^{2}},t\right)\Theta(1-\eta^{\prime}_{p})+\\\
2t\int^{2\pi}_{0}\\!d\theta^{\prime}\,\int^{\infty}_{0}\\!d\rho\,\int^{\sqrt{1-\eta_{p}^{\prime
2}}}_{0}\\!d\omega\,Q(\rho)\,\rho\,F^{\mathrm{pt}}_{1}\left(t\sqrt{\eta_{p}^{\prime
2}+\omega^{2}},t\right)\Theta(1-\eta_{p}^{\prime}),$ (53)

Where $F_{1}^{\mathrm{pt}}$ is the first collided kernel for a point source
without the step function,

$F^{\mathrm{pt}}_{1}(r,t)=\frac{e^{-t}}{4\pi
rt^{2}}(ct)\log\left[\frac{1+\eta_{p}}{1-\eta_{p}}\right],$ (54)

and $F_{2}^{\mathrm{pt}}$ is the integrand for the second on to infinite
collided solution for a point source without the step function,

$F_{2}^{\mathrm{pt}}(r,t)=\frac{1}{2\pi}\frac{e^{-t}}{4\pi
rt^{2}}\left(\frac{ct}{2}\right)^{2}\left(1-\eta_{p}^{2}\right)\sec^{2}\left(\frac{u}{2}\right)\mathrm{Re}\left[\left(\eta_{p}+i\tan\left(\frac{u}{2}\right)\right)\xi^{3}e^{\frac{ct}{2}\left(1-\eta_{p}^{2}\right)\xi}\right],$
(55)

where $\eta_{p}$ is given by Eq. (39). The source, $Q(\rho)$ is

$Q(\rho)=\exp\left(-\frac{\rho^{2}}{\sigma^{2}}\right).$ (56)

The step function, $\Theta(1-\eta_{p}^{\prime})$, causes the effective
integration domain to be erratic and the integrand to be badly behaved. As
with the uncollided case, solving $\eta_{p}^{\prime}=1$ gives the upper and
lower bounds of the integration for $\rho$ as dictated by the step function.
Now a convergent form of collided flux can be found

$\phi_{\mathrm{c}}^{\mathrm{gp}}(r,\theta,t)=2t\int^{2\pi}_{0}\\!d\theta^{\prime}\,\int^{\rho_{b}}_{\rho_{a}}\\!d\rho\,\int^{\sqrt{1-\eta_{p}^{\prime
2}}}_{0}\\!d\omega\,\int^{\pi}_{0}\\!du\,Q(\rho)\,\rho\,F^{\mathrm{pt}}_{2}\left(t\sqrt{\eta_{p}^{\prime
2}+\omega^{2}},t\right)\Theta(1-\eta^{\prime}_{p})+\\\
2t\int^{2\pi}_{0}\\!d\theta^{\prime}\,\int^{\rho_{b}}_{\rho_{a}}\\!d\rho\,\int^{\sqrt{1-\eta_{p}^{\prime
2}}}_{0}\\!d\omega\,Q(\rho)\,\rho\,F^{\mathrm{pt}}_{1}\left(t\sqrt{\eta_{p}^{\prime
2}+\omega^{2}},t\right)\Theta(1-\eta_{p}^{\prime}),$ (57)

where

$\rho_{a}=r\cos(\theta-\theta^{\prime})-\sqrt{\left[\frac{r^{2}\cos(2(\theta-\theta^{\prime}))-r^{2}}{2}+t^{2}\right]_{+}},$
(58)
$\rho_{b}=r\cos(\theta-\theta^{\prime})+\sqrt{\left[\frac{r^{2}\cos(2(\theta-\theta^{\prime}))-r^{2}}{2}+t^{2}\right]_{+}},$
(59)
$\eta^{\prime}_{p}=\frac{\sqrt{(r\cos(\theta)-\rho\cos(\theta^{\prime}))^{2}+(r\sin(\theta)-\rho\sin(\theta^{\prime}))^{2}}}{t}.$
(60)

This solution with $\sigma=0.5$ and $t_{0}=5$ is shown in Figure 7.
Integrating the source over a finite spatial range leads to a removal of the
singularity in the solution at the wavefront. This leads to a smooth solution
at all times.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 7: Cylindrical Gaussian source scalar flux solutions,
$\phi_{\mathrm{u}}^{\mathrm{gs}}+\phi_{\mathrm{c}}^{\mathrm{gs}}$, for $c=1$,
$t_{0}=5$ and $\sigma=0.5$ at several times; panel (a) also contains the
uncollided scalar flux, $\phi_{\mathrm{u}}^{\mathrm{gs}}$, denoted by a dashed
line.

##  6 Conclusion

We have presented the uncollided and collided solutions for several different
source shapes that are either pulsed or turned-on for a finite time. In many
cases we were able to produce a simple formula for the uncollided solution so
that could be used as a source in a computer code for verification purposes in
a similar manner to MMS. We believe that these solutions will be useful for
code verification for a variety of transport codes.

All of the solutions presented were checked for systematic errors with an
independent numerical solver that discretized the transport equation. A
package that evaluates all of the presented solutions along with precomputed
numerical results for $t=1,\,5,\,\mathrm{and}\,10$ are available on
GitHub222github.com/wbennett39/transport_benchmarks.. The code can be used to
modify the parameters in the solutions to produce different benchmarks.

## Acknowledgements

The authors would like to thank Minwoo Shin for providing independent
numerical solutions to the cylindrical Gaussian pulse (Sec. 5). This work was
supported by a Los Alamos National Laboratory under contract #599008,
“Verification for Radiation Hydrodynamics Simulations”.

## References

  * Alcouffe, O’Dell, and Brinkley (1990) Alcouffe, R. E., R. D. O’Dell, and F. W. Brinkley. 1990. “A First-Collision Source Method That Satisfies Discrete Sn Transport Balance.” _Nuclear Science and Engineering_ 105 (2): 198–203.
  * Ganapol et al. (2001) Ganapol, B. D., R. S. Baker, J. A. Dahl, and Raymond E. Alcouffe. 2001. “Homogeneous Infinite Media Time-Dependent Analytical Benchmarks.” https://www.osti.gov/biblio/975281.
  * Garrett and Hauck (2013) Garrett, C Kristopher, and Cory D Hauck. 2013. “A comparison of moment closures for linear kinetic transport equations: The line source benchmark.” _Transport Theory and Statistical Physics_ 42 (6-7): 203–235.
  * Harel, Burov, and Heizler (2021) Harel, Re’em, Stanislav Burov, and Shay I. Heizler. 2021. “The Time-Dependent Asymptotic PN Approximation for the Transport Equation.” _Nuclear Science and Engineering_ 195 (6): 578–597.
  * Hauck and Heningburg (2019) Hauck, Cory, and Vincent Heningburg. 2019. “Filtered discrete ordinates equations for radiative transport.” _Journal of Scientific Computing_ 80 (1): 614–648.
  * Hauck and McClarren (2013) Hauck, Cory D., and Ryan G. McClarren. 2013. “A Collision-Based Hybrid Method for Time-Dependent, Linear, Kinetic Transport Equations.” _Multiscale Modeling & Simulation_ 11 (4): 1197–1227.
  * Heizler (2010) Heizler, Shay I. 2010. “Asymptotic Telegrapher’s Equation (P1) Approximation for the Transport Equation.” _Nuclear Science and Engineering_ 166 (1): 17–35.
  * Heningburg and Hauck (2020) Heningburg, Vincent, and Cory D Hauck. 2020. “Hybrid solver for the radiative transport equation using finite volume and discontinuous galerkin.” _arXiv preprint arXiv:2002.02517_ .
  * McClarren and Lowrie (2008) McClarren, Ryan G, and Robert B Lowrie. 2008. “Manufactured solutions for the P1 radiation-hydrodynamics equations.” _Journal of Quantitative Spectroscopy and Radiative Transfer_ 109 (15): 2590–2602.
  * Monin (1956) Monin, A.S. 1956. “A statistical interpretation of the scattering of microscopic particles.” _Theory of Probability and its Applications_ 1 (3): 298–311.
  * Peng and McClarren (2021) Peng, Zhuogang, and Ryan G McClarren. 2021. “A high-order/low-order (HOLO) algorithm for preserving conservation in time-dependent low-rank transport calculations.” _Journal of Computational Physics_ 447: 110672.
  * Salari and Knupp (2000) Salari, Kambiz, and Patrick Knupp. 2000. _Code verification by the method of manufactured solutions_. Technical Report. Sandia National Lab.(SNL-NM), Albuquerque, NM (United States); Sandia ….
  * Schlachter and Schneider (2018) Schlachter, Louisa, and Florian Schneider. 2018. “A hyperbolicity-preserving stochastic Galerkin approximation for uncertain hyperbolic systems of equations.” _Journal of Computational Physics_ 375: 80–98.
  * Seibold and Frank (2014) Seibold, Benjamin, and Martin Frank. 2014. “StaRMAP—A second order staggered grid method for spherical harmonics moment equations of radiative transfer.” _ACM Transactions on Mathematical Software (TOMS)_ 41 (1): 1–28.
  * Variansyah and McClarren (2022) Variansyah, Ilham, and Ryan G McClarren. 2022. “Population Control Techniques for Time-Dependent and Eigenvalue Monte Carlo Neutron Transport Calculations.” _arXiv preprint arXiv:2202.08631_ .
  * Walters and Haghighat (2017) Walters, William J., and Alireza Haghighat. 2017. “The adaptive collision source method for discrete ordinates radiation transport.” _Annals of Nuclear Energy_ 105: 45–58.
  * Wright (2006) Wright, Matthew CM. 2006. “Green function or green’s function?” _Nature Physics_ 2 (10): 646–646.


# Viewpoints: A high-performance high-dimensional exploratory data analysis
tool

P.R. Gazis  SETI Institute, Mountain View, CA, 94043<EMAIL_ADDRESS>C.
Levit  NASA Ames Research Center, Moffett Field, CA 94035
<EMAIL_ADDRESS>M.J. Way11affiliation: NASA Ames Research Center, Space
Science Division, Moffett Field, CA 94035, USA 22affiliation: Department of
Astronomy and Space Physics, Uppsala, Sweden NASA Goddard Institute for Space
Studies, 2880 Broadway, New York, NY, 10025<EMAIL_ADDRESS>
###### Abstract

Scientific data sets continue to increase in both size and complexity. In the
past, dedicated graphics systems at supercomputing centers were required to
visualize large data sets, but as the price of commodity graphics hardware has
dropped and its capability has increased, it is now possible, in principle, to
view large complex data sets on a single workstation. To do this in practice,
an investigator will need software that is written to take advantage of the
relevant graphics hardware. The Viewpoints visualization package described
herein is an example of such software. Viewpoints is an interactive tool for
exploratory visual analysis of large, high-dimensional (multivariate) data. It
leverages the capabilities of modern graphics boards (GPUs) to run on a single
workstation or laptop. Viewpoints is minimalist: it attempts to do a small set
of useful things very well (or at least very quickly) in comparison with
similar packages today. Its basic feature set includes linked scatter plots
with brushing, dynamic histograms, normalization and outlier
detection/removal. Viewpoints was originally designed for astrophysicists, but
it has since been used in a variety of fields that range from astronomy,
quantum chemistry, fluid dynamics, machine learning, bioinformatics, and
finance to information technology server log mining. In this article, we
describe the Viewpoints package and show examples of its usage.

Multivariate Datasets, Visualization

##  1 Introduction

Analysis and visualization of extremely large and complex data sets has become
one of the more significant challenges facing the scientific community today.
Current instruments and simulations produce enormous volumes of data. These
data sets include hyperspectral images from spacecraft, multivariate data of
high dimensionality from sky surveys, time-varying three-dimensional flows
from supercomputer simulations, and complex interrelated time series from
vehicle telemetry, DNA microarrays or physiological monitoring, to name a few.

One representative example from astronomy is the Sloan Digital Sky (SDSS York
et al., 2000). The spectroscopic catalog from the SDSS Data Release 7 (DR7
Abazajian et al., 2009) contains information for over 1.5 million objects. For
each object in the joint spectroscopic and photometric catalog there are over
100 physical measurements. Properties in these catalogs include position on
the sky, spectroscopic redshift, photometric magnitude in five bands,
energy/width of up to 50 spectral lines, galaxy angular size, bulge-to-disk
ratio, etc. In addition, many other physical quantities of interest can be
derived from these data, such as absolute magnitudes, photometric redshifts,
galaxy type, average distance to nearest neighbors and other properties of the
galaxy’s local environment. This is the kind of catalog in which Viewpoints
excels on a relatively modern desktop computer.

On the other hand the DR7 photometric catalog has measured over 350 million
unique objects. This is an enormous amount of data by traditional standards
and these kinds of data volumes will continue to grow in size over the coming
decade. The actual data volume in the SDSS DR7, including pre- and
postprocessing products, is more than 50 TB, while the searchable database for
the photometric and redshift catalogs is over 3.5
TB.111http://www.sdss.org/dr7 We do not claim that Viewpoints is able to
explore a data set of this size all at once as it is beyond the memory
capabilities of a normal desktop computer. Still, Viewpoints can quickly
explore subsets of a large complex data set like this. When Viewpoints does
reach the limits of desktop hardware, one can reach for more capable, but
expensive and rare, custom-built hardware systems such as the hyperwall
(Sandstrom et al., 2003). The hyperwall can then be utilized efficiently using
the knowledge already obtained through exploration via Viewpoints.

Although advances in hardware speed and storage technology have kept up with
(and, indeed, have driven) increases in database size, the same is not true of
the tools we use to explore and understand these data (e.g., SubbaRao et al.,
2008). Modern data sets routinely outstrip the capabilities of contemporary
visualization and analysis software, and this problem can only get worse as
data volumes from future astronomy-related programs (e.g., Large Synoptic
Survey Telescope [LSST]222See http://www.lsst.org, Supernova Acceleration
Probe [SNAP]333See http://snap.lbl.gov, Panoramic Survey Telescope and Rapid
Response System [PanSTARRS]444See http://pan-starrs.ifa.hawaii.edu/public, and
Dark Energy Survey [DES]555See http://www.darkenergysurvey.org) and ever-
larger supercomputer simulations increase by orders of magnitude. See Borne
(2009) for a review. Data sets of this size and dimensionality defy
interactive analysis with standard tools such as MATLAB,666See
http://www.mathworks.com Octave,777See http://www.octave.org
Mathematica,888See http://www.wolfram.com R,999See http://www.r-project.org
S++,101010See http://www.splusplus.com xGobi,111111See
http://www.research.att.com/areas/stat/xgobi Mirage,121212See http://www.bell-
labs.com/project/mirage gnuplot,131313See http://www.gnuplot.info etc. To
address this problem, we have developed Viewpoints.

Viewpoints was inspired to a large extent by the hyperwall (Sandstrom et al.,
2003). This facility, developed by the NASA Advanced Supercomputing division
at NASA Ames Research Center, is a high-performance visualization cluster that
was specifically designed to address the problem of exploring, visualizing,
and analyzing large, complex, multidimensional data. See Wong & Bergeron
(1997) for a review on the history of multidimensional multivariate
visualization. In its original form, the hyperwall was a 50 node Linux cluster
with a 7$\times$7 array of tiled displays that could be used to display data
in a variety of formats that might range from simple plots to elaborate 3D
rendering. It has since been superseded by even larger facilities at NASA Ames
Research Center and several other institutions throughout the world.

Hyperwalls have been applied to a wide range of problems in astrophysics,
Earth science, aerodynamics, and life sciences (MacDougall & Henze, 2003;
Brown et al., 2004; Murman et al., 2004). Unfortunately, tools such as the
hyperwall require custom hardware and a team of onsite specialists at a
dedicated facility, which places them beyond the reach of most investigators
and research programs.

Viewpoints, shown in Figure 1, is intended to provide a small but powerful
subset of the hyperwall’s functionality on any working scientist’s desktop or
laptop computer.

Figure 1: Viewpoints displaying a 2$\times$2 array of scatter plots of SDSS
data. (Upper left) xy positions of $\sim$20,000 galaxies from the SDSS. Colors
indicate different classes identified by a clustering algorithm. (Upper right)
Density classes identified by a Bayesian partitioning scheme vs class
identified by a Self-Organizing Map. (Lower left) Local density vs local
gradient. (Lower right) Local density vs density classes identified by a
Bayesian partitioning scheme. (Far right) The viewpoints control panel.

Viewpoints provides a high-performance interactive multidimensional visual
exploratory data analysis environment that allows users to quickly explore and
visually analyze the most common forms of large multidimensional scientific
data. It can also be run as a “plug in” or external module for scientific
programming environments such as MATLAB.

In contrast to programs like Mirage (Ho, 2007) or VisIVO (Becciani et al.,
2010) which have a large number of supporting features, Viewpoints does a few
things very well. Its basic feature set includes linked scatter plots with
brushing (e.g., Becker & Cleveland, 1987; Davidson & Sardy, 2000; Chen, 2003;
Stump et al., 2004; Co et al., 2005; Jordan et al., 2008), dynamic histograms,
normalization, and outlier detection/removal. It also includes a limited
capability for Boolean operations with brushes analogous to those of Hurter et
al. (2009).

##  2 Overview and Capabilities

Viewpoints is a fast interactive tool for the analysis and visualization of
high-dimensional data primarily via linked scatter plots with brushing. It is
intended for general-purpose visual data analysis, interactive exploratory
data mining, and machine-assisted pattern discovery. The targeted application
is the exploration of large, complex, multivariate data sets, with
dimensionalities up to 100 or more and sample sizes up to 105–107. Viewpoints
takes advantage of hardware-accelerated graphics (graphics processing units
[GPUs]) and other capabilities of modern desktop workstations and laptops
(e.g., multithreading and memory-mapped input/output). As previously above, it
is highly interactive, and almost all of its parameters can be modified in
real time to provide immediate visual feedback.

Viewpoints can be invoked either from a desktop icon or from the command line
with a set of optional parameters. Once it is running, it is controlled by a
conventional GUI and/or configuration files saved from previous sessions.
Viewpoints reads multidimensional data in the form of ASCII files, binary
files, or simple Flexible Image Transport System (FITS) tables and displays
these data as an arbitrary number of simultaneously linked 2D and 3D scatter
plots (e.g., Comparato et al., 2007) with independently controllable viewing
transformations, normalization schemes, overplotting schemes (Zhang et al.,
2003), axis labels, and other formatting. These data can then be visualized
and explored interactively using a set of graphical, statistical, and
semiautomated pattern-recognition techniques. The user can select, flag,
and/or delete subsets of the data using brushes, and the modified data set can
be written back to disk in ASCII, binary, or FITS format.

The program’s capabilities were carefully selected using a minimalist
philosophy. We have chosen to provide a small core of essential features that
are universally useful and keep the program size small, execution speeds high,
and the learning curves shallow. Figure 1 shows an example of Viewpoints
displaying a small 2$\times$2 array of linked scatter plots. This figure is
described in detail in §5. In practice, one can simultaneously display an
arbitrary number of linked scatter plots, and we routinely use up to 30
windows to examine data sets of high dimensionality. Viewpoints has a
conventional GUI, shown in the right panel of Figure 1, with an assortment of
menus, buttons, and sliders that provide access to its features. The most
important of these capabilities can also be invoked by mouse gestures and/or
keyboard shortcuts. It has a full range of basic help facilities, such as
usage messages, help files, and tooltips. It also produces a diagnostic output
stream when invoked via the command line.

Viewpoints has a tightly integrated set of features. A partial list of these
features is provided next. Some were in the original release, others were
added in response to users requests, and a few were contributed to the open
source repository141414See http://www.assembla.com/wiki/show/viewpoints by a
small community of developers.

  1. 1.

Multiple linked scatter plots with dynamic brushing rendered via OpenGL with
update rates of greater than $10^{7}$ points $s^{-1}$ for most modern graphics
cards.

  2. 2.

The ability to brush separate selections using separate colors and/or symbols.

  3. 3.

The ability to map densely overplotted regions smoothly to different hues
and/or intensities.

  4. 4.

The ability to save the program’s complete state to an XML formatted
configuration file at any time.

  5. 5.

The ability to initialize and/or restore the program’s state from a
configuration file.

  6. 6.

The ability to append and merge additional data and create or destroy
additional plot windows on the fly.

  7. 7.

Per-axis one-dimensional equiwidth marginal histograms with dynamically
adjustable bin widths and dynamic update of the selected fraction in each bin.
Histograms are individually selectable and adjustable for both the $x$ and $y$
axes of each plot separately.

  8. 8.

Per-axis automatic normalizations based data limits, quantiles, moments, rank,
or simple Gaussianization.

  9. 9.

Elementary outlier detection and removal.

  10. 10.

The ability to change most parameters globally across all plots
simultaneously, or locally on a single plot.

  11. 11.

The ability to lock (and unlock) individual plots so they are immune from
(susceptible to) global operations.

  12. 12.

Rudimentary semi-automated browsing through high-dimensional data by permuting
all (unlocked) axes on all plots.

  13. 13.

Real-time display of quantitative data for individual axes and for the current
selection.

  14. 14.

Basic online help.

  15. 15.

A comprehensive set of error tests along with graceful error recovery if a
user attempts to read corrupt and ill-formatted data or configuration files.

  16. 16.

The ability to replace missing values with user-specified default values.

  17. 17.

Basic ability to plot time series analogous to that of Akiba & Ma (2007).

  18. 18.

Viewpoints also achieves smooth interactive performance on large data as long
as:

    1. (a)

The (binary) array corresponding to the dataset fits into real system (CPU)
memory

    2. (b)

The coordinates of all vertices being displayed fit into available graphics
(GPU) memory. For example, on a 512 Mbyte graphics card viewpoints can
simultaneously update and display (4
plots)$\times$(8x106points/plot)$\times$(3 vertices/point)$\times$ (4
bytes/vertex)= $\sim$400Megabytes

##  3 Software Design

Viewpoints uses hardware-accelerated OpenGL 1.5 on all platforms and performs
rendering via OpenGL151515http://www.opengl.org vertex buffer objects to
provide extremely high graphics performance. The package has cross-platform
availability and can be easily ported to other architectures. Table 1 has a
list of operating systems and architectures known to work to-date. Viewpoints
has a preliminary interface to MATLAB and Octave. It has the potential to work
with scripting languages such as Python161616See http://www.python.org or
Perl.171717See http://www.perl.org Most importantly, Viewpoints is
specifically engineered for use with large data sets. As noted previously, the
targeted application is the interactive analysis of large, complex,
multivariate data sets, with dimensionalities that may surpass 100 and sample
sizes that may exceed 105–107.

Table 1: Platforms Known To Run viewpoints Operating System | OS Revision | Architecture
---|---|---
Apple | Mac OS 10.3–10.6 | Intel
Apple | Mac OS 10.3–10.6 | PPC
Linux | Most Distributions with Kernels 2.4.x–2.6.x | Intel
Windows | 2000 | Intel
Windows | XP | Intel
Windows | Vista | Intel
Windows | 7 | Intel

In the interests of portability, extensibility and efficiency, Viewpoints is
written in C++. This allows us to program at a high level of abstraction and
still achieve excellent runtime performance on all target platforms. Most
importantly, C++ allows us to leverage vast amounts of code written by other
people. Identical code compiles and runs on Linux, Microsoft Windows, and
Apple OS X.

Viewpoints is a relatively small program. The compiled Linux binary is only
5.5 Mbytes.181818The compiled file contains a large number of statically
linked libraries Since it makes extensive use of libraries including some that
leverage C++ template metaprogramming (e.g.,
Blitz++191919http://www.oonumerics.org/blitz and
Boost202020http://www.boost.org), the Viewpoints source tree itself is also
rather small (approximately 40,000 lines). The libraries it depends on are
listed in Table 2.

Table 2: C++ libraries used in viewpoints Library | Description
---|---
libc++ | The standard runtime library for basic math, input/output, and memory management
STL | The standard template library for data structures and algorithms
OpenGL | The standard cross-platform library for hardware accelerated graphics
FLTK | A lightweight platform-independent open source graphical user interface (GUI) library
Blitz++ | A library of extensions that add FORTRAN 90-style array expressions (and fortran-90 style numerical performance) to c++
GSL | The GNU scientific library. Includes functionality such as statistics and linear algebra
Boost serialization | a library of extensions for saving and restoring the state of c++ objects.

##  4 GPU Usage and OpenGL Implementation

Graphics boards that were available at the time development began, such as
NVIDIA Geforce 7950, can transform and display over a billion floating-point
vertices per second and can drive two or more 4 Mpixel displays
simultaneously. This capability continues to increase. In its default mode,
Viewpoints uses GPUs explicitly to take advantage of this tremendous
processing and rendering power. This allows the package to provide smooth
high-frame-rate interactivity while manipulating multiple linked views, each
containing millions of data points. If for some reason a GPU is not available,
Viewpoints can also be instructed to use a system’s ordinary CPU, at some cost
in performance.

The current release uses the following advanced OpenGL features:

  1. 1.

GPU-resident vertex buffer objects (VBOs) for point coordinates. This
conserves main memory (by offloading relevant data to GPU memory), increases
graphics performance (by caching frequently rendered vertex data locally on
the GPU), and frees up essentially all of the CPU bandwidth and bus bandwidth
for other uses.

  2. 2.

GPU-resident VBOs for point indices. A single copy of the vertex index arrays
(which encode the selection status) is shared by all plots.

  3. 3.

Squared alpha blending with separate blending control per plot. For data with
wide variations in projected density this provides fast, adjustable mapping of
overplot density to varying hue and/or intensity.

  4. 4.

Point sprites for displaying symbols and/or smooth kernels instead of just
simple points (splatting).

  5. 5.

Antialiased points and/or point sprites of variable color, brightness, and
size per brush and/or per plot.

  6. 6.

Floating-point color and alpha values, allowing both more sensitivity and more
dynamic range for overplotting.

These design features provide dramatic graphics performance. When one of our
early beta releases was run on a 1997 MacBook Pro laptop, we measured
throughput of over $70\times 10^{6}$ vertices s-1. This consisted of six
linked plots, each with 2 million data points simultaneously updating at five
frames s-1.

Viewpoints has proved quite useful on small-sized (e.g., netbooks) and medium-
sized computer platforms with as little as 1Gbyte of RAM. We do most of our
demos on a laptop displaying up to 18-27 variables in 6-9 windows on a single
monitor. To take full advantage of this package we recommend using a high-end
workstation with high-end graphics card, two high-resolution displays, and
several gigabytes of RAM. Experience suggests that this is a good impedance
match to the visual perception of a scientist or a team of up to three
collaborators.

The best way to think about the limitations of Viewpoints with respect to
common desktop computer hardware is that the data should be able to literally
fit in RAM memory. Specifically Viewpoints can achieve smooth interactive
performance on data as long as:

  1. 1.

The (binary) array corresponding to the data set fits into real system (CPU)
memory

  2. 2.

The coordinates of all vertices being displayed fit into available graphics
(GPU) memory.

For those looking to the future please note that the relationship between the
amount of memory on a graphics card and the allowable size of the data set
appears to be linear. So the latest graphics cards with 1.5Gbyte GPU memory
should display four plots with $\sim$25$\times 10^{6}$ points per plot, with
an update rate of a few frames per second. The update rate is harder to
estimate, since “vertices per second” is no longer a common figure of merit in
GPU specs. However, this appears to scale roughly along with clock speed and
GPU memory.

When visible plots exhaust GPU memory, one can run Viewpoints with the
“–no_vbo” flag, and check the “defer redraws” box in the lower right-hand
control panel. We have run 100 million point data sets on a laptop this way
and it is quite usable.

##  5 Examples

Figure 1 shows a simple example of Viewpoints usage with a small 2$\times$2
array of linked scatter plots of data from the SDSS (York et al., 2000),
several methods of calculating local density, and a cluster detection scheme.
The upper left panel shows $xy$ coordinates of $\sim$20,000 galaxies from the
SDSS. When used this way, a scatter plot becomes a map of a region of space
that is approximately 30 Mpc across. The blue, yellow, and black colors
indicate cluster galaxies, halo galaxies, and field galaxies identified by a
partitioning scheme based on self-organizing maps (SOMs). It can be seen that
cluster galaxies (blue) are concentrated into filaments while the field
galaxies (black) are distributed more or less uniformly. The upper right panel
shows density classes identified by a Bayesian partitioning scheme versus
classes identified by the SOM-based algorithm. This panel is where the
selection capability of Viewpoints was used. The rectangle in the lower right
side of the panel is the selection box that was used to flag the halo galaxies
(yellow points). The lower left panel shows local density versus local
gradient to display the distribution of the three galaxy classes in parameter
space. The lower right panel shows local density versus density class
identified by a Bayesian partitioning scheme. This panel shows a simple
example of the histogram capability of Viewpoints. Figure 1 provides a real-
world example of the abilities of Viewpoints to quickly deduce scientifically
useful features for research purposes. The example provided here is directly
related to a recent publication by two of the authors of this article (Way,
Gazis & Scargle, 2010). All of the figures comparing the three different
clustering technologies described in (Way, Gazis & Scargle, 2010) were
originally devised using Viewpoints. Publication-quality plots were later
produced using MATLAB.

Figure 2 shows a more elaborate usage involving a 2$\times$2 array of linked
scatter plots of astrometric data from the Tycho mission (Perryman et al.,
1997; Høg et al., 2000). For clarity, the control panel has been eliminated
from this figure.

Figure 2: Viewpoints displaying a 2$\times$2 array of linked scatter plots of
Tycho data for 1 million stars. (Upper left) Right ascension vs declination.
(Upper right) Number of transits vs parallax error, with a selection box
outlined in white. (Lower left) Signal to noise ratio vs a quality flag.
(Lower right) V band magnitude plotted vs the B-V magnitude difference.

The upper left panel shows a plot of right ascension versus declination for 1
million stars. This panel illustrates the effect of overplotting with alpha
blending. The galactic plane, with its greater concentration of stars, is
plainly visible as a curving white band that meanders through the image. The
upper right panel shows a plot of the number of transits versus parallax
error. The selection box (outlined in blue) has been used to select points for
which more than 250 transits were available. The legend at the top of the
panel reports that these involved 32,503 points or 3.14% of the total. This
selection is repeated in the other panels. In the upper left panel described
previously, it appears as two parallel and roughly sinusoidal bands that mark
the part of the spacecraft’s orbit from which it was possible to obtain this
number of transits. The lower left panel shows a plot of signal-to-noise ratio
(S/N) versus a quality flag. This is included to illustrate use of a marginal
(conventional) and a conditional histogram. The latter peaks at a
significantly higher S/N ratio than the former, as would be expected, since
measurements that involve more transits should have a higher S/N. In the lower
right panel the $V$-band magnitude has been plotted against the $B--V$
magnitude difference to produce a typical color-magnitude diagram.

Figure 3 shows an example of Viewpoints used to visualize the set of solar
wind data that is provided with the Viewpoints package for demonstration
purposes.

Figure 3: Viewpoints displaying a 2$\times$3 array of linked scatter plots of
solar wind data. (Upper left) Solar wind speed versus time. (Upper middle)
Solar wind density vs time. (Upper right) Solar wind temperature vs time.
(Lower left) Density vs speed. (Lower middle) Temperature vs speed. (Lower
right) Temperature vs density.

In the top three panels, solar wind speed, density, and temperature have been
plotted versus day since 1950 to produce a conventional time series. These
show a typical structure in the solar wind, in which high-speed/low-density
plasma (visible on the left of each panel) overtakes low-speed/low-density
material (visible on the right) to produce an interaction region (visible in
the middle) that is characterized by intermediate speeds, high density, and
high temperature. These interaction regions are often subdivided into a
leading portion and trailing portion with different ranges of densities and
temperatures separated by a well-defined boundary. These are plainly visible
in this time-series plot. The selection capability of Viewpoints has been used
to select these regions as blue and yellow. The bottom three panels show
scatter plots of velocity versus density, velocity versus temperature, and
density versus temperature. It can be seen that the leading portion of this
interaction region (blue) was associated with a large range of densities and
moderate range of temperatures, while the trailing portion was associated with
a moderate range of densities and large range of temperatures.

As mentioned in the Abstract, Viewpoints has been used extensively in a
variety of fields to advance specific scientific goals. Further examples
include:

  * •

Geophysics, Relationships between earthquake occurrence frequency, location,
season, and time of day.

  * •

Information Technology, Relationships between query response time, query
source location, time of day, type of query, etc. These were causing slow
response and high system loads on commercial servers which viewpoints helped
resolve quickly.

  * •

Finance, Currency arbitrage relationships and portfolio clustering in
statistical arbitrage.

##  6 Opening up to the open source community

Two significant concerns with any software tools are long-term support and
response to the user community. Viewpoints has been released as an open source
and is on a collaborative-oriented public wiki site called Assembla.212121See
http://www.assembla.com We now have a significant and growing user community
that is successfully and enthusiastically applying Viewpoints to a huge
variety of problem spread across multiple disciplines.

##  7 Conclusions and the future of viewpoints

Viewpoints was adapted from its alpha prototype to general use with seed
funding from the NASA Applied Information Systems Research Program.222222See
http://aisrp.nasa.gov It is now available on all major computing platforms
(see Table 1) and is in use by scientists and engineers across all NASA
mission directorates, as well as by some early adopters in industry and
academia. Its operation, interface, source code, and impressive performance
are the same on all platforms. Viewpoints is now directly callable from
MATLAB. There is a Viewpoints wiki and Web site232323See
http://www.assembla.com/wiki/show/viewpoints that allows the open source
community to take part in development and affords users access to the latest
information regarding bug fixes and program enhancements. In the near future
we expect more development to come from the open source community and plan to
modify the base code to take advantages of new advances in video chip
technology as it continues to evolve.

Many thanks go to Chris Henze and Jeff Scargle for their input and ideas. The
authors would also like to thank Joe Bredekamp and his NASA Applied
Information Systems Research Program for pilot funding of Viewpoints. M.J.W.
would like to thank the Astronomy Department of Uppsala University for their
kind hospitality. This research has made use of NASA’s Astrophysics Data
System. Viewpoints has been released under the NASA Open Source Agreement and
is considered an open source license. The source code and precompiled binaries
can be freely downloaded from the Viewpoints Web site. If you use Viewpoints
in your research, please consider referencing this work.

## References

  * Abazajian et al. (2009) Abazajian, K. et al. 2009, ApJS, 182, 543
  * Akiba & Ma (2007) Akiba, H. & Kwan-Liu, M. 2007, “An Interactive Interface for Visualizing Time-Varying Multivariate Volume Data”, Proceedings of the Joint Eurographics-IEEE VGTC Symposium on Visualization, May 2007
  * Becciani et al. (2010) Becciani, U. et al. 2010, PASP, 122, 119
  * Becker & Cleveland (1987) Becker, R.A. & Cleveland, W.S. 1987, Technometrics, 29, 127
  * Borne (2009) Born, K. 2009, arXiv:0911.0505v1
  * Brown et al. (2004) Brown, J., Tobak, M., Prabhu, D. & Sandstrom, T. “Flow Topology About an Orbiter Leading Edge Cavity at STS-107 Reentry Conditions”, AIAA-2004-2285 37th AIAA Thermophysics Conference, Portland, Oregon, 2004.
  * Chen (2003) Chen, H. 2003, IEEE Symp. on Inf. Vis. (Washington: IEEE), 181
  * Co et al. (2005) Co, C.S. et al. 2005, Eurographics/IEEE-VGTC Symp. Vis. (Washington: IEEE) 279
  * Comparato et al. (2007) Comparato, M. et al. 2007, PASP, 119, 898
  * Davidson & Sardy (2000) Davidson, A.C. & Sardy, S. 2000, J. Comput. and Graph. Stat., 9, 750
  * Ho (2007) Ho, T.K. 2007, in Statistical Challenges in Modern Astronomy IV, ASP Conferences Series, Vol. 371, p.391. eds G. Jogesh Babu and Eric D. Feigelson
  * Høg et al. (2000) Høg, E. et al. 2000, A&A, 355, L27
  * Hurter et al. (2009) Hurter, C., Tissoires, B. & Conversy, S. 2009, “FromDaDy: Spreading Aircraft Trajectories Across Views to Support Iterative Queries”, IEEE Transactions on Visualization and Computer Graphics, 15, 1017
  * Jordan et al. (2008) Jordan, D.D. et al. 2008, AAS 08-245, AAS/AIAA Space Flight Mechanics Meeting, (AAS 08-245; Springfield; AAS), 499
  * MacDougall & Henze (2003) MacDougall, P. & Henze C. “Fleshing-out pharmacophores with volume rendering of molecular charge densities and hyperwall visualization technology, 226th ACS National Meeting New York, NY September 07-11, 2003.
  * Murman et al. (2004) Murman, S.M., Aftosimis, M.J. & Nemec, M. “Automated Parameter Studies Using a Cartesian Method”, AIAA 2004-5076 22nd AIAA Applied Aerodynamics Conference, Providence, RI., 2004
  * Perryman et al. (1997) Perryman, M.A.C., et al., 1997, The Hipparcos and Tycho Catalogues, (SP-1200, Noordwijk: ESA), 117
  * Sandstrom et al. (2003) Sandstrom, T.A., Henze, C., & Levit, C. 2003, Coordinated & Multiple Views in Exploratory Visualization (Piscataway: IEEE), 124
  * Stump et al. (2004) Stump, G.M. et al. 2004, Aerospace Conference (Piscataway: IEEE), 6, 3885
  * SubbaRao et al. (2008) SubbaRao, M. U., Aragón-Calvo, M.A., Chen, H.W., Quashnock, J.M., Szalay, A.S. & York, D.G. 2008, NJP, 10, 125015
  * Way, Gazis & Scargle (2010) Way, M.J., Gazis, P.R. & Scargle, J.S. 2010, ApJ, submitted (arXiv:1009.0387v1)
  * Wong & Bergeron (1997) Wong, P.C. & Bergeron, R.D. 1997, Scientific Visualization Overviews Methodologies and Techniques (Washington: IEEE), 3
  * York et al. (2000) York, D.G. et al. 2000, AJ, 120, 1579
  * Zhang et al. (2003) Zhang, J. et al. 2003, IEEE International Conference on Data Mining (Washington: IEEE)