# Scaling relations between viscosity and diffusivity in shear-thickening
suspensions

Abhinendra Singh  Department of Macromolecular Science and Engineering, Case
Western Reserve University, Cleveland, Ohio 44106, USA
<EMAIL_ADDRESS>Kuniyasu Saitoh  Department of Physics, Faculty of
Science, Kyoto Sangyo University, Kyoto 603-8555, Japan<EMAIL_ADDRESS>su.ac.jp

###### Abstract

Dense suspensions often exhibit a dramatic response to large external
deformation. The recent body of work has related this behavior to transition
from an unconstrained lubricated to a constrained frictional state. Here, we
use numerical simulations to study the flow behavior and shear-induced
diffusion of frictional non-Brownian spheres in two dimensions under simple
shear flow. We first show that _both_ viscosity $\eta$ and diffusivity
$D/\dot{\gamma}$ of the particles increase at characteristic shear stress,
which is associated with lubrication to frictional transition. Subsequently,
we propose a one-to-one relation between viscosity and diffusivity using the
length scale $\xi$ associated with the size of collective motions (rigid
clusters) of the particles. We demonstrate that $\eta$ and $D/\dot{\gamma}$
are controlled by $\xi$ in two distinct flow regimes, i.e. in the frictionless
and frictional states, where the one-to-one relation is described as a
crossover from $D/\dot{\gamma}\sim\eta$ (frictionless) to $\eta^{1/3}$
(frictional). We also confirm the proposed power laws are insensitive to the
interparticle friction and system size.

††preprint: Draft

##  I Introduction

Suspension of rigid particles dispersed in a liquid are ubiquitous in natural
Coussot (1997), geophysical Jerolmack and Daniels (2019), industrial Benbow
and Bridgewater (1993), and biological Maheshwari  _et al._ (2019) phenomena.
These suspensions often display non-Newtonian rheology, including shear
thinning, normal stress differences, shear-thickening, or even jamming Morris
(2020); Denn and Morris (2014); Guazzelli and Pouliquen (2018); Ness  _et al._
(2022); Jamali  _et al._ (2020). A typical feature observed in these systems
is the eventual divergence of the viscosity $\eta$ and a fluid-solid
transition beyond a critical volume fraction often called the jamming volume
fraction $\phi_{J}$ Bonn  _et al._ (2017). In systems with particles, the
precise value of $\phi_{J}$ depends sensitively on numerous microscopic
details as well as the history of sample Torquato  _et al._ (2000); Makse  _et
al._ (2000); Silbert  _et al._ (2002); Zhang and Makse (2005); Shundyak  _et
al._ (2007). Significantly, the jamming point is affected by the constraints
on inter-particle relative motion (sliding, rolling, and twisting) Singh  _et
al._ (2020a); Santos  _et al._ (2020) (originating from roughness Pradeep  _et
al._ (2021); Singh  _et al._ (2022); Singh (2022), interfacial chemistry James
_et al._ (2018, 2019)), particle-shape Donev  _et al._ (2004); Salerno  _et
al._ (2018), and polydispersity Hopkins  _et al._ (2013); Guy  _et al._
(2020); Srivastava  _et al._ (2021).

A striking phenomenon commonly observed in diverse systems such as cornstarch
suspensions Fall  _et al._ (2008) and other ideal suspensions Brown and Jaeger
(2009); Guazzelli and Pouliquen (2018); Singh  _et al._ (2022), wherein the
shear stress increases faster than linearly in response to the applied
deformation rate is termed as Shear-thickening (ST). At volume fractions
$\phi$ far below the frictional jamming volume fraction $\phi_{J}^{\mu}$, the
relative viscosity $\eta_{r}$ increases continuously with shear rate, and it
termed as Continuous shear-thickening (CST); here, the relative viscosity
$\eta_{r}(\phi,\dot{\gamma})/\eta_{0}$ is normalized by the solvent viscosity.
At volume fractions close enough to the jamming volume fraction
$\phi_{C}<\phi<\phi_{J}^{\mu}$, $\eta_{r}$ often increases abruptly by order
of magnitude at a critical shear rate $\dot{\gamma}_{c}$. This phenomenon is
termed as Discontinuous shear-thickening (DST) Morris (2020). At volume
fractions $\phi\geq\phi_{J}^{\mu}$, the suspension flows smoothly at low
stresses but can be driven into a reversible solid-like jammed state. This
phenomenon is called shear jamming (SJ) Peters  _et al._ (2016); Singh  _et
al._ (2019); Seto  _et al._ (2019); Singh  _et al._ (2020a). The recent body
of work Guy  _et al._ (2018); Singh  _et al._ (2020a); Pradeep  _et al._
(2021); Singh  _et al._ (2022); Singh (2022) has linked both CST and DST to
the activation of constraints that hinder relative particle motion as stress
$\sigma$ is increased. These stress-activated constraints might arise from the
Coulombic, static friction Lin  _et al._ (2015); Guy  _et al._ (2015); Clavaud
_et al._ (2017); Hsu  _et al._ (2018); Seto  _et al._ (2013); Mari  _et al._
(2014); Ness and Sun (2016); Kawasaki and Berthier (2018); Singh  _et al._
(2020a, 2018); Mari and Seto (2019) or from combined effect of hydrodynamics
and asperities Jamali and Brady (2019). At low-applied stresses below
$\sigma^{*}$, steric repulsion keeps particle surfaces apart, maintaining the
lubrication layer between particles and offering no constraints in the
tangential pairwise motion. Above $\sigma^{*}$, the stabilization is overcome,
leading to the formation of direct, frictional contacts Seto  _et al._ (2013);
Mari  _et al._ (2014, 2015a); Ness and Sun (2016); Clavaud  _et al._ (2017);
Singh  _et al._ (2020a). The mean-field models Wyart and Cates (2014); Singh
_et al._ (2018); Guy  _et al._ (2018); Singh  _et al._ (2019, 2020a) capture
this transition between lubricated, frictionless “unconstrained” state to
frictional “constrained” state predict CST at low volume fractions $\phi$, DST
at volume fractions $\phi_{C}<\phi<\phi_{J}^{\mu}$, and shear jamming (SJ) at
$\phi\geq\phi_{J}^{\mu}$.

Stress-activated constraints on relative tangential motion between particles
due to static friction lead to stabilization of load bearing force network
under external deformation Mari  _et al._ (2014); Thomas  _et al._ (2018);
Boromand  _et al._ (2018); Singh  _et al._ (2020a); E. Thomas  _et al._
(2020); Sedes  _et al._ (2020); Gameiro  _et al._ (2020); Xu  _et al._ (2020);
Edens  _et al._ (2021); Nabizadeh  _et al._ (2022). With an increase in
$\sigma$, the number of frictional contacts increases, leading to particle
motion deviating from a pure affine trajectory to pass each other under
applied shear. This implies that the particle trajectories are more correlated
at higher stresses suggesting an enhanced collective motion of particle
groups, as reported previously Ness and Sun (2016); Singh  _et al._ (2020a).
In the case of frictionless dry granular particles, dissipation due to
non–affine velocity fluctuations have been shown to be connected to the
emergence of the jamming behavior Andreotti  _et al._ (2012). As jamming
approaches from below, the rearrangement of marginally rigid, isostatic
clusters, which eventually give way to long-range elastic interactions,
dominates the flow. The non-Brownian suspensions under shear present a
particular case where these ideas can be applied to a continuously evolving
high-viscosity fluid whose vicinity to the isostatic condition can be tuned by
two independent variables: $\phi$ and $\sigma$. Here we connect the concepts
developed in dry granular systems to the rheology of sheared dense non-
Brownian frictional suspensions to derive the relationship between rheology
and microscopic variables, namely diffusion of particles and non-affine
velocity correlation length or so-called rigid clusters.

In this paper, we numerically study the rheology and diffusion of frictional
non-Brownian spheres under simple shear. We investigate the microstructure and
dynamics as the suspension shear thickens; our specific aim is to extract the
relations among viscosity, diffusivity, and correlation length of the particle
motion. In the following, we explain our numerical methods (Sec. II), show our
results of viscosity, diffusivity, and spatial correlations of the particles
(Sec. III), and discuss our findings and future work (Sec. IV).

##  II Numerical methods

We perform quasi-2D simple shear simulations with a monolayer of $N=2000$ non-
Brownian spheres immersed in a Newtonian fluid. To avoid crystallization of
the system, we randomly distribute an equal volume of small and large
particles (with radii $a_{S}=a$ and $a_{L}=1.4a$) in an $L\times L$ periodic
square box using Lees-Edwards periodic boundary condition Lees and Edwards
(1972); O’Hern  _et al._ (2003). In the limit of zero Reynolds and Stokes
numbers, the problem reduces to the overdamped limit and force balance between
hydrodynamic $\mathrm{\mathbf{F}_{H}}$ and contact $\mathrm{\mathbf{F}_{C}}$
forces on every particle

$0=\mathrm{\mathbf{F}_{H}}(\mathbf{r,u})+\mathrm{\mathbf{F}_{C}}(\mathbf{r}),$
(1)

where $\mathbf{r}$ and $\mathbf{u}$ are many body particle positions and
velocities, respectively with $\mathbf{u}\equiv\dot{\mathbf{r}}$. Similar
relation also holds for the torque balance. We regularize the lubrication
singularity at a surface separation $h=10^{-3}a$ to allow the contact between
particles Mari  _et al._ (2014). The contact between particles is modeled
using the approach of Cundall and Strack Cundall and Strack (1979) by
introducing normal and tangential springs with spring constants being $k_{n}$
and $k_{t}=0.5k_{n}$, respectively. We tune stiffness at every $\phi$ such
that the maximum overlap does not exceed 3% of the particle radius to stay
near the rigid limit Singh  _et al._ (2015). To introduce rate dependence, a
critical load model (CLM) scheme, where a critical normal force $F^{*}$ is
needed to activate friction between particles Mari  _et al._ (2014). Here, we
use interparticle friction coefficient $\mu=1$. Note that we examine the
effects of the system size $N$ and friction coefficient $\mu$ on our numerical
results in Appendix B.

To simulate simple shear flows of the system, we apply external shear stress
$\sigma$ to the system under the Lees-Edwards boundary condition Lees and
Edwards (1972). The inset to Fig. 1a displays a snapshot of our simulation,
where the system is sheared along the horizontal arrows. At constant shear
stress, the shear rate becomes the fluctuating time-dependent quantity,

$\dot{\gamma}=\frac{\sigma-\sigma_{\mathrm{C}}}{\eta_{0}(1+2.5\phi)+\eta_{\mathrm{H}}}~{},$
(2)

where $\sigma_{\mathrm{C}}$ and $\sigma_{\mathrm{H}}$ are the contact and
hydrodynamic contributions to the total stress Mari  _et al._ (2015b); Seto
_et al._ (2019). All the length scales are normalized by small particle radius
$a$. The unit scale for shear rate and stress are $\dot{\gamma}_{0}\equiv
F^{\ast}/6\pi\eta_{0}a^{2}$ and $\sigma_{0}\equiv\eta_{0}\dot{\gamma}_{0}$,
respectively. The time scale in the simulations is $t_{0}=6\pi
a\eta_{0}/F^{\ast}$. Here $F^{\ast}$ is the minimum force needed to activate
frictional contacts leading to shear-thickening behavior.

From the simulations, the particle positions, normal and tangential contact
forces, and non-contact lubrication forces are extracted. At any given
condition, the simulations are performed for $40$ strain units, and the data
are recorded in strain steps of $10^{-1}$.

##  III Results

This section shows our numerical results of the quasi-2D simple shear
simulations. First, we focus on transverse motions of the non-Brownian
particles under shear, where their diffusivity and correlation length are
extracted from the transverse component of _mean squared displacement_ (MSD)
and spatial correlation function in transverse velocities, respectively (Sec.
III.1). Further, we examine how the diffusivity and correlation length are
affected by increase in volume fraction and stress (Sec. III.2). Then, we
relate the viscosity to diffusivity by scaling arguments (Sec. III.3), where
the scaling relations are robust to the changes of friction coefficient $\mu$
and system size $N$ (Appendix B).

###  III.1 Transverse motions

We analyze diffusion of the particles under shear by the transverse component
of MSD Olsson (2010); Kharel and Rognon (2017); Lemaître and Caroli (2009);
Chattoraj  _et al._ (2011),

$\Delta(\tau)^{2}=\left\langle\frac{1}{N}\sum_{i=1}^{N}\Delta
y_{i}(\tau)^{2}\right\rangle~{}.$ (3)

Here, $\Delta y_{i}(\tau)$ is the $y$-component of particle displacement for
the duration $\tau$ and the ensemble average $\langle\dots\rangle$ is taken
over different choices of the initial time 111The MSDs defined by the _total_
non-affine displacements show quantitatively the same results (data are not
shown for brevity).. Figure 1a displays the typical transverse MSDs,
$\Delta(\tau)^{2}$, where the horizontal axis is the duration scaled by shear
rate, $\gamma\equiv\dot{\gamma}\tau$ i.e., the shear strain applied to the
system for the duration $\tau$. As can be seen, the MSD exhibits a crossover
from super-diffusion (in the sense that the slope is larger than unity) to
normal diffusion (dashed line) around $\gamma\approx 1$. This trend is
qualitatively insensitive to the packing fraction $\phi$ and stress $\sigma$
(though the MSD increases with the increase of $\phi$, as listed in the
legend). Note that, different from thermal systems under shear Ikeda  _et al._
(2012); Mandal  _et al._ (2013); Yamamoto and Onuki (1998), we do not observe
any plateaus in the MSDs. Therefore, neither _caging_ nor _sub-diffusion_ of
the particles exists in our system Hatano (2011); Heussinger  _et al._ (2010);
Tsamados (2010).

We further introduce a correlation length as a typical size of _rigid
clusters_ Reichhardt and Reichhardt (2014); Olsson and Teitel (2007);
Heussinger and Barrat (2009a); Kharel and Rognon (2017); Saitoh and Mizuno
(2016); Saitoh and Kawasaki (2020), where collective motions of the particles
are interpreted as a result of rigid body rotations of the clusters. Following
Ref. Olsson and Teitel (2007); Tighe  _et al._ (2010); Heussinger and Barrat
(2009b), we quantify the collective motions by a spatial correlation function,

$C(x)=\langle v_{y}(x_{i},y_{i})v_{y}(x_{i}+x,y_{i})\rangle~{},$ (4)

where $v_{y}(x,y)$ is the transverse velocity field, and the ensemble average
$\langle\dots\rangle$ is taken over particle positions and time (in the steady
state). Figure 1b shows the correlation functions, where the horizontal axis
($x$-axis) is scaled by the small particle radius $a$. As can be seen, the
correlation function exhibits a well-defined minimum at a characteristic
distance, $x=\xi$ (as indicated by the vertical arrow for the data of
$\phi=0.80$). Because the minimum is negative, $C(\xi)<0$, the transverse
velocities are most “anti-correlated” at the minima $x=\xi$ Olsson and Teitel
(2007). Assuming that the rigid cluster is circular, its mean diameter is
comparable in size with $\xi$ 222 In a short range of $x$, transverse
velocities tend to be aligned in the same direction, i.e.,
$v_{y}(x_{i},y_{i})v_{y}(x_{i}+x,y_{i})>0$, such that the correlation function
$C(x)$ is a positive decreasing function of $x$. If the transverse velocities
align in the opposite direction, i.e.
$v_{y}(x_{i},y_{i})v_{y}(x_{i}+x,y_{i})<0$, the correlation function becomes
negative, $C(x)<0$. Because the correlation function is minimum when the
transverse velocities are located on either side of a vortex-like structure,
typical size of rigid clusters can be defined as the distance at which $C(x)$
becomes minimum. In a long distance of $x$, due to the randomness of
transverse velocities, the correlation function eventually decays to zero
Saitoh and Mizuno (2016).. Note that this “rigid cluster” is different from
the rigid clusters extracted using various algorithms like pebble game on the
contact network Henkes  _et al._ (2016); Zhang  _et al._ (2019).

Figure 1:  (a) Transverse MSDs, $\Delta(\tau)^{2}$ (Eq. (3)), as functions of
the shear strain $\gamma\equiv\dot{\gamma}\tau$, where the dashed line
indicates the normal diffusion, $\Delta^{2}\sim\gamma$. The packing fraction
$\phi$ increases as listed in the legend. Inset: A snapshot of our MD
simulation, where the particles (circles) are sheared along the horizontal
arrows, and the lines represent force-chain networks. (b) Spatial correlation
functions of the transverse velocities $C(x)$ (Eq. (4)). The vertical arrow
indicates the minimum of $C(x)$ for the data of $\phi=0.80$. The symbols are
as in (a).

###  III.2 Viscosity, diffusivity, and size of rigid clusters

Figure 2a displays $\eta_{r}$ as functions of the shear rate $\dot{\gamma}$
for packing fraction $\phi$ as listed in the legend. For $\phi\leq 0.77$, the
shear-thickening is continuous between two rate-independent values, i.e., CST
is observed. The shear thickening becomes discontinuous for
$\phi\rightarrow\phi_{c}\simeq 0.78$ and $\eta_{r}$ exhibits the
characteristic “S-shape” for $\phi\geq 0.78$. For $\phi=0.80$, the system at
the highest stress, $\sigma/\sigma_{0}=100$, is eventually _shear jammed_ such
that the viscosity cannot be defined. Note that the simulations presented here
are performed under fixed shear stress $\sigma$, while the S-shape indicates a
discontinuous jump (or drop) of $\eta_{r}$ in rate-controlled experiments
Boersma  _et al._ (1991), and simulations Mari  _et al._ (2015b).

As shown in Fig. 1a, the diffusion of the particles becomes normal diffusion,
where $\Delta^{2}\sim\gamma$ (dashed line) for $\gamma>1$. To quantify the
normal diffusion of the particles, we introduce diffusivity as 333We define
the diffusivity [Eq. (5)] as the slope of the MSD [Eq. (3)] in the normal
diffusive regime $\gamma=\dot{\gamma}\tau>1$, where we take sample averages of
$\Delta(\tau)^{2}/2\tau$ as $D/\dot{\gamma}\equiv<\Delta(\gamma)^{2}/2\gamma>$
in the range between $1<\gamma<40$.

$D=\lim_{\tau\rightarrow\infty}\frac{\Delta(\tau)^{2}}{2\tau}~{}.$ (5)

Figure 2b shows diffusivity scaled by the shear rate, $D/\dot{\gamma}$, as
functions of the shear rate $\dot{\gamma}$. The diffusivity of particles
exhibits behavior across the thickening transition analogous to rheology. We
observe a rate-independent plateau in the frictionless and frictional limits,
with increased diffusivity as the suspension shear thickens. Note that
diffusivity $D/\dot{\gamma}$ cannot be defined at the highest stress
($\sigma/\sigma_{0}=100$) if $\phi=0.80$, where the system is shear jammed so
that no particle can diffuse (in the transverse $y$-direction).

To further probe particle scale dynamics, we analyze the typical size of rigid
clusters $\xi$ using the correlation function, $C(x)$. Figure 2c plots $\xi$
as functions of scaled shear rate $\dot{\gamma}/\dot{\gamma}_{0}$ for various
values of $\phi$. Again, its dependence on $\phi$ and $\dot{\gamma}$ is quite
similar to those of $\eta_{r}$ and $D/\dot{\gamma}$ (Figs. 2a and b): $\xi$
increases as the suspension shear thickens. For $\phi=0.80$ the system is
shear jammed (and thus does not induce collective motions of the particles)
Kharel and Rognon (2017, 2018) at the highest stress
($\sigma/\sigma_{0}=100$).

In Appendix A, we analyze the _coordination number_ $z$ and _anisotropic
fabric_ $A$ of the particles Ness and Sun (2016). The coordination number $z$
increase with increasing stress. At the same time, the contact network becomes
less anisotropic as the suspension shear thickens (see Appendix A). These
observations are consistent with the findings of Ness & Sun Ness and Sun
(2016).

Figure 2:  Double logarithmic plots of (a) the relative viscosity $\eta_{r}$,
(b) diffusivity (Eq. (5)) over the shear rate $D/\dot{\gamma}$, and (c) size
of rigid clusters $\xi$ as functions of the shear rate $\dot{\gamma}$.
Different symbols represent different values of the packing fraction $\phi$ as
listed in the legend of (a).

###  III.3 Scaling arguments

In the previous section, we have shown that the dependence of diffusivity of
particles $D/\dot{\gamma}$ on $\phi$ and $\dot{\gamma}/\dot{\gamma}_{0}$ is
similar to that of $\eta_{r}(\phi,\dot{\gamma}/\dot{\gamma}_{0})$. Motivated
by this, below we derive relation between the rheology and diffusion of the
non-Brownian particles, which is reminiscent of the Stokes-Einstein relation
for Brownian particles. To reveal such a relation, we connect the viscosity
$\eta_{r}$ with diffusivity $D$ based on the scaling arguments in Refs.
Heussinger (2013); Maiti  _et al._ (2016); Kharel and Rognon (2017, 2018);
Saitoh and Kawasaki (2020).

Let us represent the transverse speed of the particles by characteristic
length and time scales as Heussinger (2013); Maiti  _et al._ (2016)

$v_{y}\sim\xi\dot{\gamma}~{}.$ (6)

To examine Eq. (6), we plot $\langle v_{y}^{2}\rangle$ as a function of
$\xi\dot{\gamma}$, where

$\langle
v_{y}^{2}\rangle=\left\langle\frac{1}{N}\sum_{i=1}^{N}v_{iy}^{2}\right\rangle$
(7)

is the mean squared transverse speed ($\langle\dots\rangle$ is the time
average in a steady state). Figure 3 confirms the scaling, Eq. (6), over the
whole ranges of $\phi$ and $\dot{\gamma}$, where the dashed line indicates
$\langle v_{y}^{2}\rangle\sim(\xi\dot{\gamma})^{2}$.

To extract the dependence of $\eta_{r}$ on the length scale $\xi$ we use the
energy balance equation Heussinger (2013); Maiti  _et al._ (2016). As the
system reaches a steady state, the injection of energy by shear
$\dot{\gamma}\sigma$ and the rate of energy dissipation $\Gamma$ (due to
dissipative forces) needs to be in balance as $\dot{\gamma}\sigma\sim\Gamma$
Woldhuis  _et al._ (2015). If the shear stress (or equivalently shear rate) is
low enough, the system is in the frictionless state where the drag force is
more significant for energy dissipation than the contact damping so that
$\Gamma\approx\Gamma_{\mathrm{visc}}\sim v_{y}^{2}$ 444 The dissipation rate
due to the drag force is proportional to the square of velocity fluctuation,
$\Gamma_{\mathrm{visc}}\propto(\delta v)^{2}$ Maiti  _et al._ (2016). Because
the mean flow is parallel to the $x$-axis, the velocity fluctuation is
comparable in size with the transverse speed, $\delta v\sim v_{y}$, and thus
$\Gamma\approx\Gamma_{\mathrm{visc}}\sim v_{y}^{2}$.. At higher stresses, as
the frictional contacts form, interparticle friction becomes the dominant
source of energy dissipation, where $\Gamma\approx\Gamma_{\mathrm{fric}}\sim
v_{y}^{3}/\dot{\gamma}$ 555 The dissipation rate due to the friction scales as
$\Gamma_{\mathrm{fric}}\propto\xi\delta v$ Maiti  _et al._ (2016)..
Substituting Eq. (6) into $\Gamma_{\mathrm{visc}}$ and
$\Gamma_{\mathrm{fric}}$, we find scaling relations of the dissipation rate as
Maiti  _et al._ (2016)

$\Gamma\sim\begin{cases}\xi^{2}\dot{\gamma}^{2}&(\mathrm{{frictionless}})\\\
\xi^{3}\dot{\gamma}^{2}&(\mathrm{{frictional}})\end{cases}~{}.$ (8)

If we rewrite the energy injection as
$\dot{\gamma}\sigma=\dot{\gamma}^{2}\eta_{r}$, the energy balance
($\dot{\gamma}\sigma\sim\Gamma$) tells us
$\eta_{r}\sim\Gamma/\dot{\gamma}^{2}$ so that the viscosity depends on the
length scale as

$\eta_{r}\sim\begin{cases}\xi^{2}&(\mathrm{{frictionless}})\\\
\xi^{3}&(\mathrm{{frictional}})\end{cases}~{}.$ (9)

In addition, it was recently pointed out that rotations of rigid clusters
govern the diffusion of the particles under shear, where the size of rigid
clusters describes the diffusivity as $D\sim a\xi\dot{\gamma}$ Kharel and
Rognon (2017, 2018); Saitoh and Kawasaki (2020). At low stresses, the
suspension is in the frictionless state; the relative particle motion is
correlated to only a few particles radii (Fig. 2c), i.e. $\xi\sim a$, meaning
that the diffusivity is simply given by $D\sim\xi^{2}\dot{\gamma}$. As
frictional contacts form with increasing stress, the correlation length can be
as large as $\xi\sim 50a$ (Fig. 2c). This presents a length scale associated
with particle dynamics that is distinct from the particle size, i.e. $\xi\gg
a$. Additionally diffusivity scales as $D\sim\xi\dot{\gamma}$ Kharel and
Rognon (2017, 2018); Saitoh and Kawasaki (2020). To summarize, scaling
relations of the diffusivity are given by

$D\sim\begin{cases}\xi^{2}\dot{\gamma}&(\mathrm{{frictionless}})\\\
\xi\dot{\gamma}&(\mathrm{{frictional}})\end{cases}~{}.$ (10)

From Eqs. (9) and (10), we can derive a connection between the viscosity and
diffusivity as

$D/\dot{\gamma}\sim\begin{cases}\eta_{r}&(\mathrm{{frictionless}})\\\
\eta_{r}^{1/3}&(\mathrm{{frictional}})\end{cases}~{}.$ (11)

To validate Eq. (11), we plot $D/\dot{\gamma}$ vs $\eta_{r}$ in Fig. 4, where
the dashed and solid lines indicate the relations,
$D/\dot{\gamma}\sim\eta_{r}$ and $\eta_{r}^{1/3}$, respectively.

In Appendix B, we show that one-to-one relations, Eq. (11), are robust to the
changes of friction coefficient $\mu$ and system size $N$. Therefore, it is
universal that the viscosity and diffusivity of the non-Brownian particles
under shear are governed by the size of rigid clusters and are related with
each other. Because the one-to-one relations are insensitive to the system
size, our results are distinct from numerical results of sheared amorphous
solids Lemaître and Caroli (2009); Chattoraj  _et al._ (2011), where the
diffusivity over the shear rate becomes linear in the system length as
$D/\dot{\gamma}\sim L$ in the quasi-static limit. In Refs. Lemaître and Caroli
(2009); Chattoraj  _et al._ (2011), the interparticle friction is not
introduced so that S-shaped transitions of $\eta_{r}$ and $D/\dot{\gamma}$
(from the frictionless to frictional states) are not expected. Thus, the one-
to-one relations, Eq. (11), are specific to our frictional non-Brownian
particles.

Figure 3:  Scatter plots of $\langle v_{y}^{2}\rangle$ and $\xi\dot{\gamma}$,
where the dashed line represents the scaling, $\langle
v_{y}^{2}\rangle\sim(\xi\dot{\gamma})^{2}$. The packing fraction $\phi$
increases as listed in the legend.  Figure 4:  Scatter plots of the
diffusivity over the shear rate $D/\dot{\gamma}$ and relative viscosity
$\eta_{r}$. The dashed line indicates the scaling in the frictionless state,
i.e. $D/\dot{\gamma}\sim\eta_{r}$, while the solid line is that in the
frictional state, i.e. $D/\dot{\gamma}\sim\eta_{r}^{1/3}$. The packing
fraction $\phi$ increases as listed in the legend.

##  IV Discussion

In this study, we have numerically investigated the rheology, diffusion, and
spatial correlations of rigid frictional non-Brownian suspension under shear.
We analyzed the transverse component of MSD and spatial correlation function
of the transverse velocities to extract the diffusivity over the shear rate
$D/\dot{\gamma}$ and size of rigid clusters $\xi$, respectively. We found that
the dependence of $D/\dot{\gamma}$ and $\xi$ on the packing fraction $\phi$
and shear rate $\dot{\gamma}$ are analogous to that of the relative viscosity
$\eta_{r}(\dot{\gamma})$. They tend to plateau at low and high shear rates
(stresses) and increase between these plateaus as the suspension shear
thickens. We have derived one-to-one relations between rheology and diffusion
using extensive numerical data for rheology, diffusivity, and correlation
lengths. This, in spirit, is similar to the Stokes-Einstein relation for
Brownian particles. Additionally, making use of the scaling arguments
suggested in Refs. Heussinger (2013); Maiti  _et al._ (2016), we introduced
the power-law scaling $\eta_{r}\sim\xi^{\alpha}$ with the exponents,
$\alpha=2$ and $3$, for the frictionless and frictional states, respectively.
We also explained the scaling relations of the diffusivity as
$D\sim\xi^{2}\dot{\gamma}$ (frictionless) and $\xi\dot{\gamma}$ (frictional).
Then, from these scaling relations, we have found the one-to-one relations as
$D/\dot{\gamma}\sim\eta_{r}$ and $\eta_{r}^{1/3}$ for the frictionless and
frictional states, respectively. These relations were numerically validated
and found to be robust to the changes of the friction coefficient $\mu$ and
system size $N$. It is surprising that the one-to-one relations, Eq. (11), are
not altered by the changes of $\mu$ and $N$ because the shear-thickening is
sensitive to the interparticle friction Seto  _et al._ (2013) and the
diffusivity is believed to be system-size dependent Lemaître and Caroli
(2009); Chattoraj  _et al._ (2011). Our results suggest that the typical size
of rigid clusters $\xi$ exhibits a similar dependence on the parameters
($\phi$, $\mu$, and $N$) with $\eta_{r}$ and $D/\dot{\gamma}$, and the
viscosity and diffusivity of frictional non-Brownian particles are governed by
the length scale $\xi$.

###  IV.1 Comparison with previous studies

In Refs. Heussinger (2013); Maiti  _et al._ (2016), the scaling relations,
Eqs. (6) and (8), have been numerically confirmed for the 2D system of
inertial suspensions under simple shear. We also numerically validated Eq. (6)
and used Eq. (8) to derive the one-to-one relations, Eq. (11). The main
difference between our model and that in Refs. Heussinger (2013); Maiti  _et
al._ (2016) is the absence (or presence) of particle inertia. Therefore, the
scaling relations, Eqs. (6) and (8), seem to be independent of the particle
mass though it is left to future work to examine whether the one-to-one
relations, Eq. (11), are altered by the particle inertia.

In Ref. Ness and Sun (2016), the sudden increases of $D/\dot{\gamma}$, $\xi$,
coordination number $z$, and anisotropic fabric $A$ around DST have already
been reported. However, the connection between rheology and diffusion was not
investigated, and the effects of friction coefficient $\mu$ and system size
$N$ were not studied. In addition, it was found in molecular dynamics
simulations of sheared amorphous solids Lemaître and Caroli (2009); Chattoraj
_et al._ (2011) that the diffusivity over the shear rate becomes linear in the
system length as $D/\dot{\gamma}\sim L$ if the packing fraction is far above
the jamming transition density and the shear rate is sufficiently small.
Different from the numerical model used in Refs. Lemaître and Caroli (2009);
Chattoraj  _et al._ (2011), we use the frictional non-Brownian particles and
the correlation length $\xi$ does not extend to the system length, i.e.,
$\xi<L$. As a result, we did not observe any finite-size effects in the
diffusivity.

###  IV.2 Future works

In our simulations, we employed the model for shear-thickening suspensions to
investigate the connection between their rheology and diffusion under shear.
However, another drastic change of the viscosity is known as _discontinuous
shear thinning_ Irani  _et al._ (2019), where cohesive forces are introduced
between the frictionless particles. Therefore, it is important to analyze the
relation between rheology and diffusion of soft cohesive particles under shear
and clarify whether the diffusion is “suppressed” by discontinuous shear
thinning or if cohesion would lead to shear localization Singh  _et al._
(2014). In addition, recent studies of the rheology of non-spherical particles
Nath and Heussinger (2019); Marschall  _et al._ (2019); Keta and Olsson (2020)
suggest that the influence of particle shapes on one-to-one relations should
be examined in the future. Moreover, extending the current work to three
dimensions Olsson (2019); Oyama  _et al._ (2019) is crucial to the practical
applications of this work. Confinement that can slow down the dynamics Li  _et
al._ (2020); Singh  _et al._ (2020b) is another essential aspect to
investigate.

###### Acknowledgements.

We thank A. Ikeda, H. Mizuno, and S. Yamanaka for fruitful discussions. A. S.
would like to acknowledge Case Western Reserve University for the start-up
funding for the project. K.S. was financially supported by JSPS KAKENHI from
JSPS (Grant Numbers 18K13464, 20H01868, and 21H01006) and Inamori Research
Grants 2021.

## References

  * Coussot (1997) P. Coussot,  _Mudflow Rheology and Dynamics_ (CRC Press, 1997).
  * Jerolmack and Daniels (2019) D. J. Jerolmack and K. E. Daniels, Nature Reviews Physics 1, 716 (2019).
  * Benbow and Bridgewater (1993) J. Benbow and J. Bridgewater,  _Paste flow and extrusion_ (Oxford University Press(UK), 1993).
  * Maheshwari  _et al._ (2019) A. J. Maheshwari, A. M. Sunol, E. Gonzalez, D. Endy, and R. N. Zia, Physical Review Fluids 4, 110506 (2019).
  * Morris (2020) J. F. Morris, Annual Review of Fluid Mechanics 52, 121 (2020).
  * Denn and Morris (2014) M. M. Denn and J. F. Morris, Annu. Rev. Chem. Biomol. Eng. 5 (2014).
  * Guazzelli and Pouliquen (2018) É. Guazzelli and O. Pouliquen, J. Fluid Mech. 852 (2018).
  * Ness  _et al._ (2022) C. Ness, R. Seto, and R. Mari, Annual Review of Condensed Matter Physics 13, 97 (2022).
  * Jamali  _et al._ (2020) S. Jamali, E. Del Gado, and J. F. Morris, Journal of Rheology 64, 1501 (2020).
  * Bonn  _et al._ (2017) D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Rev. Mod. Phys. 89, 035005 (2017).
  * Torquato  _et al._ (2000) S. Torquato, T. M. Truskett, and P. G. Debenedetti, Phys. Rev. Lett. 84, 2064 (2000).
  * Makse  _et al._ (2000) H. A. Makse, D. L. Johnson, and L. M. Schwartz, Phys. Rev. Lett. 84, 4160 (2000).
  * Silbert  _et al._ (2002) L. E. Silbert, D. Ertas, G. S. Grest, T. C. Halsey, and D. Levine, Phys. Rev. E 65, 031304 (2002).
  * Zhang and Makse (2005) H. P. Zhang and H. A. Makse, Phys. Rev. E 72, 011301 (2005).
  * Shundyak  _et al._ (2007) K. Shundyak, M. van Hecke, and W. van Saarloos, Phys. Rev. E 75, 010301(R) (2007).
  * Singh  _et al._ (2020a) A. Singh, C. Ness, R. Seto, J. J. de Pablo, and H. M. Jaeger, Phys. Rev. Lett. 124, 248005 (2020a).
  * Santos  _et al._ (2020) A. P. Santos, D. S. Bolintineanu, G. S. Grest, J. B. Lechman, S. J. Plimpton, I. Srivastava, and L. E. Silbert, Physical Review E 102, 032903 (2020).
  * Pradeep  _et al._ (2021) S. Pradeep, M. Nabizadeh, A. R. Jacob, S. Jamali, and L. C. Hsiao, Physical Review Letters 127, 158002 (2021).
  * Singh  _et al._ (2022) A. Singh, G. L. Jackson, M. van der Naald, J. J. de Pablo, and H. M. Jaeger, Physical Review Fluids 7, 054302 (2022).
  * Singh (2022) A. Singh, Science Talks 3, 100028 (2022).
  * James  _et al._ (2018) N. M. James, E. Han, R. A. L. de la Cruz, J. Jureller, and H. M. Jaeger, Nat. Mater. 17, 965 (2018).
  * James  _et al._ (2019) N. M. James, C.-P. Hsu, N. D. Spencer, H. M. Jaeger, and L. Isa, J. Phys. Chem. Lett. 10, 1663 (2019).
  * Donev  _et al._ (2004) A. Donev, I. Cisse, D. Sachs, E. A. Variano, F. H. Stillinger, R. Connelly, S. Torquato, and P. M. Chaikin, Science 303, 990 (2004).
  * Salerno  _et al._ (2018) K. M. Salerno, D. S. Bolintineanu, G. S. Grest, J. B. Lechman, S. J. Plimpton, I. Srivastava, and L. E. Silbert, Physical Review E 98, 050901 (2018).
  * Hopkins  _et al._ (2013) A. B. Hopkins, F. H. Stillinger, and S. Torquato, Physical Review E 88, 022205 (2013).
  * Guy  _et al._ (2020) B. M. Guy, C. Ness, M. Hermes, L. J. Sawiak, J. Sun, and W. C. Poon, Soft matter 16, 229 (2020).
  * Srivastava  _et al._ (2021) I. Srivastava, S. A. Roberts, J. T. Clemmer, L. E. Silbert, J. B. Lechman, and G. S. Grest, Physical Review Research 3, L032042 (2021).
  * Fall  _et al._ (2008) A. Fall, N. Huang, F. Bertrand, G. Ovarlez, and D. Bonn, Phys. Rev. Lett. 100, 018301 (2008).
  * Brown and Jaeger (2009) E. Brown and H. M. Jaeger, Phys. Rev. Lett. 103, 086001 (2009).
  * Peters  _et al._ (2016) I. R. Peters, S. Majumdar, and H. M. Jaeger, Nature 532, 214 (2016).
  * Singh  _et al._ (2019) A. Singh, S. Pednekar, J. Chun, M. M. Denn, and J. F. Morris, Phys. Rev. Lett. 122, 098004 (2019).
  * Seto  _et al._ (2019) R. Seto, A. Singh, B. Chakraborty, M. M. Denn, and J. F. Morris, Gran. Matt. 21, 82 (2019).
  * Guy  _et al._ (2018) B. M. Guy, J. Richards, D. Hodgson, E. Blanco, and W. C. K. Poon, Phys. Rev. Lett. 121, 128001 (2018).
  * Lin  _et al._ (2015) N. Y. C. Lin, B. M. Guy, M. Hermes, C. Ness, J. Sun, W. C. K. Poon, and I. Cohen, Phys. Rev. Lett. 115, 228304 (2015).
  * Guy  _et al._ (2015) B. M. Guy, M. Hermes, and W. C. K. Poon, Phys. Rev. Lett. 115, 088304 (2015).
  * Clavaud  _et al._ (2017) C. Clavaud, A. Bérut, B. Metzger, and Y. Forterre, Proc. Natl. Acad. Sci. U.S.A. , 5147 (2017).
  * Hsu  _et al._ (2018) C.-P. Hsu, S. N. Ramakrishna, M. Zanini, N. D. Spencer, and L. Isa, Proc. Nat. Acad. Sci. (2018).
  * Seto  _et al._ (2013) R. Seto, R. Mari, J. F. Morris, and M. M. Denn, Phys. Rev. Lett. 111, 218301 (2013).
  * Mari  _et al._ (2014) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, J. Rheol. 58, 1693 (2014).
  * Ness and Sun (2016) C. Ness and J. Sun, Soft Matter 12, 914 (2016).
  * Kawasaki and Berthier (2018) T. Kawasaki and L. Berthier, Phys. Rev. E 98, 012609 (2018).
  * Singh  _et al._ (2018) A. Singh, R. Mari, M. M. Denn, and J. F. Morris, J. Rheol. 62, 457 (2018).
  * Mari and Seto (2019) R. Mari and R. Seto, Soft Matter 15, 6650 (2019).
  * Jamali and Brady (2019) S. Jamali and J. F. Brady, Physical review letters 123, 138002 (2019).
  * Mari  _et al._ (2015a) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Proc. Natl. Acad. Sci. U.S.A. 112, 15326 (2015a).
  * Wyart and Cates (2014) M. Wyart and M. E. Cates, Phys. Rev. Lett. 112, 098302 (2014).
  * Thomas  _et al._ (2018) J. E. Thomas, K. Ramola, A. Singh, R. Mari, J. F. Morris, and B. Chakraborty, Phys. Rev. Lett. 121, 128002 (2018).
  * Boromand  _et al._ (2018) A. Boromand, S. Jamali, B. Grove, and J. M. Maia, Journal of Rheology 62, 905 (2018).
  * E. Thomas  _et al._ (2020) J. E. Thomas, A. Goyal, D. Singh Bedi, A. Singh, E. Del Gado, and B. Chakraborty, Journal of Rheology 64, 329 (2020).
  * Sedes  _et al._ (2020) O. Sedes, A. Singh, and J. F. Morris, Journal of Rheology 64, 309 (2020).
  * Gameiro  _et al._ (2020) M. Gameiro, A. Singh, L. Kondic, K. Mischaikow, and J. F. Morris, Physical Review Fluids 5, 034307 (2020).
  * Xu  _et al._ (2020) Q. Xu, A. Singh, and H. M. Jaeger, Journal of Rheology 64, 321 (2020).
  * Edens  _et al._ (2021) L. E. Edens, E. G. Alvarado, A. Singh, J. F. Morris, G. K. Schenter, J. Chun, and A. E. Clark, Soft Matter 17, 7476 (2021).
  * Nabizadeh  _et al._ (2022) M. Nabizadeh, A. Singh, and S. Jamali, Physical Review Letters 129, 068001 (2022).
  * Andreotti  _et al._ (2012) B. Andreotti, J.-L. Barrat, and C. Heussinger, Phys. Rev. Lett. 109, 105901 (2012).
  * Lees and Edwards (1972) A. W. Lees and S. F. Edwards, J. Phys. C: Solid State Phys. 5, 1921 (1972).
  * O’Hern  _et al._ (2003) C. S. O’Hern, L. E. Silbert, A. J. Liu, and S. R. Nagel, Phys. Rev. E 68, 011306 (2003).
  * Cundall and Strack (1979) P. A. Cundall and O. D. L. Strack, Geotechnique 29, 47 (1979).
  * Singh  _et al._ (2015) A. Singh, V. Magnanimo, K. Saitoh, and S. Luding, New J. Phys 17, 043028 (2015).
  * Mari  _et al._ (2015b) R. Mari, R. Seto, J. F. Morris, and M. M. Denn, Phys. Rev. E 91, 052302 (2015b).
  * Olsson (2010) P. Olsson, Phys. Rev. E 81, 040301(R) (2010).
  * Kharel and Rognon (2017) P. Kharel and P. Rognon, Phys. Rev. Lett. 119, 178001 (2017).
  * Lemaître and Caroli (2009) A. Lemaître and C. Caroli, Phys. Rev. Lett. 103, 065501 (2009).
  * Chattoraj  _et al._ (2011) J. Chattoraj, C. Caroli, and A. Lemaître, Phys. Rev. E 84, 011501 (2011).
  * Note (1) The MSDs defined by the _total_ non-affine displacements show quantitatively the same results (data are not shown for brevity).
  * Ikeda  _et al._ (2012) A. Ikeda, L. Berthier, and P. Sollich, Phys. Rev. Lett. 109, 018301 (2012).
  * Mandal  _et al._ (2013) S. Mandal, V. Chikkadi, B. Nienhuis, D. Raabe, P. Schall, and F. Varnik, Phys. Rev. E 88, 022129 (2013).
  * Yamamoto and Onuki (1998) R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 (1998).
  * Hatano (2011) T. Hatano, J. Phys. Conf. Series 319, 012011 (2011).
  * Heussinger  _et al._ (2010) C. Heussinger, L. Berthier, and J.-L. Barrat, Eur. Phys. Lett. 90, 20005 (2010).
  * Tsamados (2010) M. Tsamados, Eur. Phys. J. E 32, 165 (2010).
  * Reichhardt and Reichhardt (2014) C. Reichhardt and C. J. O. Reichhardt, Soft Matter 10, 2932 (2014).
  * Olsson and Teitel (2007) P. Olsson and S. Teitel, Phys. Rev. Lett. 99, 178001 (2007).
  * Heussinger and Barrat (2009a) C. Heussinger and J.-L. Barrat, Phys. Rev. Lett. 102, 218303 (2009a).
  * Saitoh and Mizuno (2016) K. Saitoh and H. Mizuno, Soft Matter 12, 1360 (2016).
  * Saitoh and Kawasaki (2020) K. Saitoh and T. Kawasaki, Front. Phys. 8, 99 (2020).
  * Tighe  _et al._ (2010) B. P. Tighe, E. Woldhuis, J. J. C. Remmers, W. van Saarloos, and M. van Hecke, Phys. Rev. Lett. 105, 088303 (2010).
  * Heussinger and Barrat (2009b) C. Heussinger and J.-L. Barrat, Phys. Rev. Lett. 102, 218303 (2009b).
  * Note (2) In a short range of $x$, transverse velocities tend to be aligned in the same direction, i.e., $v_{y}(x_{i},y_{i})v_{y}(x_{i}+x,y_{i})>0$, such that the correlation function $C(x)$ is a positive decreasing function of $x$. If the transverse velocities align in the opposite direction, i.e. $v_{y}(x_{i},y_{i})v_{y}(x_{i}+x,y_{i})<0$, the correlation function becomes negative, $C(x)<0$. Because the correlation function is minimum when the transverse velocities are located on either side of a vortex-like structure, typical size of rigid clusters can be defined as the distance at which $C(x)$ becomes minimum. In a long distance of $x$, due to the randomness of transverse velocities, the correlation function eventually decays to zero Saitoh and Mizuno (2016).
  * Henkes  _et al._ (2016) S. Henkes, D. A. Quint, Y. Fily, and J. M. Schwarz, Phys. Rev. Lett. 116, 028301 (2016).
  * Zhang  _et al._ (2019) S. Zhang, L. Zhang, M. Bouzid, D. Z. Rocklin, E. Del Gado, and X. Mao, Physical review letters 123, 058001 (2019).
  * Boersma  _et al._ (1991) W. H. Boersma, P. J. M. Baets, J. Laven, and H. N. Stein, J. Rheol. 35, 1093 (1991).
  * Note (3) We define the diffusivity [Eq. (5)] as the slope of the MSD [Eq. (3)] in the normal diffusive regime $\gamma=\dot{\gamma}\tau>1$, where we take sample averages of $\Delta(\tau)^{2}/2\tau$ as $D/\dot{\gamma}\equiv<\Delta(\gamma)^{2}/2\gamma>$ in the range between $1<\gamma<40$.
  * Kharel and Rognon (2018) P. Kharel and P. Rognon, Eur. Phys. Lett. 124, 24002 (2018).
  * Heussinger (2013) C. Heussinger, Phys. Rev. E 88, 050201(R) (2013).
  * Maiti  _et al._ (2016) M. Maiti, A. Zippelius, and C. Heussinger, Eur. Phys. Lett. 115, 54006 (2016).
  * Woldhuis  _et al._ (2015) E. Woldhuis, V. Chikkadi, M. S. van Deen, P. Schall, and M. van Hecke, Soft Matter 11, 7024 (2015).
  * Note (4) The dissipation rate due to the drag force is proportional to the square of velocity fluctuation, $\Gamma_{\mathrm{visc}}\propto(\delta v)^{2}$ Maiti  _et al._ (2016). Because the mean flow is parallel to the $x$-axis, the velocity fluctuation is comparable in size with the transverse speed, $\delta v\sim v_{y}$, and thus $\Gamma\approx\Gamma_{\mathrm{visc}}\sim v_{y}^{2}$.
  * Note (5) The dissipation rate due to the friction scales as $\Gamma_{\mathrm{fric}}\propto\xi\delta v$ Maiti  _et al._ (2016).
  * Irani  _et al._ (2019) E. Irani, P. Chaudhuri, and C. Heussinger, Phys. Rev. Fluids 4, 074307 (2019).
  * Singh  _et al._ (2014) A. Singh, V. Magnanimo, K. Saitoh, and S. Luding, Physical Review E 90, 022202 (2014).
  * Nath and Heussinger (2019) T. Nath and C. Heussinger, Eur. Phys. J. E 42, 157 (2019).
  * Marschall  _et al._ (2019) T. Marschall, Y.-E. Keta, P. Olsson, and S. Teitel, Phys. Rev. Lett. 122, 188002 (2019).
  * Keta and Olsson (2020) Y.-E. Keta and P. Olsson, Phys. Rev. E 102, 052905 (2020).
  * Olsson (2019) P. Olsson, Phys. Rev. Lett. 122, 108003 (2019).
  * Oyama  _et al._ (2019) N. Oyama, H. Mizuno, and K. Saitoh, Phys. Rev. Lett. 122, 188004 (2019).
  * Li  _et al._ (2020) J. Li, X. Jiang, A. Singh, O. G. Heinonen, J. P. Hernández-Ortiz, and J. J. de Pablo, The Journal of Chemical Physics 152, 204109 (2020).
  * Singh  _et al._ (2020b) A. Singh, J. Li, X. Jiang, J. P. Hernández-Ortiz, H. M. Jaeger, and J. J. de Pablo, Physics of Fluids 32, 053307 (2020b).
  * Singh  _et al._ (2013) A. Singh, V. Magnanimo, and S. Luding, in  _AIP Conference Proceedings_ , Vol. 1542 (American Institute of Physics, 2013) pp. 682–685.

##  Appendix A Viscosity, diffusivity, and size of rigid clusters as functions
of the external shear stress

In this appendix, we show additional data of the relative viscosity
$\eta_{r}$, diffusivity over the shear rate $D/\dot{\gamma}$, and size of
rigid clusters $\xi$. Since we apply the external shear stress $\sigma$ to the
system, the shear rate $\dot{\gamma}$ is “measured” in simulations. We plot
$\eta_{r}$, $D/\dot{\gamma}$, and $\xi$ as functions of $\sigma$ which is the
control parameter. The evolution of _contact networks_ is presented using
coordination number $z$ and anisotropic fabric $A$ Ness and Sun (2016).

Figures 5a-c are double-logarithmic plots of (a) $\eta_{r}$, (b)
$D/\dot{\gamma}$, and (c) $\xi$ as functions of $\sigma$, where the packing
fraction $\phi$ increases as listed in the legend of (a). Increasing $\sigma$,
we observe that all the quantities start to increase around a characteristic
shear stress, $\sigma=\sigma_{0}=0.3$. In these figures, the data of the
densest packing ($\phi=0.80$) are truncated at $\sigma=5\sigma_{0}$ because
the system is shear-jammed, i.e. does not flow, for $\sigma>5\sigma_{0}$ so
that $\eta_{r}$, $D/\dot{\gamma}$, and $\xi$ cannot be defined. Figures 5d and
e are semi-logarithmic plots of (d) $z$ and (e) $A$ as functions of $\sigma$.
As can be seen, both $z$ and $A$ increase as the suspension shear-thickens,
implying that as frictional contacts form, the contact network rearranges and
becomes less anisotropic.

We also display $z$ and $A$ as functions of the shear rate $\dot{\gamma}$
(which is measured in simulations) in Figs. 6a and b, respectively.

Figure 5:  Double-logarithmic plots of (a) the relative viscosity $\eta_{r}$,
(b) diffusivity over the shear rate $D/\dot{\gamma}$, and (c) size of rigid
clusters $\xi$ as functions of the external shear stress $\sigma$. Semi-
logarithmic plots of (d) the coordination number $z$ and (e) anisotropic
fabric $A$ as functions of $\sigma$. Here, the packing fraction $\phi$
increases as listed in the legend of (a).  Figure 6:  Semi-logarithmic plots
of (a) $z$ and (b) $A$ as functions of the shear rate $\dot{\gamma}$, where
the symbols are as in Fig. 5.

##  Appendix B Dependence on the friction coefficient and system size

In this appendix, we examine the influence of the friction coefficient $\mu$
and system size (the number of particles) $N$ on the scaling relations, Eqs.
(6) and (11).

Figure 7 displays scatter plots of (a) $\langle v_{y}^{2}\rangle$ and
$\xi\dot{\gamma}$, and (b) $D/\dot{\gamma}$ and $\eta_{r}$, where the system
size is given by $N=2000$ and we vary $\mu$ as listed in the legend of (a).
Even though the friction coefficient affects the onset of DST, by controlling
the $\phi_{J}$ Singh  _et al._ (2013, 2018). As can be seen, the scaling
relations (dashed and solid lines) well describe the data over the whole
ranges of the packing fraction $\phi$ and shear rate $\dot{\gamma}$.
Therefore, Eqs. (6) and (11) are not affected by the strength of friction
coefficient $\mu$.

Figure 8 shows the scatter plots of (a) $\langle v_{y}^{2}\rangle$ and
$\xi\dot{\gamma}$, and (b) $D/\dot{\gamma}$ and $\eta_{r}$, where $\mu=1$ and
we increase $N$ as listed in the legend of (a). As can be seen, all the data
are well described by the scaling relations (dashed and solid lines) so that
Eqs. (6) and (11) hold even if we change the system size $N$.

Figure 7:  Scatter plots of (a) $\langle v_{y}^{2}\rangle$ and
$\xi\dot{\gamma}$, and (b) $D/\dot{\gamma}$ and $\eta_{r}$, where the system
size is $N=2000$ and the friction coefficient $\mu$ increases as listed in the
legend of (a) (“infinity” means $\mu=\infty$). The lines in (a) and (b) are as
in Figs. 3 and 4, respectively.  Figure 8:  Scatter plots of (a) $\langle
v_{y}^{2}\rangle$ and $\xi\dot{\gamma}$, and (b) $D/\dot{\gamma}$ and
$\eta_{r}$, where the friction coefficient is $\mu=1.0$ and the system size
$N$ increases as listed in the legend of (a). The lines in (a) and (b) are as
in Figs. 3 and 4, respectively.


# PyFAI: a Python library for high performance azimuthal integration on GPU

Jérôme Kieffer∗†, Giannis Ashiotis† * Corresponding author:
<EMAIL_ADDRESS>European Synchrotron Radiation Facility, Grenoble,
FranceCopyright © 2014 Jérôme Kieffer et al. This is an open-access article
distributed under the terms of the Creative Commons Attribution License, which
permits unrestricted use, distribution, and reproduction in any medium,
provided the original author and source are credited.
http://creativecommons.org/licenses/by/3.0/

###### Abstract

The pyFAI package has been designed to reduce X-ray diffraction images into
powder diffraction curves to be further processed by scientists. This
contribution describes how to convert an image into a radial profile using the
Numpy package, how the process was accelerated using Cython. The algorithm was
parallelised, needing a complete re-design to benefit from massively parallel
devices like graphical processing units or accelerators like the Intel Xeon
Phi using the PyOpenCL library.

###### Index Terms:

X-rays, powder diffraction, SAXS, HPC, parallel algorithms, OpenCL.

##  1 Introduction

The Python programming language is widely adopted in the scientific community
and especially in crystallography, this is why a convenient azimuthal
integration routine, one of the basic algorithms in crystallography, was
requested by the synchrotron community. The advent of pixel-detectors with
their very high speed (up to 3000 frames per second) imposed strong
constraints in speed that most available programs ([FIT2D], [SPD], …), written
in Fortran or C, could not meet.

The [pyFAI] project started in 2011 and aims at providing a convenient
_Pythonic_ interface for azimuthal integration, so that any diffractionist can
adapt it to the type of experiment he is interested in. This contribution
describes how one of the most fundamental algorithms used in crystallography
has been implemented in Python and how it was accelerated to match the readout
speeds of today’s fastest detectors.

After describing a typical experiment and explaining what is measured and how
it must be transformed (section 2), section 3 describes how the algorithm can
be vectorised using [NumPy] and sped up with [Cython]. Section 4 highlights
the accuracy enhancement recently introduced while section 5 focuses on the
parallelisation of the azimuthal integration task on many-core systems like
Graphical Processing Units (GPU) or on accelerators via [PyOpenCL]. In section
6, serial and parallel implementations using [OpenMP] and [OpenCL] from
various vendors and devices are benchmarked.

##  2 Description of the experiment

X-rays are electromagnetic waves, similar to visible light, except for their
wavelengths which are much shorter, typically of the size of inter-atomic
distances, making them a perfect probe to analyse atomic and molecular
structures. X-rays can be elastically scattered (i.e. re-emitted with the same
energy) by the electron cloud surrounding atoms. When atoms are arranged
periodically, as in a crystal, scattered X-rays interfere in a constructive
way when the difference of their optical paths is a multiple of the
wavelength: $2dsin(\theta)=n\lambda$. In this formula, known as _Bragg’s law_
, _d_ is the distance between crystal plans, $\theta$ is the incidence angle,
$\lambda$ is the wavelength and n is an integer. An X-ray beam crossing a
powder-like sample made of many randomly oriented small crystals is then
scattered along multiple concentric cones. In a powder diffraction experiment,
one aims at measuring the intensity of X-rays as a function of the cone
aperture, averaged along each ring. This transformation is called "azimuthal
integration" as it is an averaging of the signal along the azimuthal angle.

Figure 1: Debye-Scherrer cones obtained from diffraction of a monochromatic
X-Ray beam by a powder of crystallised material. (Credits: CC-BY-SA Klaus-
Dieter Liss)

##  3 Azimuthal integration

While pyFAI addresses the needs of both mono- and bi-dimensional integration
in different spaces (real or reciprocal), this contribution focuses on the
algorithmic and implementation part of the method. The coordinates in which
pyFAI operates can be one of the following:

  * •

$r=\sqrt{x^{2}+y^{2}}$

  * •

$\chi=tan^{-1}(y/x)$

  * •

$2\theta=tan^{-1}(r/d)$

  * •

$q=4\pi sin({2\theta}/2)/\lambda$

The pyFAI library was designed to feature a _Pythonic_ interface and work
together with [FabIO] for image reading (or [H5Py] for HDF5 files). The
following snippet of code explains the basic usage of the library:

import fabio, pyFAIdata = fabio.open("Pilatus1M.edf").dataai =
pyFAI.load("Pilatus1M.poni")tth, I = ai.integrate1d(data, 1000,
unit="2th_deg",\ method="numpy")Output variables’ space (_r_ , _q_ or
$2\theta$) and units can be chosen with the _unit_ keyword. The _method_
keyword is used to choose one of the available algorithms for the integration.
These algorithms will be described in this contribution. However, the
experiment reported here will be limited to 1D full azimuthal integration,
with a planar detector, orthogonal to the incoming beam. In this case the
conics described by the beam on the detector are concentric circles. The
generic geometry used in pyFAI has already been described in [pyFAI_ocl].

###  3.1 Test case

To let the reader have an idea of the scale of the problem and the
performances needed, we will work on the simulated image of gold powder
diffracting an X-ray beam of wavelength = 1.0e-10m (the intensity of all rings
is the same). The detector, which has a pixel size of 1e-4m (2048x2048
pixels), is placed at a distance of 0.1 m from the sample, orthogonal to the
incident beam, which coincides with the centre of the rings. Figure 2
represents the input diffraction image (upper part) and the integrated profile
along the azimuthal angle (lower part). The radial unit in this case is simply
the radius calculated from $r=\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$, while
crystallographers would have preferred $2\theta$ or the scattering vector
length _q_.

Figure 2: Simulated powder diffraction image (top) and integrated profile
(bottom).

###  3.2 Naive implementation

The initial step of any implementation is calculating the radius array, from
the previous formula. Using Numpy’s slicing feature one can extract all pixels
which are between r1 and r2 and average out their values:

def azimint_naive(data, npt, radius): rmax = radius.max() res =
numpy.zeros(npt) for i in range(npt): r1 = rmax * i / npt r2 = rmax * (i+1) /
npt mask_r12 = numpy.logical_and((r1 <= radius), (radius < r2)) values_r12 =
data[mask_r12] res[i] = values_r12.mean() return resThe slicing operation
takes tens of milliseconds and needs to be repeated thousands of times for a
single image, making each integration last 40 seconds, something that is
unacceptably slow.

###  3.3 Numpy histograms

The naive formulation in section 3.2 can be re-written using histograms. The
_mean_ call can be replaced with the ratio of the sum of all values divided by
the number of contributing pixels:

values_r12.mean() = values_r12.sum() / mask_r12.sum()The denominator,
_mask_r12.sum()_ , can be obtained from the histogram of _radius_ values and
the numerator, _values_r12.sum()_ from the histogram of _radius_ weighted by
the _data_ in the image:

def azimint_hist(data, npt, radius): histu = np.histogram(radius, npt)[0]
histw = np.histogram(radius, npt, weights=data)[0] return histw / histuThis
implementation takes about 800ms which is much faster than the loop written in
Python, but can be optimised by reading the radius array from central memory
only once.

###  3.4 Cython implementation

Histograms were re-implemented using [Cython] to generate simultaneously both
the weighted and the unweighted histograms with a single memory read of the
radius array. The better use of the CPU cache decreases the integration time
significantly (down to 150ms on a single core).

####  3.4.1 OpenMP support in Cython

To accelerate further the code we decided to parallelise the [Cython] code
using [OpenMP]. While the implementation was fast, the results we obtained
were wrong (by a few percent) due to write conflicts, not protected by
atomic_add operations. Apparently the use of atomic operation is still not yet
possible in [Cython] (summer 2014). Multi-threaded histogramming was made
possible by having several threads running simultaneously, each working on a
separate histogram, which implies the allocation of much more memory for
output arrays.

Implement. |  Exec. time (ms)
---|---
loop + mean |  44000
np.histogram |  829
Cython 1_th |  149
Cython 2_th |  81
Cython 4_th |  59
Cython 8_th |  41
Cython 16_th |  48
TABLE II: Azimuthal integration time for a 4 Mpix image measured on two Xeon
E5520 (2x 4-core hyper-threaded at 2.2 GHz)

The gains in performance obtained by this method (see table I) were minor,
especially when using more than 2 threads, illustrating the limits of the
paralellisation scheme. The only way to go faster is to start thinking in
parallel from the beginning and re-design the algorithm so that it works
natively with lots of threads. This approach is the one taken by [OpenCL],
where thousands of threads are virtually running in parallel, and is described
in paragraph 5.

##  4 Pixel splitting

Pixel splitting is what occurs when a pixel of the detector spans over more
than one of the bins of the histogram. When this happens, the contribution to
each of the bins involved is assumed to be proportional to the area of the
pixel segment that falls into that bin. The goal behind the addition of extra
complexity to the code is that the results obtained this way ought to be less
noisy than the case where pixel splitting is ignored. This becomes more
apparent when the number of pixels falling into each bin is small like for
example for 2D integration. Figure 3 presents the results of such an
integration, performed using histograms on the top image, i.e. without pixel
splitting. Some high frequency patterns are visible near the beam center on
the left-hand side of this figure. The bottom image was produced using pixel
splitting and is unharmed by such defects, which are related to low
statistics. Note that for 2D integration, this transformation looks like an
interpolation, but interpolation neither guarantees the conservation of the
signal $\sum{image}=\sum{weighted\ histogram}$ nor that of the pixels
$\sum{unweighted\ histogram}=number\ of\ pixels$.

Figure 3: Bi-dimensional azimuthal integration of the gold diffraction image
using (bottom) or not (top) pixel splitting

###  4.1 Bounding Box

The first way pixel splitting was implemented was with a bounding box like in
Fit2D [FIT2D]. In this case we are working with an abstraction of the pixel.
This is represented by a rectangular box circumscribing the actual pixel, with
two sides parallel to the radial axis and the other two of unit length.
Presently, instead of calculating the contribution of each segment of the
pixel based on its area, we use the area of the bounding box segment instead.
This greatly simplifies the algorithm’s flow, providing good performance.

The algorithm loops over all the pixels of the detector, adding their
contributions to the appropriate bins. When the whole pixel falls into only
one bin, there is no pixel splitting and the algorithm proceeds as in the case
of the simple histogram. If the pixel spans over more than one bin, the
contribution to the two outermost bins (left and right) is calculated first
and added to them. Then, the remaining contribution is evenly distributed
among the “internal” bins (if any). Finally, the ratio of the two histograms
is calculated and returned.

The trade-off of using this simplistic pixel splitting is an overestimation of
the pixel size, hence a slight blurring of the signal.

###  4.2 Full Pixel Splitting

In an effort to farther improve the quality of the results of the azimuthal
integration, another pixel-splitting scheme was devised, in which no
abstraction takes place and the pixel splitting works using the area of the
actual pixel segments (assuming they are straight lines). This introduces some
additional complexity to the calculations, making the process a bit slower.

As before, the algorithm first has to check if pixel splitting occurs. In the
case it does not, the pixel is processed like in the case of the simple
histogram. Otherwise the pixel is split according to the following steps.
Firstly, a function for each of the lines that make up the sides of the pixel
being processed is defined by calculating the slope and the point of
intersection. The area of the pixel is also required. Next, the algorithm
loops over the bins that the pixel spans over and proceeds to integrate the
four functions that were previously defined over the bounds of each bin.
Taking the absolute value of the sum of all these integrals will yield the
area of the pixel segment. Now, the contributions to the histograms are
calculated using these areas. The difficult part here was the definition of
the limits of each of the integrals in a way that would not hinder the
performance by adding many conditionals.

###  4.3 Discussion on the statistics

Using either of the two pixel splitting algorithms results in some side
effects that the user should be aware of: The fact that pixels contributing to
neighbouring bins in the histogram creates some cross-correlation between
those bins, affecting, this way, the statistics of the results in a
potentially unwanted manner [Stat].

##  5 More parallelisation

For faster execution, one solution is to use many-core systems, such as
Graphical Processing Units (GPUs) or accelerators, like the Xeon-Phi from
Intel. Those chips allocate more silicon for computing (arithmetic logic units
- ALUs) and less to branch prediction, memory pre-fetching and cache
coherency, in comparison with CPUs. Our duties as programmers is to write the
code that maximises the usage of ALUs without relying on pre-fetcher and other
commodities offered by normal processors.

Typical GPUs have tens (to hundreds) of compute units able to schedule and run
dozens of threads simultaneously (in a Single Instruction Multiple Data way).
OpenCL allows the execution of the same code on processors, graphics cards or
accelerators (see table III) but the memory access pattern is important in
order to make the best use of them. Finally, OpenCL uses just-in-time (JIT)
compilation, which looks very much like Python interpreted code when
interfaced with [PyOpenCL] (thanks to the compilation speed and the caching of
the generated binary).

Vendor / driver |  Intel |  AMD |  AMD |  Nvidia |  Nvidia |  Intel
---|---|---|---|---|---|---
Model |  2xE5-2667 |  2xE5-2667 |  V7800 |  Tesla K20 |  GeForce 750Ti |  Phi 5110
Type |  CPU |  CPU |  GPU |  GPU |  GPU |  ACC
Compute Unit |  12 |  12 |  18 |  13 |  5 |  4x69
Compute Element/CU |  4:AVX |  1 |  80 |  4x8:Warp |  4x8:Warp |  16:AVX512
Core frequency |  2900 MHz |  2900 MHz |  700 MHz |  705 MHz |  1100 MHz |  1052
Mem. Bandwidth |  102 GB/s |  102 GB/s |  128 GB/s |  208 GB/s |  88 GB/s |  320 GB/s
TABLE IV: Few OpenCL devices we have tested our code on.

###  5.1 Parallel algorithms

Parallelisation of complete algorithms consists, most of the time, in their
decomposition into parallel blocks. There are a few identified parallel
building blocks like:

  * •

Map: apply the same function on all elements of a vector

  * •

Scatter: write multiple outputs from a single input, needs atomic operation
support

  * •

Gather: write a single output from multiple inputs

  * •

Reduction: single result from a large vector input, like an inner product

  * •

Scan: apply subsequently an operation to all preceding elements on an vector
like np.cumsum

  * •

Sort: There are optimised sorter for parallel implementation.

These parallel building blocks will typically be one individual kernel or a
few, since kernel execution synchronises the global memory in OpenCL. Parallel
algorithmics deal with how to assemble those blocks to implement the required
features.

###  5.2 Parallel azimuthal integration

Azimuthal integration, like histogramming, is a scatter operation, and hence
requires the support of atomic operations (in our case with double precision
floats). As Cython does not (yet) support atomic operations, enabling OpenMP
parallelisation results in a module that, while being functional, gives the
wrong results (we measured 2% errors on 8 cores)

To overcome this limitation, instead of looking at where input pixels go to in
the output curve, we focus on where the output bin comes from in the input
image. This transformation is called a “scatter to gather” transformation and
requires atomic operations. In our case, it was implemented as a single
threaded [Cython] module.

The correspondence between pixels and output bins can be stored in a look-up
table (LUT) together with the pixel weight (ratio of areas) making the
integration look like a simple (if large and sparse) matrix vector product.
The LUT size depends on whether pixels are split over multiple bins and in
order to exploit the sparse structure, both the index and the weight of each
pixel have to be stored.

By making this change we switched from a “linear read / random write” forward
algorithm to a “random read / linear write” backward algorithm which is more
suitable for parallelisation. For optimal memory access patterns, the array of
the LUT may be transposed depending on the underlying hardware (CPU vs GPU).

####  5.2.1 Optimisation of the sparse matrix multiplication

The compressed sparse row (CSR) sparse matrix format was introduced to reduce
the size of the data stored in the LUT. This algorithm was implemented both in
[Cython]-[OpenMP] and [OpenCL]. Our CSR representation contains _data_ ,
_indices_ and _indptr_ (row index pointer) so it is fully compatible with the
_scipy.sparse.csr.csr_matrix_ constructor from [SciPy]. This representation is
a _struct of array_ which is better suited to GPUs (strided memory access)
while LUT is an _array of struct_ , known to be better adapted to CPU (better
use of cache and pre-fetching). The CSR approach presents a double benefit:
first, it reduces the size of the storage needed, as compared to the LUT, by a
factor two to three, and gives the opportunity of working with larger images
on the same hardware. Secondly, the CSR implementation in [OpenCL] is using an
algorithm based on multiple parallel reductions where all threads within a
workgroup are collaborating to calculate the content of a single bin. This
makes it very well suited to run on many-core systems where hundreds to
thousands of simultaneous threads are available.

####  5.2.2 About numerical precision

Knowing the tight energy constraints, the future of high performance computing
depends on the capability of programs to use the suitable precision for their
calculations. As our detectors provide a sensitivity of 12 to 20 bits/pixel,
performing all calculations in double precision (with 52 bits mantissa) might
seem excessive, the 24 bits mantissa of single precision float being a better
choice for the task (with no precision drop). Moreover, GPU devices provide
much more computing power in single precision than in double. This factor
varies from 2 on high-end professional GPUs like Nvidia Tesla to 24 on most
consumer grade devices.

When using [OpenCL] for GPUs we used compensated arithmetics (or [Kahan]
summation), to reduce the error accumulation in the histogram summation (at
the cost of more operations). This allows numerically accurate results to be
obtained even on cheap consumer grade hardware with the use of single
precision floating point arithmetic (32 bits). Double precision operations are
currently limited to high-price / high-performance GPUs, optimised exactly for
that purpose. The additional cost of Kahan summation (4x more arithmetic
operations) is hidden by smaller data types, a higher number of single
precision units and the fact that GPUs are usually limited by the memory
bandwidth anyway.

The performances of the parallel azimuthal integration can reach 750 MPixel/s
on recent computers with a mid-range graphics card. On multi-socket servers
featuring high-end GPUs like Tesla cards, the performances are equivalent, but
with the added benefit of working with multiple detectors simultaneously.

##  6 Benchmarks

We present the results from several benchmark tests done using the different
algorithm options available in PyFAI. All benchmarks were performed using the
same bounding box pixel splitting scheme and the resulting integrated profiles
are of equivalent quality. Execution speed has been measured using the
_timeit_ module, averaged over 10 iterations (best of 3). The processing is
performed on 1, 2, 4, 6, 12 and 16 Mpixel images, with pixel ranges of either
16 or 32 bits (int or uint), taken from actual diffraction experiments, which
are part of the pyFAI test suite.

One small note on the benchmarks that follow. The casting for the 12 Mpixel
image was done by one thread on the CPU. That is why the processing time of
the 16 Mpixel image appears to be shorter than that of the 12 Mpixel one.

###  6.1 Choice of the algorithm

The LUT contains pairs of an index and a coefficient, hence it is an _array of
struct_ pattern which is known to make best use of CPU caches. On the
contrary, the CSR sparse matrix representation is a _struct of array_ which is
better adapted to GPU. As we can see in figure 4, both LUT and CSR outperform
the serial code, and both behave similarly: the penalty of the _array of
struct_ in CSR is counter-balanced by the smaller chunk of data to be
transferred from central memory to CPU.

Figure 4: Comparison of azimuthal integration speed obtained using serial
implementation versus parallel implementations with LUT and CSR sparse matrix
representation on two Intel Xeon E2667.

###  6.2 OpenMP vs OpenCL

The gain in portability obtained by the use of OpenCL does not mean a
sacrifice in performance when the code is run on a CPU, as we can see in
figure 5: the OpenCL implementation outperforms the OpenMP one, in all the
different CPUs it was tested on. This could be linked to the better use of
SIMD vector units by OpenCL. The dual Xeon E5520 (a computer from 2009),
running at only 2.2 GHz shows pretty good performances compared to more recent
computers when using OpenMP: it was the only one with activated hyper-
threading.

Figure 5: Comparison of the azimuthal integration speed between the OpenMP and
OpenCL implementations.

The choice of the OpenCL driver on CPU affects the performance of PyFAI
(figure 6): on the Intel Xeon E5-1607 (Ivy bridge core), the Intel driver
clearly outperforms the AMD driver. This can be attributed to new SIMD
instructions (AVX), supported by the Intel driver but not by the AMD one. On
the older Intel Xeon E-5520 (Nehalem core) which lacks those extensions, the
difference in speed is much less.

Figure 6: The effects of OpenCL driver selection on performance on different
generations of CPUs.

###  6.3 GPUs and Xeon Phi

Figure 7 compares the integration speed of the LUT and CSR implementation on
two GPUs. The CSR implementation, thanks to the multiple collaborative
parallel reductions, runs much faster on all the devices used, compared to the
LUT one. Another benefit of the CSR implementation when it comes to GPUs is
its lower memory usage. The ATI GPU used in this study features only 1 GB of
memory usable by OpenCL, limiting the processable size of the system. This is
the reason the benchmarks stop before reaching the largest image size. 4
Mpixel images are the largest images processable with the LUT implementation,
but 12 Mpixel images are processable using the CSR one.

Figure 7: Comparison of the azimuthal integration speed between the LUT and
CSR implementations on GPUs.

In figure 8, we have gathered the results from all of the many-core devices
available to us, including several GPUs as well as Intel’s Xeon Phi. As one
can see, Xeon Phi (from 2012) matches the performance of the AMD GPU from
2010. What is surprising though, is how well the consumer grade Nvidia GeForce
750Ti performs in comparison to high-end _Kepler_ cards (Titan, Tesla K20)
costing only a fraction of their price.

Figure 8: Comparison of the performances for several many-core accelerators:
GPUs and Xeon Phi.

###  6.4 Kernel timings

As stated previously, the benchmark tests were performed using the _timeit_
module from Python on the last line of the code snippet described in section
3. One may wonder what is the actual time spent in which part of the OpenCL
code and how much is the Python overhead. This analysis has been done using
the profiling tools of OpenCL which measured the execution of every action put
in queue. To be able to perform the azimuthal integration, the image is first
transfered to the device (GPU), then casted from integer to float. All pixel-
wise correction (dark current subtraction, flat field normalization, solid-
angle and polarization factor correction) are applied in a single pass over
each pixel of the image. Output arrays are initialised to zero, by a separate
kernel (memset) before the actual sparse-matrix-dense-vector multiplication.
Finally the three output buffers are retrieved from the device.

Table V shows the execution time measured on the GeForce Titan (controlled by
a pair of Xeon 5520). The first entry in the table is the total execution time
at the Python level, as measured by _timeit_ : 2 ms, while the second is the
sum of all of the execution times measured by the OpenCL profiler: 1.4 ms,
which highlights how little the Python overhead can be (<40%). The most time-
consuming part of the whole process is by far the memory transfer of the image
(H->D meaning Host to Device, 0.8ms). All vendors are currently working on an
unified memory space, which will be available for OpenCL 2.0, which will
reduce the time spent in transfers and simplify programming. Finally the
azimuthal integration takes up only 0.4 ms, that is, 20% of the total run
time. If one focuses only on the timing of the integration kernel, then he
would wrongly conclude that pyFAI is able to match the speed of the fastest
detectors. For example, the 2 ms of processing time for a 1 Mpixel image of 32
bit integers, correspond to a processing rate of 2 GB/s, while our fastest
storage solutions (solid-state drives) are currently only able to provide half
of that.

Python total |  2.030ms
---|---
OpenCL |  1.445ms
H->D image |  0.762ms
cast |  0.108ms
memset |  0.009ms
correction |  0.170ms
integrate |  0.384ms
D->H ratio |  0.004ms
D->H u_hist |  0.004ms
D->H w_hist |  0.004ms
TABLE VI: OpenCL profiling of the integration of a Pilatus 1M image (981x1043
pixels of signed 32 bits integers) on a GeForce Titan, running on a dual Xeon
5520.

###  6.5 Configuration and Drivers used

The computer hosting the two Intel Xeon E5-2667 (2x6 cores each, 2.9 GHZ,
without hyper-threading, 8x8 GB of RAM) is a Dell PowerEdge R720 with both a
Tesla K20 and an Intel Xeon phi accelerator, running Debian 7. The computer
hosting the two Intel Xeon E5520 (2x4cores, 2.27 GHz, hyper-threaded, 6x2 GB
of RAM) is a Dell T7500 workstation with two Nvidia GPUs: GeForce 750Ti and
Titan, running Debian 7. The computer hosting the Intel Xeon E5-1607
(1x4cores, 3.0 GHz, without hyper-threading, 2x4 GB of RAM) is a Dell T3610
workstation with two GPUs: Nvidia GeForce 750Ti and AMD FirePro V7800, running
Debian 8/Jessie.

In addition to the Debian operating system, specific OpenCL drivers were
installed:

  * •

Intel OpenCL drivers V4.4.0-117 + MPSS stack v3.2.3

  * •

AMD APP drivers 14.4

  * •

Nvidia CUDA drivers 340.24-2

##  7 Project description

PyFAI is open-source software released under the GPL license available on
GitHub (https://github.com/kif/pyFAI). PyFAI depends on Python v2.6 or v2.7
and [NumPy]. In order to be able to read images from various X-ray detectors,
pyFAI relies on the [FabIO] library. Optional [OpenCL] acceleration is
provided by [PyOpenCL]. Graphical applications for calibration and integration
rely on [matplotlib], [PyQt] and SciPy [SciPy] for image processing. A C
compiler is needed to build the [Cython] code from the related sources. PyFAI
is packaged and available in common Linux distributions like Debian and Ubuntu
but it is also tested and functional under Windows and MacOSX.

##  8 Conclusions

This contribution shows how one of the most central algorithm in
crystallography has been implemented in Python, optimised in Cython and ported
to many-core architectures using PyOpenCL. A 15x speed-up factor has been
obtained by switching from binary code to the OpenCL code running on GPUs
(400x vs NumPy). Some of the best performances were obtained on a mid-range
consumer grade Nvidia GeForce 750Ti thanks to the new _Maxwell_ generation
chip running as fast as high-end graphics based on the _Kepler_ architecture
(like the Titan), and literally outperforming both AMD GPUs and the Xeon-Phi
accelerator card. Programming CPUs in parallel is as easy as programming GPUs
via the use of PyOpenCL interfaced with Python.

##  9 Acknowledgements

Claudio Ferrero (head of the Data Analysis Unit) and Andy Götz (head of the
Software Group) are acknowledged for supporting the development of pyFAI. The
porting of pyFAI to OpenCL would have not been possible without the financial
support of LinkSCEEM-2 (RI-261600), granting the contracts of Dimitris
Karkoulis who started the GPU porting, Zubair Nawaz who ported image
distortion and one of the authors (G. Ashiotis) who is working on CSR, pixel
splitting and other algorithms. Finally, the authors would like to acknowledge
their colleagues involved in the development of the library, especially Aurore
Deschildre and Frédéric Picca. The authors would like to thank all X-ray beam-
lines promoting pyFAI and providing resources to further develop it: ESRF
BM01, ID02, ID11, ID13, ID15, ID16, ID21, ID23, BM26, ID29, BM29 and ID30; and
also in other institutes like Soleil, Petra3, CEA, APS who provide feedback,
bug reports and patches to the library.

## References

  * [Cython] S. Behnel, R. Bradshaw, C. Citro, L. Dalcin, D.S. Seljebotn and K. Smith. _Cython: The Best of Both Worlds_ Comput. Sci. Eng., 13(2):31-39, 2011.
  * [FabIO] E. B. Knudsen, H. O. Sorensen, J. P. Wright, G. Goret and J. Kieffer. _FabIO: easy access to two-dimensional X-ray detector images in Python_ , J. Appl. Cryst., 46:537-539, 2013.
  * [FIT2D] A. Hammersley, O. Svensson, M. Hanfland, A. Fitch and D. Hausermann. _Two-dimensional detector software_ , High Press. Res., 14:235–248, 1996.
  * [H5Py] A. Collette. _Python and HDF5_ ISBN 978-1-4493-6783-1, (2013)
  * [Kahan] W. Kahan. _Pracniques: Further Remarks on Reducing Truncation Errors_ , Commun. ACM,8(1):40-, Jan. 1965
  * [matplotlib] J. D. Hunter. _Matplotlib: A 2D Graphics Environment_ , Comput. Sci. Eng., 9(3):90-95, 2007.
  * [NumPy] T. E. Oliphant. _Python for Scientific Computing_ , Comput. Sci. Eng., 9(3):10-20, 2007.
  * [OpenCL] J.E. Stone, D. Gohara and G. Shi. _OpenCL: A Parallel Programming Standard for Heterogeneous Computing Systems_ , Comput. Sci. Eng., 12(3):66-73, 2010.
  * [OpenMP] OpenMP Architecture Review Board. _OpenMP Application Program Interface Version 3.0_ , 2008.
  * [pyFAI] J. Kieffer and D. Karkoulis. _PyFAI, a versatile library for azimuthal regrouping_ , Journal of Physics: Conference Series, 425:202012, 2013.
  * [pyFAI_ocl] J. Kieffer and J.P. Wright. _PyFAI: a Python library for high performance azimuthal integration on GPU_ , Powder Diffraction, 28S2:1945-7413, 2013.
  * [PyOpenCL] A. Klöckner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov and A. Fasih. _PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation_ Parallel Computing, 38(3):157-174, 2012.
  * [PyQt] Mark Summerfield. _Rapid GUI Programming with Python and Qt: The Definitive Guide to PyQt_ , ISBN 0132354187 (2007).
  * [SciPy] E. Jones, T. E. Oliphant and P. Peterson, _SciPy: Open source scientific tools for Python_ , 2001.
  * [SPD] P. Bösecke. _Reduction of two-dimensional small- and wide-angle X-ray scattering data_ , J. Appl. Cryst., 40:s423–s427, 2007.
  * [Stat] X. Yang, P. Juhás and S. J. L. Billinge. _On the estimation of statistical uncertainties on powder diffraction and small-angle scattering data from two-dimensional X-ray detectors_ , J. Appl. Cryst., 47:1273-1283, 2014.