Also at ]School of Engineering and Information Technology, University of New
South Wales, Canberra, Australia

Also at ]Centre for Audio, Acoustics and Vibration, University of Technology
Sydney, Sydney, Australia

Also at ]Centre for Audio, Acoustics and Vibration, University of Technology
Sydney, Sydney, Australia

# Acoustic radiation force and radiation torque beyond particles:
Effects of non-spherical shape and Willis coupling

Shahrokh Sepehrirahnama<EMAIL_ADDRESS>Centre for Audio,
Acoustics and Vibration, University of Technology Sydney, Sydney, Australia
Sebastian Oberst  [  Centre for Audio, Acoustics and Vibration, University of
Technology Sydney, Sydney, Australia  Yan Kei Chiang  [  School of Engineering
and Information Technology, University of New South Wales, Canberra, Australia
David Powell  [  School of Engineering and Information Technology, University
of New South Wales, Canberra, Australia

###### Abstract

Acoustophoresis deals with the manipulation of sub-wavelength scatterers in an
incident acoustic field. The geometric details of manipulated particles are
often neglected by replacing them with equivalent symmetric geometries such as
spheres, spheroids, cylinders or disks. It has been demonstrated that
geometric asymmetry, represented by Willis coupling terms, can strongly affect
the scattering of a small object, hence neglecting these terms may miss
important force contributions. In this work, we present a generalized
formalism of acoustic radiation force and radiation torque based on the
polarizability tensor, where Willis coupling terms are included to account for
geometric asymmetry. Following Gorkov’s approach, the effects of geometric
asymmetry are explicitly formulated as additional terms in the radiation force
and torque expressions. By breaking the symmetry of a sphere along one axis
using intrusion and protrusion, we characterize the changes in the force and
torque in terms of partial components, associated with the direct and Willis
Coupling coefficients of the polarizability tensor. We investigate in detail
the cases of standing and travelling plane waves, showing how the equilibrium
positions and angles are shifted by these additional terms. We show that while
the contributions of asymmetry to the force are often negligible for small
particles, these terms greatly affect the radiation torque. Our presented
theory, providing a way of calculating radiation force and torque directly
from polarizability coefficients, shows that in general it is essential to
account for shape of objects undergoing acoustophoretic manipulation, and this
may have important implications for applications such as the manipulation of
biological cells.

Acoustic bianisotropy, Primary Force, Primary Torque, Non-spherical Scattering

††preprint: APS/123-QED

##  I Introduction

Acoustic radiation forces play a key role in the field of acoustic particle
manipulation, also known as acoustophoresis [1, 2, 3, 4, 5, 6, 7, 8]. Acoustic
sorting, separation, levitation and other similar applications have been
developed to manipulate particulate phase in a fluid, e.g. to migrate
biological cells to certain locations using incident plane waves, by inducing
a radiation force field [3, 7, 9, 10, 11, 12, 13, 14, 8, 15, 16], as
schematically shown in Fig. 1. In the design of such applications, it is
customary to treat the objects in the host fluid as spheres or other simple
geometries and neglect the details of their shapes. In most cases, this leads
to a design based solely on acoustic radiation force and neglecting the
radiation torque. When a sub-wavelength object is treated as a particle, its
dynamic equilibrium and force balance is independent of its shapes and its
rotation is neglected. This assumption can be reasonable for small objects
with approximately spherical shape; however, non-spherical objects lacking
radial symmetry may be better approximated as rigid bodies than point-like
particles.

Figure 1: Principle of acoustophoresis with acoustic radiation force and
torque exerted on an object with arbitrary shape, for any choice of incident
pressure wave and the angle of incidence. The components of the acoustic
radiation force and torque and the notation for their normalized contrast
factors, $\mathbf{Q}$ and $\mathbf{Z}$, respectively, are illustrated.

Acoustic radiation stress consists of radiation pressure, which is related to
radiated momentum through scattering fields, and Reynolds stress arising from
surface oscillation. The primary radiation force acting on a scatterer
corresponds to the incident-scattering portion of the time-averaged radiation
stresses [17, 18, 19, 2, 20, 21, 22, 23, 24, 25]. In the case of multiple
objects, pairwise secondary radiation forces, also known as acoustic
interaction forces, emerge from the scattering-scattering part of the
radiation stresses [26, 21, 22, 23, 24, 25]. The acoustic radiation torque is
the resultant moment due to such stresses with respect to the centroid of the
object [27, 28]. Acoustic radiation force and torque contribute to the dynamic
equilibrium state of the scatterers in combination with other forces such as
hydrodynamic drag or gravitational force. The magnitudes of acoustic radiation
force and torque are proportional to the incident energy density and their
direction is determined from their normalized value, also referred to as the
acoustic contrast factor [18, 2, 24], as shown in Fig. 1.

Current theories of acoustic radiation force and torque use either the
radiation stresses at the surface of the object [17, 18, 29, 30, 31], or
convert the surface integral of radiation stresses in the far-field to a
volume integral using the Divergence theorem [19, 2, 24]. The advantage of the
latter approach is that the object can be replaced by a set of acoustic
multipole sources to evaluate the volume integral, capturing the essential
geometric features and material properties of the object with the simplest
possible model. Acoustic radiation force and torque have been developed
analytically or modelled numerically for objects with spheroidal, cylindrical
and disk shapes [31, 27, 32, 33, 34, 35, 36, 37, 38]. However, these shapes
exhibit a high degree of symmetry (at least axial symmetry), and are not fully
representative of three-dimensional objects with arbitrary shape. Although the
numerical approaches based on surface or volume integrals can be applied to
objects of arbitrary shape, they do not directly show how the geometric
asymmetry impacts on the force and torque. For acoustic scattering from small
particles, it has been shown that asymmetry can be accounted for by
incorporating Willis coupling terms into the multipole tensor [39, 40, 41, 42,
43]. However, to date there has been no investigation of the role that these
terms may play in the acoustic radiation force and torque.

In this work, we derive a rigorous mathematical formulation of acoustic
radiation force and acoustic radiation torque that includes the shape
complexity as represented by the Willis coupling terms. We make use of the
far-field approach, also referred to as the Gorkov approach, and show how the
particle responds to both the incident pressure and velocity fields. We
present a general formalism, but consider in detail the cases of plane
travelling and standing waves. Starting from a symmetric shape, we show how
introducing geometric protrusions or instructions controls the Willis
coupling, and hence the additional radiation force and torque terms.

##  II Theory

###  II.1  Scattering of sub-wavelength objects and polarizability

The acoustic wave propagation in a lossless fluid is governed by the wave
equation , which is expressed in terms of acoustic pressure $p$ as follows,

$\displaystyle\nabla^{2}p=\frac{1}{c_{f}^{2}}\partial_{tt}p,$ (1)

where $\nabla^{2}$ represents the Laplacian operator, $c_{f}$ denotes the
speed of sound in the fluid medium and $\partial_{t}=\frac{\partial}{\partial
t}$. The acoustic density $\rho$ and velocity $\mathbf{v}$ are related to the
pressure as follows,

$\displaystyle
p=c_{f}^{2}\rho,\quad\partial_{t}\rho=-\rho_{f}\nabla\cdot\mathbf{v},\quad
p=-\rho_{f}\partial_{t}\phi,\quad\mathbf{v}=\nabla\phi,$ (2)

where $\kappa_{f}=\rho_{f}c_{f}^{2}$ is the fluid compressibility, $\rho_{f}$
is the mean fluid density, and $\phi$ denotes the velocity potential. Acoustic
pressure, density and velocity fields are time-harmonic,

$\displaystyle p=p(\mathbf{x})e^{-j\omega
t},\quad\mathbf{v}=\mathbf{v}(\mathbf{x})e^{-j\omega
t},\quad\rho=\rho(\mathbf{x})e^{-j\omega t},$ (3)

where $t$, $\omega$ and $\mathbf{x}$ denote the time, the angular frequency
and the position vector, respectively.

In the Rayleigh limit, the monopole-dipole approximation of the scattering
field of a sphere is given as [19, 2, 24]

$\displaystyle\begin{split}\phi_{s}&\approx-\frac{a_{s}^{3}}{3\rho_{f}}f_{1}\partial_{t}\rho_{i}\frac{e^{jkr}}{r}-\frac{a_{s}^{3}}{2}\nabla\cdot\Big{(}f_{2}\mathbf{v}_{i}\frac{e^{jkr}}{r}\Big{)},\\\
f_{1}&=1-\frac{\kappa_{s}}{\kappa_{f}},\qquad
f_{2}=\frac{2\rho_{s}-2\rho_{f}}{2\rho_{s}+\rho_{f}},\end{split}$ (4)

where $\rho_{i}$ and $\mathbf{v}_{i}$ denote the value of incident density and
velocity fields, respectively, $a_{s}$ denotes the sphere radius, $r$ is the
radial distance measured from the center, $\rho_{s}$ and $\kappa_{s}$ denote
the density and compressibility of the scatterer, respectively, and
$k=\omega/c_{f}$ denotes the wave number. To generalize this approach to
arbitrary small objects, including those exhibiting Willis coupling, we employ
the method of multipole moments, up to dipole accuracy, and express the
scattered pressure as [44, 40]

$\displaystyle\begin{split}p_{s}&\approx-\omega^{2}\Omega
MG+\omega^{2}\pmb{\nabla}\cdot\Big{[}\Omega\mathbf{D}G\Big{]},\quad
G=G\big{(}kr\big{)}=\frac{e^{jkr}}{4\pi r},\\\
\phi_{s}&=\frac{p_{s}}{j\rho_{f}\omega},\end{split}$ (5)

where $G$ denote the Green’s function in 3D domain, and $M$ and $\mathbf{D}$
denote the volumetric monopole and dipole moments, respectively. The term
$a_{s}^{3}$ from Eq. (4) is generalized to the volume of the object, denoted
$\Omega$. Furthermore, a characteristic size, denoted by $a$, is required for
non-spherical shapes to calculate the Rayleigh index of $ka$. For this study,
we propose

$\displaystyle a$ $\displaystyle=3\frac{\Omega}{\Gamma},$ (6)

where $\Gamma$ is the outer surface area. This allows us to incorporate the
volume and surface area as two measures of three-dimensional geometries in our
numerical analysis. A factor of $3$ is included to normalize $a$ to the radius
$a_{s}$ for the case of a spherical object.

The scattering moments are expressed in terms of incident pressure $p_{i}$ and
velocity $\mathbf{v}_{i}$ fields through the polarization tensor
$\pmb{\alpha}$ as follows,

$\displaystyle\begin{split}\begin{bmatrix}M\\\
\mathbf{D}\end{bmatrix}=\frac{1}{\Omega}\begin{pmatrix}\alpha_{pp}&\pmb{\alpha}_{pv}^{T}\\\
\pmb{\alpha}_{vp}&\pmb{\alpha}_{vv}\end{pmatrix}\begin{bmatrix}p_{i}\\\
\mathbf{v}_{i}\end{bmatrix},\qquad\pmb{\alpha}_{pv}=\begin{bmatrix}\alpha_{pv}^{x}\\\
\alpha_{pv}^{y}\\\ \alpha_{pv}^{z}\end{bmatrix},\\\
\pmb{\alpha}_{vp}=\begin{bmatrix}\alpha_{vp}^{x}\\\ \alpha_{vp}^{y}\\\
\alpha_{vp}^{z}\end{bmatrix},\qquad\pmb{\alpha}_{vv}=\begin{pmatrix}\alpha_{vv}^{xx}&\alpha_{vv}^{xy}&\alpha_{vv}^{xz}\\\
\alpha_{vv}^{yx}&\alpha_{vv}^{yy}&\alpha_{vv}^{yz}\\\
\alpha_{vv}^{zx}&\alpha_{vv}^{zy}&\alpha_{vv}^{zz}\end{pmatrix}\end{split}$
(7)

where $\pmb{\alpha}_{\rho v}$ and $\pmb{\alpha}_{v\rho}$ denote the Willis
coupling coefficients, $\pmb{\alpha}_{vv}$ denotes direct-dipole polarization
tensor, and superscript $T$ denotes the transpose operator. The Cartesian form
of these sub-tensors are given in Eq. (7). The entries of $\pmb{\alpha}$ are
often normalized as follows [40],

$\displaystyle\begin{bmatrix}-\sqrt{3}M\\\ jk\mathbf{D}\end{bmatrix}$
$\displaystyle=\frac{\overline{\pmb{\alpha}}}{\Omega}\begin{bmatrix}\frac{1}{\sqrt{3}}p_{i}\\\
\rho_{f}c_{f}\mathbf{v}_{i}\end{bmatrix},\quad\overline{\pmb{\alpha}}=\begin{pmatrix}-3\alpha_{pp}&\frac{-\sqrt{3}}{\rho_{f}c_{f}}\pmb{\alpha}_{pv}^{T}\\\
jk\sqrt{3}\pmb{\alpha}_{vp}&\frac{jk}{\rho_{f}c_{f}}\pmb{\alpha}_{vv}\end{pmatrix}.$
(8)

By employing the reciprocity principle for Green’s function, it has been
proven that [40]

$\displaystyle\overline{\pmb{\alpha}}=\overline{\pmb{\alpha}}^{T-}=\begin{pmatrix}\overline{\alpha}_{pp}&-\overline{\pmb{\alpha}}_{vp}^{T}\\\
-\overline{\pmb{\alpha}}_{pv}&\overline{\pmb{\alpha}}_{vv}\end{pmatrix}.$ (9)

This yields the relation between the Willis coupling coefficients
$\overline{\pmb{\alpha}}_{vp}=-\overline{\pmb{\alpha}}_{pv}$, which are later
required for the formulation of the acoustic radiation force and the acoustic
radiation torque. By substituting Eq. (7) into Eq. (5) and comparing the
results with Eq. (4), the monopole and dipole moments for a spherical
scatterer become [2, 24]

$\displaystyle\begin{bmatrix}M\\\
\mathbf{D}\end{bmatrix}_{sphere}=\frac{1}{\Omega}\begin{pmatrix}\Omega\rho_{f}\kappa_{f}f_{1}&\mathbf{0}^{T}\\\
\mathbf{0}&-\frac{3j\Omega}{2\omega}\rho_{f}f_{2}\mathbf{I}\end{pmatrix}\begin{bmatrix}p_{i}\\\
\mathbf{v}_{i}\end{bmatrix},$ (10)

where $\mathbf{0}$ and $\mathbf{I}$ denote the zero column vector and the
identity matrix of size three, respectively. The polarizabilty tensor can be
expressed as sum of two tensors, as follow,

$\displaystyle\pmb{\alpha}$
$\displaystyle=\pmb{\alpha}_{sym}+\pmb{\alpha}_{asym}=\begin{pmatrix}\alpha_{pp}&\mathbf{0}^{T}\\\
\mathbf{0}&\pmb{\alpha}_{vv}\end{pmatrix}+\begin{pmatrix}0&\pmb{\alpha}_{pv}^{T}\\\
\pmb{\alpha}_{vp}&0\times\mathbf{I}\end{pmatrix}$ (11)

where $\pmb{\alpha}_{sym}$ denotes the tensor that only includes the direct
polarizability coefficients, and $\pmb{\alpha}_{asym}$ denotes the tensor of
polarizability arising from pure Willis coupling effect. We will utilize this
decomposition later to characterize the role of direct and Willis coupling
coefficients in determining the acoustic radiation force and radiation torque.
Subscripts $sym$ and $asym$ refers to direct and Willis coupling
polarizability hereinafter.

Considering the Green’s function identity $(\nabla^{2}+k^{2})G=-\delta(r)$,
one can write,

$\displaystyle(\nabla^{2}-\frac{1}{c_{f}^{2}}\partial_{tt})\phi_{s}\approx-\frac{j\omega}{\rho_{f}}\Omega
M\delta+\frac{j\omega}{\rho_{f}}\nabla\cdot\Big{[}\Omega\mathbf{D}\delta\Big{]},$
(12)

where $\delta$ denotes the Dirac delta impulse, and $r$ denotes the distance
from the center of the smallest sphere enclosing the scatterer. We will use
Eq.(12), which gives the approximation of the scattered field by a set of
monopole and dipole sources, to derive the acoustic radiation force and
torque.

###  II.2  Acoustic Radiation Force in the Rayleigh Limit

Using Gorkov’s far-field approach [19, 45, 2, 24], as shown in Supplementary
Notes I, the acoustic radiation force acting on a sub-wavelength particle can
be expressed as follows,

$\displaystyle\mathbf{F}$
$\displaystyle=-\rho_{f}\int_{\Omega_{\infty}}\Big{\langle}\mathbf{v}_{i}\big{[}\nabla^{2}-\frac{1}{c_{f}^{2}}\partial_{tt}\big{]}\phi_{s}\Big{\rangle}d\Omega,$
(13)

where $\Omega_{\infty}$ denotes the unbounded volume of the fluid domain,
$\langle\cdot\rangle$ denotes the time-averaging operator over one wave
period. For any two harmonically varying fields $\mathcal{F}$ and
$\mathcal{G}$, the time-averaged product is
$\langle\mathcal{F}\mathcal{G}\rangle=\frac{1}{2}\Re{\big{[}\mathcal{F}\mathcal{G}^{\ast}\big{]}}$
with $\Re$ denoting the real part of a complex quantity and $\ast$ denoting
the complex conjugation operator. By substituting Eq. (12) into (13) and
making use of the properties of the Dirac delta impulse [2, 24], the force
expression changes to

$\displaystyle\mathbf{F}$ $\displaystyle=\Big{\langle}j\omega\Omega
M\mathbf{v}_{i}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\Omega\mathbf{D}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}.$
(14)

Substituting Eq. (7) into (14), the radiation force expression expands in
terms of incident fields to

$\displaystyle\begin{split}\mathbf{F}=&\Big{\langle}j\omega\alpha_{pp}p_{i}\mathbf{v}_{i}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\pmb{\alpha}_{pv}\cdot\mathbf{v}_{i}\mathbf{v}_{i}\Big{\rangle}_{r=0}\\\
&+\Big{\langle}j\omega
p_{i}\pmb{\alpha}_{vp}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\pmb{\alpha}_{vv}\mathbf{v}_{i}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}.\end{split}$
(15)

By using $\pmb{\nabla}p_{i}=j\omega\rho_{f}\mathbf{v}_{i}$, $\langle
j\omega\mathcal{F}\mathcal{G}\rangle=-\langle\big{(}\partial_{t}\mathcal{F}\big{)}\mathcal{G}\rangle=\langle\mathcal{F}\big{(}\partial_{t}\mathcal{G}\big{)}\rangle=-\langle\mathcal{F}j\omega\mathcal{G}\rangle$
and rearranging the terms, the force expression reads

$\displaystyle\begin{split}\mathbf{F}=&-\Big{\langle}\frac{\alpha_{pp}}{\rho_{f}}p_{i}\pmb{\nabla}p_{i}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\pmb{\alpha}_{vv}\mathbf{v}_{i}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}\\\
&-\Big{\langle}\frac{1}{\rho_{f}}\pmb{\alpha}_{pv}\cdot\mathbf{v}_{i}\pmb{\nabla}p_{i}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\pmb{\alpha}_{vp}\cdot\big{[}p_{i}\pmb{\nabla}\mathbf{v}_{i}\big{]}\Big{\rangle}_{r=0}.\end{split}$
(16)

Considering, from (8) and Eq. (9) , the relation between Willis coupling
coefficients $\pmb{\alpha}_{pv}=j\omega\rho_{f}\pmb{\alpha}_{vp}$, the force
expression simplifies further to

$\displaystyle\begin{split}\mathbf{F}=&-\Big{\langle}\frac{\alpha_{pp}}{2\rho_{f}}\pmb{\nabla}\big{[}p_{i}^{2}\big{]}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\pmb{\alpha}_{vv}\mathbf{v}_{i}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}\\\
&+\Big{\langle}\frac{1}{\rho_{f}}\pmb{\alpha}_{pv}\cdot\big{[}p_{i}\pmb{\nabla}\mathbf{v}_{i}-\mathbf{v}_{i}\pmb{\nabla}p_{i}\big{]}\Big{\rangle}_{r=0}.\end{split}$
(17)

This general expression of the radiation force holds true for an object of
arbitrary shape and of sub-wavelength size.

####  II.2.1 Acoustic radiation force of a spherical object

Due to the symmetry, the Willis coupling terms $\pmb{\alpha}_{pv}$ and
$\pmb{\alpha}_{vp}$ become zero. Furthermore,
$\pmb{\alpha}_{vv}=\alpha_{vv}\mathbf{I}$ ; hence, the second term on the
right-hand side of Eq. (17) changes to

$\displaystyle\Big{\langle}j\omega\pmb{\alpha}_{vv}\mathbf{v}_{i}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}=\Big{\langle}j\omega\frac{\alpha_{vv}}{2}\pmb{\nabla}\big{[}v_{i}^{2}\big{]}\Big{\rangle}_{r=0},$
(18)

where $v_{i}$ is the magnitude of the velocity vector and
$v_{i}^{2}=\mathbf{v}_{i}\cdot\mathbf{v}_{i}$. Finally, Eq. (17) becomes

$\displaystyle\begin{split}\mathbf{F}=&-\Big{\langle}\frac{\alpha_{pp}}{2\rho_{f}}\pmb{\nabla}\big{[}p_{i}^{2}\big{]}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\frac{\alpha_{vv}}{2}\pmb{\nabla}\big{[}v_{i}^{2}\big{]}\Big{\rangle}_{r=0}.\end{split}$
(19)

Substituting Eq.(10) into Eq. (19), one could derive the Gorkov force
potential, which is the basis of radiation force potential theory [19, 2, 45,
24], $\mathbf{F}=-\pmb{\nabla}G$ with

$\displaystyle G=\frac{4\pi
a^{3}}{3}\Big{\langle}\frac{f_{1}}{2}\kappa_{f}p_{i}^{2}-\frac{3f_{2}}{4}\rho_{f}v_{i}^{2}\Big{\rangle}.$
(20)

####  II.2.2 Nonspherical scatterer with rotational and mirror symmetry

Next, we look at non-spherical geometries which still hold axial or mirror
symmetries, such as prolate (elongated) and oblate (flattened) spheroids.
Again, due to mirror symmetry, the Willis Coupling coefficients become zero.
The three terms on the diagonal of $\pmb{\alpha}_{vv}$ are no longer all equal
and Eq. (18) is invalid for an incident wave with an arbitrary 3D wavefront.
However, for incident plane waves normal to the object’s planes of symmetry,
Eq. (18) can be employed to derive the force potential since the force only
acts in the direction of incidence. And, for nonspherical scatterers with
rotational and mirror symmetry, the force expression becomes

$\displaystyle\begin{split}\mathbf{F}=&-\Big{\langle}\frac{\alpha_{pp}}{2\rho_{f}}\pmb{\nabla}\big{[}p_{i}^{2}\big{]}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\pmb{\alpha}_{vv}\mathbf{v}_{i}\cdot\pmb{\nabla}\mathbf{v}_{i}\Big{\rangle}_{r=0}.\end{split}$
(21)

The expression in Eq. (21) shows that the Gorkov force potential theory is no
longer applicable owing to
$\pmb{\alpha}_{vv}\mathbf{v}_{i}\cdot\pmb{\nabla}\mathbf{v}_{i}\neq\pmb{\alpha}_{vv}\nabla\big{[}\mathbf{v}_{i}\cdot\mathbf{v}_{i}\big{]}$,
despite the rotational or mirror symmetry of the scatterer.

####  II.2.3 Acoustic radiation force for objects of arbitrary shape

A lack of symmetry in the shape of scatterer with respect to the incident
field results in Willis coupling. Since generic shapes can always be found in
real-life, it is important to investigate the role of shape complexity and
asymmetry. For instance red and white blood cells and worm-like bacteria in
bio-acoustophoretic applications or bianisotropic metamaterials in
acoustic/photonic beam forming, wave manipulation and holography show natural
or engineered asymmetries, which are yet to be investigated in the context of
acoustic radiation force and radiation torque. Equation (17) applies to this
general case as long as the characteristic length of the object is within the
Rayleigh limit. Finally, it is evident from Eq. (17) that Gorkov potential for
acoustic radiation force is inapplicable to the general case of objects with
arbitrary shapes; hence, the force is required to be calculated directly from
Eq. (17).

###  II.3 Acoustic Radiation Torque

Acoustic radiation torque $\mathbf{T}$ is obtained from the radiation stresses
using the far-field approach, as shown in Supplementary Notes I [19, 2, 46,
27, 28], as follows,

$\displaystyle\mathbf{T}$
$\displaystyle=-\rho_{f}\int_{\Omega_{\infty}}\mathbf{x}\times\Big{\langle}\mathbf{v}_{i}\big{[}\nabla^{2}-\frac{1}{c_{f}^{2}}\partial_{tt}\big{]}\phi_{s}\Big{\rangle}d\Omega,$
(22)

where $\mathbf{x}$ denote the position vector. Substituting Eq. (12) into Eq.
(22) and using $\mathbf{x}\times\mathbf{n}=\mathbf{0}$, the torque expression
becomes

$\displaystyle\mathbf{T}$
$\displaystyle=\Big{\langle}j\omega\Omega\mathbf{D}\times\mathbf{v}_{i}\Big{\rangle}_{r=0}.$
(23)

This expression is the same as Eq. (11) in reference [47], in which the
radiation force and the radiation torque were derived from their canonical
momentum and spin densities. Finally, using Eq. (7), we obtain the radiation
torque for an arbitrarily shaped scatterer, as follows,

$\displaystyle\mathbf{T}$ $\displaystyle=\Big{\langle}j\omega
p_{i}\pmb{\alpha}_{vp}\times\mathbf{v}_{i}\Big{\rangle}_{r=0}+\Big{\langle}j\omega\big{[}\pmb{\alpha}_{vv}\mathbf{v}_{i}\big{]}\times\mathbf{v}_{i}\Big{\rangle}_{r=0}.$
(24)

This formulation of radiation torque not only gives the term corresponding to
the spin density, which is proportional to
$\big{\langle}\mathbf{v}_{i}\times\mathbf{v}_{i}\big{\rangle}$ [47], but also
shows the role of Willis coupling effects distinctively. Supplementary Notes
II outlines the details of how to calculate the polarizability tensor from the
scattering of incident standing waves in three dimensions using Boundary
Element Method.

##  III Results

###  III.1 Case of Standing Plane Wave

The incident pressure, velocity and their derivative fields are expressed as

$\displaystyle\begin{split}p_{i}&=P_{a}\cos(kz)e^{-j\omega
t},\quad\mathbf{v}_{i}=\frac{jP_{a}}{\rho_{f}c_{f}}\sin(kz)e^{-j\omega
t}\mathbf{e}_{z},\\\ \pmb{\nabla}p_{i}&=-P_{a}k\sin(kz)e^{-j\omega
t}\mathbf{e}_{z},\quad\pmb{\nabla}\mathbf{v}_{i}=\frac{jP_{a}k}{\rho_{f}c_{f}}\cos(kz)e^{-j\omega
t}\mathbf{e}_{z}\mathbf{e}_{z},\end{split}$ (25)

Without loss of generality, the propagation direction is assumed to be along
the $z$-axis, denoted by $\mathbf{e}_{z}$. Substituting Eq. (25) into Eq. (17)
and time averaging, the force and torque acting on an arbitrarily shaped
object become

$\displaystyle\begin{split}\mathbf{F}\cdot\mathbf{e}_{z}=&F_{sym}+F_{asym},\\\
F_{sym}=&\frac{kE_{i}}{\rho_{f}}\Big{(}\frac{\Re{\big{[}\alpha_{pp}\big{]}}}{\kappa_{f}}-kc_{f}\Im{\big{[}\pmb{\alpha}_{vv}\mathbf{e}_{z}\big{]}}\cdot\mathbf{e}_{z}\Big{)}\sin(2kz),\\\
F_{asym}=&\frac{kE_{i}}{\rho_{f}}c_{f}\Im{\big{[}2\pmb{\alpha}_{pv}\cdot\mathbf{e}_{z}\big{]}}\cos(2kz),\\\
\mathbf{T}=&\mathbf{T}_{sym}+\mathbf{T}_{asym},\\\
\mathbf{T}_{sym}=&-\frac{kE_{i}}{\rho_{f}}c_{f}\Im{\big{[}2\pmb{\alpha}_{vv}\mathbf{e}_{z}\big{]}}\times\mathbf{e}_{z}\sin^{2}(kz),\\\
\mathbf{T}_{asym}=&\frac{kE_{i}}{\rho_{f}}\Big{(}\frac{1}{\kappa_{f}}\Re{\big{[}\pmb{\alpha}_{vp}\big{]}}\times\mathbf{e}_{z}\Big{)}\sin(2kz),\end{split}$
(26)

where $E_{i}=\frac{P_{a}^{2}}{4\rho_{f}c_{f}^{2}}$ denotes the acoustic energy
density of the incident wave, and $\Im$ denotes the imaginary part of a
complex variable. Subscripts $sym$ and $asym$, previously defined in Eq. (11),
refer to the contributions from direct and Willis coupling coefficients,
respectively. $F_{asym}$ and $\mathbf{T}_{asym}$ denote partial force and
torque terms that arise from the Willis coupling representation of shape
complexity. The spatial dependence of partial force $F_{asym}$ is $\cos(2kz)$,
which gives the stable zero-force location with negative force gradient at
$\lambda/8$. Compared to $\mathbf{F}_{sym}$, which is classically referred to
as acoustic radiation force, with $\sin(2kz)$ leading to the prediction of
acoustic traps at pressure/velocity nodes under plane standing waves for sub-
wavelength spherical and spheroidal particles [19, 18, 31, 45, 33, 24], the
location of zero net force is shifted by up to $\lambda/8$ along the wave
direction. The effects of geometrical complexity on the primary radiation
force is the largest at pressure and velocity nodal planes, where
$\sin\big{(}2kz\big{)}=0$. Furthermore, the actual location of stable zero-
force for acoustic traps in a plane standing wave, considering the
contribution of Willis coupling, are shifted from the nodal locations as a
results of the additional force induced by the Willis Coupling effects.
However, this shift depends on how large the Willis coupling effect is. This
result implies that it is possible to obtain anomalous force and torque fields
by engineering the Willis coupling coefficients through shape manipulation.
Finally, it is noted that changing the object symmetry also changes the
$\alpha_{pp}$ and $\pmb{\alpha}_{vv}$, since a portion of the scattered energy
goes to the Willis coupling effect.

###  III.2 Case of Travelling Plane Wave

For a travelling wave in the $z$-direction, the incident pressure, velocity
and their derivative fields are expressed as

$\displaystyle\begin{split}p_{i}&=P_{a}e^{jkz}e^{-j\omega
t},\mathbf{v}_{i}=\frac{P_{a}}{\rho_{f}c_{f}}e^{jkz}e^{-j\omega
t}\mathbf{e}_{z},\\\ \pmb{\nabla}p_{i}&=-jkP_{a}e^{jkz}e^{-j\omega
t}\mathbf{e}_{z},\pmb{\nabla}\mathbf{v}_{i}=\frac{jP_{a}k}{\rho_{f}c_{f}}e^{jkz}e^{-j\omega
t}\mathbf{e}_{z}\mathbf{e}_{z},\end{split}$ (27)

Substituting Eq.(27) into (17), the force and torque under a plane travelling
wave becomes

$\displaystyle\begin{split}\mathbf{F}=&\frac{2E_{i}}{\rho_{f}}\Big{(}-\Im{\Big{[}\frac{\alpha_{pp}}{\kappa_{f}}\Big{]}}\mathbf{e}_{z}+kc_{f}\Re{\Big{[}\pmb{\alpha}_{vv}\mathbf{e}_{z}\Big{]}}\Big{)},\\\
\mathbf{T}=&-\frac{2kE_{i}}{\rho_{f}}\Big{(}\frac{1}{\kappa_{f}}\Im{\Big{[}\pmb{\alpha}_{vp}\Big{]}}\times\mathbf{e}_{z}+c_{f}\Im{\Big{[}\pmb{\alpha}_{vv}\mathbf{e}_{z}\Big{]}}\times\mathbf{e}_{z}\Big{)},\end{split}$
(28)

This expression shows that the direct contribution of the Willis Coupling
terms $\pmb{\alpha}_{pv}$ and $\pmb{\alpha}_{vp}$ to the force under traveling
wave is zero. However, the force generally has also transverse components in
$x$\- and $y$\- direction due to the $\pmb{\alpha}_{vv}\mathbf{e}_{z}$ term,
despite the incident wave’s one-dimensional propagation line.

###  III.3 Numerical results

The shapes of objects are constructed by adding a tail-like attachment
(protrusion) to or by creating a hole (intrusion) in a sphere, to engineer a
geometrical anti-symmetry in one direction. The circular tail/hole as shown in
Fig. 2 are considered to generate a non-spherical and axisymmetric shape. The
apparent symmetry is further reduced by considering the cross-section of the
attachment being a square but having the same edge length as the diameter of
the circular one; however, they both exhibit zero Willis coupling in the
normal to the length directions due to the mirror symmetry. The later design
was used to investigate the shape effects over a given range of size within
the Rayleigh limit, which in practice is reasonably $ka<0.3$ [48].

Figure 2: Schematic of non-spherical shapes for the study Willis coupling with
engineered anit-symmetry along the $z$-axis by adding, $(+)$, or removing,
$(-)$, material at one side. Square and circle cross-sections are considered,
as shown in panels (a) and (b), respectively, to reduce the spherical symmetry
to mirror symmetry in $x-$and $y$-directions, and axisymmetric respectively.
Results for the ones with circle cross section are provided in Supplemental
Notes III. For the study of size effects, the design with the tail and square
cross-section is considered for the given scaling ratios in panel (c), the
largest being 10:1 compared to the reference size.

We consider a standing plane wave and assume that objects are sound hard and
immovable, allowing us to focus only on the effects of scatterer’s exterior
shape. An incident pressure field with $10$ mm wavelength in air with
$c_{f}=343.140$ $ms^{-1}$ and $\rho_{f}=1.204$ $kgm^{-3}$ is considered. The
acoustic radiation force and torque results are normalized, as follows,

$\displaystyle\mathbf{Q}$
$\displaystyle=\frac{\mathbf{F}}{E_{i}k\Omega},\qquad\mathbf{Z}=\frac{\mathbf{T}}{E_{i}\Omega},$
(29)

where $\mathbf{Q}$ and $\mathbf{Z}$ are dimensionless vector quantities, as
shown in Fig. 1. For one dimensional problems such as a sphere in plane waves,
these quantities reduce to a scalar value, which has been referred to as force
contrast factor an torque contrast factor, respectively, and used to determine
the direction of the force and torque [18, 49, 50, 24, 21, 51, 52, 53].

(a) Real part of non-zero polarizability coefficients

(b) Imaginary part of non-zero polarizability coefficients

Figure 3: Changes of (a) real and (b) imagnary parts of the non-zero
polarizability coefficients, given in Eq. (7), with respect to the deviation
from spherical shape by adding a tail or hole with square cross-section.

To study the effect of asymmetry in the axial direction $z$, the objects are
perturbed from sphere by adding a tail or hole with square cross section. The
non-zero polarizability coefficients are shown in Fig. 3, which indicate that
the real parts of Willis coupling coefficients are negligible for the small
object sizes considered here. For shapes with $z$-symmetry, i.e. the reference
sphere and the case with a hole all the way through, both real and imaginary
parts of the Willis coupling coefficients become zero, marked by the red
circles on the horizontal axis in Fig. 3. The polarizability coefficients are
larger for the cases with the tail, added material volume ($+$), than those
with the hole, subtracted material volume ($-$). Considering Eq. (26) which
shows the relation between polarizability coefficients and the acoustic
radiation force and torque, it is expected that addition of a tail produces
larger magnitude of the torque than the cases with a hole. Finally, a
comparison between the cases with square and circular cross-sections is given
in Supplemental Notes III, showing that the polarizability coefficients are
smaller for the square cross section, for any given tail length or depth hole
in the studied range of $0$ to $2a_{s}$.

(a) Direct-polarization force $F_{sym}$

(b) Willis-coupling force $F_{asym}$

(c) Direct-polarization torque $T_{sym}$

(d) Willis-coupling torque $T_{asym}$

Figure 4: Normalized force and torque values for the case of a sphere with a
rectangular tail/hole, oriented at $\pi/4$ in the $y$-plane. Panels (a) and
(c) shows the parts corresponding to the direct polarizability coefficients,
while (b) and (d) are those arising from the Willis coupling coefficients.

To investigate the radiation force and torque further, normalized values
$Q_{x}$, $Q_{z}$ and $Z_{y}$ are shown in Fig. 4, for the object oriented at
$\pi/4$ with respect to the incidence direction, in the $y$-plane. This means
that the radiation force has the same components in the $x$ and $z$
directions, and the torque has only one component in the $y$ direction. The
spatial dependence of the force and torque components are according to Eq.
(26). The normalized values of the partial force $F_{sym}$, due to direct-
polarization coefficients, in Fig. 4(a), show that the magnitude change is
negligible compared to the reference sphere. However, it was observed that the
cases with a tail experiences almost similar forces as the reference sphere.
For those with a hole, the force increases with increasing hole depth, except
for the case of through-hole, which shows a sudden force reduction. In
contrast, the Willis-coupling force $F_{asym}$ varies more significantly as
the tail length or the depth hole increases. The only exception is the case of
through hole, $l_{h}/a_{s}=2$, which gives zero Willis-coupling force due to
$z$-symmetry. Nonetheless, the Willis-coupling force $F_{asym}$ is at least
two orders of magnitude smaller than the direct-polarization force $F_{sym}$;
hence, it could be neglected for estimating the radiation force for practical
applications.

The results of torque contrast factor $Z_{y}$, in Fig.4(c)-4(d), shows that
radiation torque is more sensitive to the deviation from spherical shape. The
changes of $T_{sym}$ are larger for the case with the tail, experiencing a
negative torque that aligns the tail with the direction of the incident wave
vector. Those with a hole are subjected to a positive torque that tends to
align the hole in the normal to the wave vector direction. The case with a
through hole is an exception as it experiences a large negative torque,
similar to the cases with a tail. These results imply that the partial torque
$T_{sym}$, due to direct polarization effect, tends to bring the object to the
orientation with the smallest cross-section normal to the incident wave
direction. Similar relation between object orientation and the radiation
torque was observed for prolate and oblate spheroids [33].

However, as shown in Fig. 4(d), these non-spherical objects are also subjected
to the Willis-coupling torque $T_{asym}$, as expressed in Eq. (26). Unlike the
large difference between partial radiation forces, this partial torque is just
smaller than the $T_{sym}$ by less than an order of magnitude. It was found
that $T_{asym}$ opposes and reinforces the $T_{sym}$ before and after the
pressure node at $z/\lambda=0.25$, respectively. Therefore, the orientation of
the objects under the action of the radiation torque needs to be determined by
accounting for both components. These observations are particularly of
interest as it captures not only the contribution of Willis coupling, but also
the significance of calculating radiation torques in capturing the essence of
shape complexity, even for sub-wavelength objects with $ka\ll 1$.

(a) Real part of non-zero polarizability coefficients

(b) Imaginary part of non-zero polarizability coefficients

Figure 5: Variation of (a) real and (b) imagnary parts of the non-zero
polarizability coefficients in terms of size for the case of adding a tail
with square cross section.

The changes of polarizability response in terms of size were investigated for
the case of a tail with square cross section and $l_{t}/a_{s}=2$. A scaling
ratio from 1:1 to 1:10 was applied, indicating an increase of size by one
order of magnitude, as shown by the $ka$ values in Fig. 5. It was found that
both real and imaginary parts of the non-zero polarization coefficients
increases with increasing size ratio. Real part of Willis-coupling
coefficients $\alpha_{pv}^{z}$ increases by at least seven orders of
magnitude, while the direct polarizability coefficients grow by three. For the
imaginary parts in Fig. 5(b), the growth of Willis-coupling is around four
orders of magnitude while it is six orders for the direct polarizability
coefficients. These changes of polarizability coefficients influences the
radiation force and torque, as can be seen in Eq. (26). Moreover, our results
of $\alpha_{pv}^{z}$ and $\alpha_{vp}^{z}$ indicate that the reciprocity
principle for the acoustic polarization, as expressed in Eq. (9), was
satisfied within the computational margin of accuracy over the given size
range.

(a) Radiation force components $F_{sym}$ and $F_{asym}$

(b) Radiation torque components $T_{sym}$ and $T_{asym}$

Figure 6: Variation of normalized radiation force and radiation torque in
terms of size for the case of a sphere with a rectangular tail, oriented at
$\pi/4$ in the $y$-plane. Solid lines indicate $F_{sym}$ and $T_{sym}$ and the
other parts are shown by dashed lines.

Inspecting Eq. (26), we expect the changes of polarizability coefficients with
size to influence the radiation force and torque, and our results are shown in
Fig. 6. Although the Willis coupling increased for larger objects, its
contribution to the radiation force is still negligible. The radiation force
is less sensitive to size increase as the spread of forces is rather narrow,
indicating almost the same magnitude for the given size range. $T_{asym}$
increased to one third of $T_{sym}$, which implies a greater contribution from
the Willis coupling for larger objects. Both parts of the radiation torque
increase faster for $ka>0.05$, which shows a nonlinear change with respect to
the size factor. These results demonstrate that the effects of size and shape
on the radiation torque are more prominent than the radiation force.

##  IV Discussion

In the presented formulation, the assumption of a sub-wavelength scatterer,
$ka<1$, led to the monopole-dipole approximation of the scattering field. For
$ka>1$, a more accurate approximation including the quadrupole and higher
order multipole moments is required to obtain the analytical expressions of
the acoustic radiation force and radiation torque. This could be achieved by
incorporating the Willis coupling factors into the partial-wave expansion
series and adjusting the scattering coefficients accordingly.

We also observed that the normalized force, also called contrast factor, shows
far less variation across the presented range of $ka$ values (Rayleigh index),
when $k\Omega$ is used instead of the choice of $\pi a^{2}$. The difference is
because $k\Omega\propto\pi a^{2}\times ka$. In previous studies mainly focused
on spheres or spheroids, the extra $ka$ factor was always applied to the
contrast factor, leading to a size dependence which bears no helpful
information [18, 49, 22, 53, 24]. Therefore, it is concluded that the
normalisation of radiation force and torque with the respect to volume, as
expressed in Eq. (29), better indicates the dependence on the shape and size
features when it comes to non-spherical objects.

The assumptions of sound-hard immovable scatterers were made to focus this
study on the effects of shape and size. However, including the effects of
material properties and wave refraction is straightforward, since these can be
readily incorporated into the multipole moments of a scatterer. Our choice of
numerical approach for solving scattering problem was Boundary Element Method,
which was combined by the analytical multipole translation and rotation to
obtain the multipole moments from the four reference cases of plane standing
wave. The alternative approach is to use Finite Element Method, which is
available in commercial software packages, and use the polarizability
retrieval techniques [44, 41], based on surface integral of scattering
pressure with appropriate spherical harmonics as weight functions.
Nevertheless, it is noted that the discretisation of the scatterer’s surface
for either of these methods requires extra care to ensure a uniform
distribution of element size and aspect ratio across the surface. Our results
were obtained after a numerical convergence test to achieve the optimal
element size. The BEM scattering results was verified by comparing against FEM
simulation of pressure on a fictitious spherical surface at $4\lambda$
distance, with relative difference of less than $1\%$, for non-spherical
objects.

The advantage of using the presented novel formulation is its ability to
calculate acoustic radiation force and acoustic radiation torque straightaway
from the polarizability tensor, as a measure of acoustic transfer function and
being independent of the incident wave field. If non-planar acoustic waves
such as Bessel or Gaussian beams are of interest, one could simply calculate
the the radiation force and torque from the incident field values at the
centroid of the object. Moreover, the effect of the angle of incidence is
easily included by rotation of the polarizability tensor. For sub-wavelength
objects, our explicit force and torque expressions can be implemented into
commercial software packages such as COMSOL for multi-physics simulations of
acoustofluidic processes. Our results show for the first time that geometrical
complexities have no effects on the radiation force that is induced by a plane
traveling wave. However, in the case of two counter-propagating waves, which
gives the plane standing wave scenario, the radiation force and radiation
torque will be influenced by such shape effects. This is an indication of the
non-linear nature of radiation force and radiation torque, for which the wave
superposition principle becomes inapplicable.

##  V Conclusion

The theory of acoustic radiation force and torque was revisited to incorporate
the Willis coupling arising from the shape complexities of an arbitrary object
in incident acoustic fields. The mathematical expression of the acoustic
radiation force and radiation torque were provided in terms of polarizability
coefficients, applicable for any choice of incident pressure field. As
examples, we derived these expressions for plane standing and travelling waves
to characterize the impact of Willis coupling and changes to the
polarizability tensor as compared to the case of simple sphere. We found that
radiation torque is more sensitive to the shape effects than radiation force.
As the size factor $ka$ increases, the contribution of the Willis coupling
effect becomes larger for the radiation torque, while it remains negligible to
the radiation force. It is concluded that knowing the acoustic radiation
torque becomes imperative when the geometry of the object is accounted for
even within the Rayleigh limit. Furthermore, there will be locations in a
standing wave with stable angular balance and zero radiation torque. The
presented work outlined the significance of the acoustic radiation force and
radiation torque in the analysis and design of acoustophoresis processes
involving non-spherical objects.

## Acknowledgement

This research is supported by Australian Research Council Discovery Project
DP200101708.

## References

  * Bruus [2011] H. Bruus, Acoustofluidics 1: Governing equations in microfluidics, Lab on a Chip 11, 3742 (2011).
  * Bruus [2012a] H. Bruus, Acoustofluidics 7: The acoustic radiation force on small particles, Lab on a Chip 12, 1014 (2012a).
  * Lenshof  _et al._ [2012] A. Lenshof, C. Magnusson, and T. Laurell, Acoustofluidics 8: Applications of acoustophoresis in continuous flow microsystems, Lab on a Chip 12, 1210 (2012).
  * Bruus [2012b] H. Bruus, Acoustofluidics 10: scaling laws in acoustophoresis, Lab on a Chip 12, 1578 (2012b).
  * Dual  _et al._ [2012] J. Dual, P. Hahn, I. Leibacher, D. Möller, T. Schwarz, and J. Wang, Acoustofluidics 19: Ultrasonic microrobotics in cavities: devices and numerical simulation, Lab on a Chip 12, 4010 (2012).
  * Wiklund  _et al._ [2012] M. Wiklund, R. Green, and M. Ohlin, Acoustofluidics 14: Applications of acoustic streaming in microfluidic devices, Lab on a Chip 12, 2438 (2012).
  * Hartono  _et al._ [2011] D. Hartono, Y. Liu, P. L. Tan, X. Y. S. Then, L.-Y. L. Yung, and K.-M. Lim, On-chip measurements of cell compressibility via acoustic radiation, Lab on a Chip 11, 4072 (2011).
  * Mohapatra  _et al._ [2018] A. R. Mohapatra, S. Sepehrirahnama, and K.-M. Lim, Experimental measurement of interparticle acoustic radiation force in the rayleigh limit, Physical Review E 97, 053105 (2018).
  * Augustsson  _et al._ [2012] P. Augustsson, C. Magnusson, M. Nordin, H. Lilja, and T. Laurell, Microfluidic, label-free enrichment of prostate cancer cells in blood based on acoustophoresis, Analytical chemistry 84, 7954 (2012).
  * Garcia-Sabaté  _et al._ [2014] A. Garcia-Sabaté, A. Castro, M. Hoyos, and R. González-Cinca, Experimental study on inter-particle acoustic forces, The Journal of the Acoustical Society of America 135, 1056 (2014).
  * Bernassau  _et al._ [2014] A. Bernassau, P. Glynne-Jones, F. Gesellchen, M. Riehle, M. Hill, and D. Cumming, Controlling acoustic streaming in an ultrasonic heptagonal tweezers with application to cell manipulation, Ultrasonics 54, 268 (2014).
  * Antfolk  _et al._ [2015] M. Antfolk, C. Magnusson, P. Augustsson, H. Lilja, and T. Laurell, Acoustofluidic, label-free separation and simultaneous concentration of rare tumor cells from white blood cells, Analytical chemistry 87, 9322 (2015).
  * Marzo  _et al._ [2015] A. Marzo, S. A. Seah, B. W. Drinkwater, D. R. Sahoo, B. Long, and S. Subramanian, Holographic acoustic elements for manipulation of levitated objects, Nature communications 6, 8661 (2015).
  * Wijaya  _et al._ [2016] F. B. Wijaya, A. R. Mohapatra, S. Sepehrirahnama, and K.-M. Lim, Coupled acoustic-shell model for experimental study of cell stiffness under acoustophoresis, Microfluidics and Nanofluidics 20, 69 (2016).
  * Marzo and Drinkwater [2019] A. Marzo and B. W. Drinkwater, Holographic acoustic tweezers, Proceedings of the National Academy of Sciences 116, 84 (2019).
  * Polychronopoulos and Memoli [2020] S. Polychronopoulos and G. Memoli, Acoustic levitation with optimized reflective metamaterials, Scientific reports 10, 1 (2020).
  * King [1934] L. V. King, On the acoustic radiation pressure on spheres, Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 147, 212 (1934).
  * Yosioka and Kawasima [1955] K. Yosioka and Y. Kawasima, Acoustic radiation pressure on a compressible sphere, Acta Acustica united with Acustica 5, 167 (1955).
  * Gorkov [1962] L. P. Gorkov, On the forces acting on a small particle in an acoustic filed in an ideal fluid, Soviet Physics - Doklady 6, 773 (1962).
  * Sepehrirahnama  _et al._ [2015a] S. Sepehrirahnama, K.-M. Lim, and F. S. Chau, Numerical analysis of the acoustic radiation force and acoustic streaming around a sphere in an acoustic standing wave, Physics Procedia 70, 80 (2015a).
  * Silva and Bruus [2014] G. T. Silva and H. Bruus, Acoustic interaction forces between small particles in an ideal fluid, Physical Review E 90, 063007 (2014).
  * Lopes  _et al._ [2016] J. H. Lopes, M. Azarpeyvand, and G. T. Silva, Acoustic interaction forces and torques acting on suspended spheres in an ideal fluid, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 63, 186 (2016).
  * Sepehrirahnama  _et al._ [2016] S. Sepehrirahnama, F. S. Chau, and K.-M. Lim, Effects of viscosity and acoustic streaming on the interparticle radiation force between rigid spheres in a standing wave, Physical Review E 93, 023307 (2016).
  * Sepehrirahnama and Lim [2020a] S. Sepehrirahnama and K.-M. Lim, Generalized potential theory for close-range acoustic interactions in the rayleigh limit, Physical Review E 102, 043307 (2020a).
  * Sepehrirahnama and Lim [2020b] S. Sepehrirahnama and K.-M. Lim, Acoustophoretic agglomeration patterns of particulate phase in a host fluid, Microfluidics and Nanofluidics 24, 1 (2020b).
  * Doinikov [2001] A. A. Doinikov, Acoustic radiation interparticle forces in a compressible fluid, Journal of Fluid Mechanics 444, 1 (2001).
  * Fan  _et al._ [2008] Z. Fan, D. Mei, K. Yang, and Z. Chen, Acoustic radiation torque on an irregularly shaped scatterer in an arbitrary sound field, The Journal of the Acoustical Society of America 124, 2727 (2008).
  * Zhang and Marston [2011] L. Zhang and P. L. Marston, Acoustic radiation torque and the conservation of angular momentum (l), The Journal of the Acoustical Society of America 129, 1679 (2011).
  * Sepehrirahnama  _et al._ [2015b] S. Sepehrirahnama, K.-M. Lim, and F. S. Chau, Numerical study of interparticle radiation force acting on rigid spheres in a standing wave, The Journal of the Acoustical Society of America 137, 2614 (2015b).
  * Doinikov [1994a] A. A. Doinikov, Acoustic radiation pressure on a compressible sphere in a viscous fluid, Journal of Fluid Mechanics 267, 1 (1994a).
  * Doinikov [1994b] A. A. Doinikov, Acoustic radiation pressure on a rigid sphere in a viscous fluid, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 447, 447 (1994b).
  * Foresti  _et al._ [2012] D. Foresti, M. Nabavi, and D. Poulikakos, On the acoustic levitation stability behaviour of spherical and ellipsoidal particles, Journal of Fluid Mechanics 709, 581 (2012).
  * Wijaya and Lim [2015] F. B. Wijaya and K.-M. Lim, Numerical calculation of acoustic radiation force and torque acting on rigid non-spherical particles, Acta Acustica united with Acustica 101, 531 (2015).
  * Mitri [2015a] F. Mitri, Acoustic radiation force on a rigid elliptical cylinder in plane (quasi) standing waves, Journal of Applied Physics 118, 214903 (2015a).
  * Wei  _et al._ [2004] W. Wei, D. B. Thiessen, and P. L. Marston, Acoustic radiation force on a compressible cylinder in a standing wave, The Journal of the Acoustical Society of America 116, 201 (2004).
  * Xie and Wei [2004] W. Xie and B. Wei, Dynamics of acoustically levitated disk samples, Physical Review E 70, 046611 (2004).
  * Garbin  _et al._ [2015] A. Garbin, I. Leibacher, P. Hahn, H. Le Ferrand, A. Studart, and J. Dual, Acoustophoresis of disk-shaped microparticles: A numerical and experimental study of acoustic radiation forces and torques, The Journal of the Acoustical Society of America 138, 2759 (2015).
  * Wijaya  _et al._ [2018] F. B. Wijaya, S. Sepehrirahnama, and K.-M. Lim, Interparticle force and torque on rigid spheroidal particles in acoustophoresis, Wave Motion 81, 28 (2018).
  * Sieck  _et al._ [2017] C. F. Sieck, A. Alù, and M. R. Haberman, Origins of willis coupling and acoustic bianisotropy in acoustic metamaterials through source-driven homogenization, Physical Review B 96, 104303 (2017).
  * Quan  _et al._ [2018] L. Quan, Y. Ra’di, D. L. Sounas, and A. Alù, Maximum willis coupling in acoustic scatterers, Physical Review Letters 120, 254301 (2018).
  * Jordaan  _et al._ [2018] J. Jordaan, S. Punzet, A. Melnikov, A. Sanches, S. Oberst, S. Marburg, and D. A. Powell, Measuring monopole and dipole polarizability of acoustic meta-atoms, Applied Physics Letters 113, 224102 (2018).
  * Melnikov  _et al._ [2019] A. Melnikov, Y. K. Chiang, L. Quan, S. Oberst, A. Alù, S. Marburg, and D. Powell, Acoustic meta-atom with experimentally verified maximum willis coupling, Nature communications 10, 1 (2019).
  * Chiang  _et al._ [2020] Y. K. Chiang, S. Oberst, A. Melnikov, L. Quan, S. Marburg, A. Alù, and D. A. Powell, Reconfigurable acoustic metagrating for high-efficiency anomalous reflection, Physical Review Applied 13, 064067 (2020).
  * Su and Norris [2018] X. Su and A. N. Norris, Retrieval method for the bianisotropic polarizability tensor of willis acoustic scatterers, Physical Review B 98, 174305 (2018).
  * Settnes and Bruus [2012] M. Settnes and H. Bruus, Forces acting on a small particle in an acoustical field in a viscous fluid, Physical Review E 85, 016327 (2012).
  * Maidanik [1958] G. Maidanik, Torques due to acoustical radiation pressure, The Journal of the Acoustical Society of America 30, 620 (1958).
  * Toftul  _et al._ [2019] I. Toftul, K. Bliokh, M. I. Petrov, and F. Nori, Acoustic radiation force and torque on small particles as measures of the canonical momentum and spin densities, Physical review letters 123, 183901 (2019).
  * Sepehrirahnama  _et al._ [2015c] S. Sepehrirahnama, F. S. Chau, and K.-M. Lim, Numerical calculation of acoustic radiation forces acting on a sphere in a viscous fluid, Physical Review E 92, 063309 (2015c).
  * Hasegawa and Yosioka [1969] T. Hasegawa and K. Yosioka, Acoustic-radiation force on a solid elastic sphere, The Journal of the Acoustical Society of America 46, 1139 (1969).
  * Marston  _et al._ [2006] P. L. Marston, W. Wei, and D. B. Thiessen, Acoustic radiation force on elliptical cylinders and spheroidal objects in low frequency standing waves, in  _AIP Conference Proceedings_ , Vol. 838 (AIP, 2006) pp. 495–499.
  * Silva [2011] G. T. Silva, Off-axis scattering of an ultrasound bessel beam by a sphere, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 58, 298 (2011).
  * Mitri [2009] F. Mitri, Acoustic radiation force of high-order bessel beam standing wave tweezers on a rigid sphere, Ultrasonics 49, 794 (2009).
  * Mitri [2015b] F. Mitri, Acoustic radiation force on oblate and prolate spheroids in bessel beams, Wave Motion 57, 231 (2015b).


# Towards a Benchmarking Suite for Kernel Tuners

Jacob O. Tørring1, Ben van Werkhoven2, Filip Petrovič3, Floris-Jan Willemsen2,
Jiří Filipovič3 and Anne C. Elster1 {jacob.torring<EMAIL_ADDRESS>{b.vanwerkhoven<EMAIL_ADDRESS>{fillo<EMAIL_ADDRESS>1 Norwegian University of Science and Technology (NTNU), Trondheim, Norway
2 Netherlands eScience Center, Amsterdam, Netherlands
3 Masaryk University, Brno, Czech Republic

###### Abstract

As computing system become more complex combining CPUs and GPUs, it is
becoming harder and harder for programmers to keep their codes optimized as
the hardware gets updated. Autotuners try to alleviate this by hiding as many
architecture-based optimization details as possible from the end-user, so that
the code can be used efficiently across different generations of systems.
Several autotuning frameworks have emerged, but a comparative analysis between
these related works is scarce, owing to the significant manual effort required
to port a tunable kernel from one tuner another.

In this article we introduce a new benchmark suite for evaluating the
performance of optimization algorithms used by modern autotuners targeting
GPUs. The suite contains tunable GPU kernels that are representative of real-
world applications, allowing for comparisons between optimization algorithms
and the examination of code optimization, search space difficulty, and
performance portability. Our framework facilitates easy integration of new
autotuners and benchmarks by defining a shared problem interface.

Our benchmark suite is evaluated based on five characteristics: convergence
rate, local minima centrality, optimal speedup, Permutation Feature Importance
(PFI), and performance portability. The results show that optimization
parameters greatly impact performance and the need for global optimization.
The importance of each parameter is consistent across GPU architectures,
however, the specific values need to be optimized for each architecture.

Our portability study highlights the crucial importance of autotuning each
application for a specific target architecture. The results reveal that simply
transferring the optimal configuration from one architecture to another can
result in a performance ranging from 58.5% to 99.9% of the optimal
performance, depending on the GPU architecture. This highlights the importance
of autotuning in modern computing systems and the value of our benchmark suite
in facilitating the study of optimization algorithms and their effectiveness
in achieving optimal performance for specific target architectures.

###### Index Terms:

autotuning, benchmarking

##  I Introduction

As computers have become more advanced in recent decades, there has been a
significant increase in their complexity. Central Processing Units (CPUs)
still form the core of modern computers, but we have seen a growing use of
accelerators like Graphics Processing Units (GPUs) and co-processors to
improve efficiency and performance. These accelerators can be highly
effective, but they can also make the task of optimizing programs for
performance increasingly difficult.

New data centers and many of the world’s top supercomputers, such as the
Top500 systems, have increasingly relied on High-Performance Computing (HPC)
systems that use a combination of CPUs and GPUs [1]. These heterogeneous
systems can make code optimization a challenging task, as the architectures of
the different components change rapidly.

To make programs run efficiently on these systems, programmers in computer
science and related fields spend considerable efforts in optimization [2].
However, as the complexity of these systems grows, the task of understanding
how all the system components interact becomes increasingly complex. To
address this, analytical methods such as compiler optimization are used to
understand the systems and turn that knowledge into rules for optimizing
programs. These rules can modify large chunks of code without changing the
semantics of the program, but they can also contain heuristics that may not
always lead to the optimal solution.

As architectures are constantly evolving and different, a heuristic-based rule
cannot always generalize well to the entire system. In these cases, empirical
methods which search for the optimal solution through trial and error can be
used, this is called autotuning.

One of the main consideration when developing such solutions is to hide as
many architecture-based optimization details as possible from the end-user, so
that the code can be used efficiently across different generations of systems.
The goal is to provide easy-to-use libraries and APIs that enable developers
to write code that runs well on different systems without having to understand
the intricacies of each system’s architecture.

In recent years, several studies have been conducted that present advancements
in optimization algorithms for autotuning [3]. Despite this, comparative
analysis between related works is scarce, owing to the significant manual
effort required to porting a tunable kernel from one tuner another. In order
to effectively study the performance of optimization algorithms for
autotuning, as well as to facilitate comparisons between optimization
algorithms implemented in different tuners, it is necessary to have a
benchmark suite that is compatible with all of these tuners.

Existing benchmark suites, such as Rodinia [4], SHOC [5], PolyBenchGPU [6],
are not suitable for this purpose as they are not tunable. They typically have
hardcoded thread block dimensions, but even worse, the code is often written
with assumptions that the block dimensions, parallelizations, and amount of
work per thread will never change. Furthermore, in many cases, the hardcoded
numbers are directly dependent on the input problems, which are often
unrealistically small for modern GPUs. As a result, tuning these codes
requires extensive modifications.

Earlier attempts at creating tunable benchmark suites [7, 8] are limited
regarding studying the effectiveness of optimization algorithms. The benchmark
suite from Petrovič et al. [7] only supports one tuning framework, while the
optimization parameters in Sund et al. [8] have limited performance impact.

The "BAT 2.0" benchmark suite aims to stimulate autotuning research by
offering a set of tunable kernels that are representative of those used in
various real-world applications, such as machine learning, image processing,
astrophysics, thermal modeling, microscopy, geospatial information systems,
and radio astronomy.

This benchmark suite facilitates for comparisons between optimization
algorithms from different autotuners by providing a standardized problem
interface for both the autotuners and benchmarks. The benchmark suite provides
general configuration space and kernel handler classes providing for easy
integration towards Optuna [9], SMAC3 [10], KernelTuner [11], KTT[7] as well
as our own basic reference tuner. This enables the study of code optimization
effectiveness, search space difficulty, performance portability, and more.
With the creation of this benchmark suite, researchers can now investigate key
questions about code optimization, search space, and performance portability.

##  II Background

GPUs were originally designed for fast processing of graphics, but they have
also been found to be effective accelerators for parallelizable AI and
computational science workloads. After their programming environments
facilitated using them as General Purpose GPUs (GPGPUs) to speed up a variety
of scientific computations, they are now a key component in many of the
world’s largest computing clusters. The optimization workload running on GPU
architectures is thus important for a range of fields beyond video game
graphics.

###  II-A Performance Portability

Relative performance111performance of a code compared to peak performance of
hardware can vary greatly across different hardware platforms, even if a
program functionally behaves the same on all of them. A problem configuration
that runs well on one GPU may perform poorly on another GPU. When writing
high-performance software, it is essential to utilize the hardware as
efficiently as possible. One solution to this is to use libraries like ATLAS
[12], which automatically tunes its configuration based on the executing
hardware. This can avoid the manual effort of finding the optimal
configuration for every hardware platform. FFTW [13] also used a similar self-
tuning approach to optimize its fast Fourier transforms.

We can measure how portable these configurations are by finding the optimal
configurations for our target platforms and then examine the relative
performance of these optimal configurations on other systems. We can thus
analyze how sensitive the configurations are to platform changes and how large
the relative performance differences are between different architectures.

###  II-B Analysing characteristics of benchmarks

####  II-B1 Feature importance

Feature importance in ML is a technique to identify influential features of a
dataset on a model’s outcome [14]. Different methods compute feature
importance, e.g. permutation importance, feature importances from tree-based
models, LASSO, etc. These methods quantify feature importance to understand a
dataset’s characteristics, identify redundant/irrelevant features, and guide
feature selection. The feature importance scores express the dataset’s
characteristics and the contribution of each feature to the model’s
performance.

In this study, we use Permutation Feature Importance (PFI) to evaluate feature
importance. PFI measures a model’s performance decrease when a feature’s
values are shuffled to understand the feature’s importance. We calculate PFI
by training a Catboost [15] Regression model on the original dataset,
shuffling each feature’s values, retraining the model, and comparing the
original and shuffled dataset’s performance metric. The PFI score for each
feature is the difference between the two. PFI helps identify important
features for a model’s performance and can detect multicollinearity among
features to prevent overfitting and unstable models.

####  II-B2 Proportion of Centrality metric

The proportion of centrality metric introduced by Schoonhoven et al. [3] is a
way to quantify the difficulty of GPU tuning. It is based on the concept of
the fitness flow graph (FFG), which contains all points in the search space
and creates a directed edge to a neighboring point if the neighbor has lower
fitness. This means that a random walk across the FFG mimics the behavior of a
randomized first-improvement local search algorithm. The expected proportion
of arrivals of each minimum then gives a metric for weighting reachability of
each minimum. The likelihood of arrival per local minima is computed using the
PageRank node centrality, which was originally used to determine the relevance
of a webpage. The PageRank values are the values of the dominant right
eigenvector of the adjacency matrix of a directed graph G, rescaled such that
each column adds up to 1. The metric is a measure of difficulty, it considers
how likely a certain subset of "suitably good" local minima are to be visited
by a local search algorithm relative to the rest. This subset is defined by
the optimal fitness and the proportion p, taking the set of nodes consisting
of local minima with fitness less than $(1+p)f_{opt}$ for minimization
problems, otherwise $(1-p)f_{opt}$.

##  III Related work

In the field of autotuning, the most relevant prior work includes the
benchmark suite developed by Petrovič et al. [7], Polybench-GPU [6], and Sund
et al. [8]. However, these benchmark suites have issues that limit their
usefulness. The benchmark in PolyBench-GPU has small search spaces, ranging
from 116 to 725 different possible configurations, which would fall below our
threshold for an interesting autotuning study as many real-world applications
have much larger search spaces. The benchmark suite from Petrovič et al. only
supports a single autotuning framework, while the optimization parameters in
the benchmark suite by Sund et al. have limited performance impact. To address
this issue, we have selected benchmarks with larger and more interesting
search spaces for this new version of BAT, and also added several new
benchmarks that meet our criteria. These benchmarks all give significant
speedups, with the performance of the optimal configurations significantly
varying between different systems.

For other related works, the Collective Knowledge framework (CK) developed by
Fursin et al. [16] offers a more generalized approach to benchmarking and
reproducibility. In the area of Hyperparameter Optimization (HPO) several
benchmark suites have been developed, such as the high-dimensional HPO
benchmark suite [17], LassoBench [18], and HPOBench [19].

##  IV Benchmarks

###  IV-A GEMM

Generalized dense matrix-matrix multiplication (GEMM) is part of the BLAS
linear algebra specification, and is one of the most widely-used GPU kernels.
The GEMM kernel included in BAT is from CLBlast [20], a tunable OpenCL BLAS
library. GEMM implements the multiplication of two matrices, $A$ and $B$:

$C=\alpha A\cdot B+\beta C$

where $\alpha$ and $\beta$ are scalars and $C$ is the output matrix. The
CLBlast GEMM kernel is tunable with the parameters shown in Table I. `MWG` and
`NWG` control the amount of work assigned to each thread block. `MDIMC` and
`NDIMC` describe the size of thread block, while `MDIMA` and `MDIMB` control
shared memory usage, `VWM` and `VWN` are the vector widths used for loading
from and storing to global memory, and `SA` and `SB` enables or disables the
use of shared memory for elements in $A$ and $B$.

Table I: Tunable parameters – GEMM kernel in BAT. Parameter | Values | #
---|---|---
`MWG` |  $\\{$16, 32, 64, 128$\\}$ | 4
`NWG` |  $\\{$16, 32, 64, 128$\\}$ | 4
`MDIMC` |  $\\{$8, 16, 32$\\}$ | 3
`NDIMC` |  $\\{$8, 16, 32$\\}$ | 3
`MDIMA` |  $\\{$8, 16, 32$\\}$ | 3
`NDIMB` |  $\\{$8, 16, 32$\\}$ | 3
`VWM` |  $\\{$1, 2, 4, 8$\\}$ | 4
`VWN` |  $\\{$1, 2, 4, 8$\\}$ | 4
`SA` |  $\\{$0, 1$\\}$ | 2
`SB` |  $\\{$0, 1$\\}$ | 2

###  IV-B N-body

The N-body kernel computes gravitational forces between N bodies, typically
applied in astrophysical simulations. The N-body kernel in BAT was created by
Petrovič et al. for use in KTT [7], as a tunable implementation of the code
sample from the CUDA SDK. The N-body kernel follows a simple quadratic scheme
where the forces between all pairs of bodies are computed every iteration. As
such, the kernel is very compute intensive.

The tunable parameters for the N-body kernel in BAT are shown in Table II. The
inner loop unroll factor parameters determine the degree to which partial loop
unrolling is applied for various loops in the kernel. The
`outer_unroll_factor` controls the amount of work allocated to each thread.
The `use_soa` parameter specifies whether the input bodies are stored in an
array of structures or a structure of arrays. `local_mem` enables or disables
the use of shared memory as a software managed cache. `vector_type` is to
control the number of elements loaded from memory in one instruction.

Table II: Tunable parameters – Nbody kernel in BAT. Parameter |  Values | #
---|---|---
`block_size` |  $\\{$64, 128, 256, 512$\\}$ | 4
`outer_unroll_factor` |  $\\{$1, 2, 4, 8$\\}$ | 4
`inner_unroll_factor1` |  $\\{$0, 1, 2, 4, 8, 16, 32$\\}$ | 7
`inner_unroll_factor2` |  $\\{$0, 1, 2, 4, 8, 16, 32$\\}$ | 7
`use_soa` |  $\\{$0, 1$\\}$ | 2
`local_mem` |  $\\{$0, 1$\\}$ | 2
`vector_type` |  $\\{$1, 2, 4$\\}$ | 3

###  IV-C Hotspot

The Hotspot kernel included in BAT is based on the Hotspot kernel in the
Rodinia Benchmark suite [4]. The kernel is part of a thermal simulation
application used to estimate processor temperature based on processor
architecture and simulated power currents. The kernel iteratively solves a
series of differential equations. The kernel inputs are the power and initial
temperatures, the output is a grid of average temperature values spanning the
chip.

To simplify the indexing scheme and increase the tunability of the kernel, we
have re-implemented the Hotspot kernel in Rodinia from scratch. The main
difference of our implementation with that of Rodinia is that our kernel can
be used with any thread block dimension, can arbitrarily vary the amount of
work per thread, and vary the extent to which temporal tiling is applied.

The tunable parameters for the Hotspot kernel in BAT are shown in Table III.
`block_size_x` and `block_size_y` describe the thread block dimensions in x
and y, the kernel uses at least 32 and at most 1024 threads. `tile_size_x` and
`tile_size_y` control the number of output elements computed by each thread in
the x and y dimensions. `temporal_tiling_factor` is the number of iterations
of the stencil operation performed by a single kernel launch, for more details
on the temporal tiling optimization see Hijma et al. [2]. `sh_power` enables
or disables the use of shared memory as a cache for storing the input power
currents. `blocks_per_sm` is used in the `__launch__bounds()` directive in
CUDA to hint the compiler to aim for a certain occupancy when running the
kernel, effectively this optimization encourages the compiler to decrease
register usage in the kernel.

Table III: Tunable parameters – Hotspot kernel in BAT. Parameter |  Values | #
---|---|---
`block_size_x` |  $\\{1,2,4,8,32n\mid 32n\in[32,1024]\\}$ | 37
`block_size_y` |  $\\{$1, 2, 4, 8, 16, 32$\\}$ | 6
`tile_size_x` |  $\\{$n $\mid$ n $\in$ [1, 10] $\\}$ | 10
`tile_size_y` |  $\\{$n $\mid$ n $\in$ [1, 10] $\\}$ | 10
`temporal_tiling_factor` |  $\\{$n $\mid$ n $\in$ [1, 10] $\\}$ | 10
`loop_unroll_factor_t` |  $\\{$n $\mid$ n $\in$ [1, 10] $\\}$ | 10
`sh_power` |  $\\{$0, 1$\\}$ | 2
`blocks_per_sm` |  $\\{$0, 1, 2, 3, 4$\\}$ | 5

###  IV-D Pnpoly

Pnpoly (Point-in-polygon) kernel is used by Goncalves et al. [21] as part of a
geospatial database system for massive point clouds obtained through airborne
LiDAR. The kernel is used to query all points within a certain outline, for
example points on highways or all points within a city. Pnpoly has been used
as a benchmark kernel for autotuning in several studies [22, 3]. However, the
Pnpoly kernel in BAT includes only the GPU kernel of the full GPU-enabled
database operator.

The tunable parameters of the Pnpoly kernel in BAT are listed in Table IV.
`block_size_x` is simply the number of threads per block. `tile_size` the
amount of points processed by each thread. `between_method` selects the
algorithm to use to see if a point lies between two other points. Similarly,
`use_method` selects the algorithm that is used to keep track of whether the
evaluated point is inside or outside of the polygon.

Table IV: Tunable parameters – Pnpoly kernel in BAT. Parameter |  Values | #
---|---|---
`block_size_x` |  $\\{32n|32n\in[32,1024]\\}$ | 31
`tile_size` |  $\\{1,2n|2n\in[2,20]\\}$ | 11
`between_method` |  $\\{0,1,2,3\\}$ | 4
`use_method` |  $\\{0,1,2\\}$ | 3

###  IV-E Convolution

Van Werkhoven et al. [23] have implemented an optimized and highly-tunable
GPU-accelerated library for 2D Convolution operations, which has become a
commonly used benchmark in autotuning [24, 7, 25, 3].

A convolution operation computes a linear combination of weights and a range
of the input image for each output pixel. A 2D convolution of an input image
$I$ of size $(w\times h)$ and a convolution filter $F$ of size $(F_{w}\times
F_{h})$ computes an output image $O$ of size $((w-F_{w})\times(h-F_{h}))$:

$O(x,y)=\sum\limits_{j=0}^{F_{h}}\sum\limits_{i=0}^{F_{w}}I(x+i,y+j)\times
F(i,j)$

The tunable parameters of the Convolution kernel in BAT are listed in Table V.
`block_size_x` and `block_size_y` describe the thread block dimensions,
`tile_size_x` and `tile_size_y` the number of output pixels processed by each
thread in the x and y dimensions. `use_padding` controls whether or not to use
the padding scheme in shared memory that is used to avoid shared memory bank
conflicts as described in Van Werkhoven et al. [23]. Padding is only
significant when the `block_size_x` is not a multiple of number of memory
banks in shared memory. Finally, `read_only` controls whether or not to load
input elements from global memory through read-only cache.

Table V: Tunable parameters – Convolution kernel in BAT. Parameter |  Values | #
---|---|---
`block_size_x` |  $\\{$1, 2, 4, 8, 16, 32, 48, 64, 80, 96, 112, 128$\\}$ | 12
`block_size_y` |  $\\{1,2,4,8,16,32\\}$ | 6
`tile_size_x` |  $\\{1,2,3,4,5,6,7,8\\}$ | 8
`tile_size_y` |  $\\{1,2,3,4,5,6,7,8\\}$ | 8
`use_padding` |  $\\{0,1\\}$ | 2
`read_only` |  $\\{0,1\\}$ | 2

###  IV-F Expdist

The Expdist kernel is part of a localization microscopy applications that
implements a template-free particle fusion algorithm by combining many
different observations into a single super-resolution reconstruction [26]. The
expdist kernel is used as part of the registration process where the kernel is
called repeatedly to quantify the registration of two particles. The distance
between two particles $t$ and $m$, given registration $M$, is computed as
follows:

$D=\sum\limits_{i=1}^{K_{t}}\sum\limits_{j=1}^{K_{m}}\textrm{exp}\left(-\frac{\|\vec{x}_{t,i}-M(\vec{x}_{m,j})\|^{2}}{2\sigma^{2}}\right)$

The kernel operates directly on the individual localizations ($\vec{x_{t}}$
and $\vec{x_{m}}$) in each particle rather than pixelated images and takes the
uncertainties in the localizations ($\sigma$) into account. The algorithm is
quadratic in the number of localizations per particle and is as such very
compute intensive.

The tunable parameters used in the ExpDist kernel are shown in Table VI. The
kernel supports two main implementations that are controlled by the
`use_column` parameter. When `use_column` is set to 1, the kernel reduces the
number of thread blocks used to perform the computation by using a fixed
number of thread blocks in the y dimension, set by `n_y_blocks`.
`use_shared_mem` use shared memory selects the way in which shared memory is
used.

Table VI: Tunable parameters – ExpDist kernel in BAT. Parameter |  Values | #
---|---|---
`block_size_x` |  {32, 64, 128, 256, 512, 1024} | 6
`block_size_y` |  {1, 2, 4, 8, 16, 32} | 6
`tile_size_x` |  {1, 2, 3, 4, 5, 6, 7, 8} | 8
`tile_size_y` |  {1, 2, 3, 4, 5, 6, 7, 8} | 8
`use_shared_mem` |  {0, 1, 2} | 3
`loop_unroll_factor_x` |  {1, 2, 3, 4, 5, 6, 7, 8} | 8
`loop_unroll_factor_y` |  {1, 2, 3, 4, 5, 6, 7, 8} | 8
`use_column` |  {0, 1} | 2
`n_y_blocks` |  {1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024} | 11

###  IV-G Dedispersion

The Dedispersion kernel in BAT originates from the AMBER pipeline for the
detection of single pulse astronomical transients [27]. Dedispersion is the
process of reverting the dispersion of a radio signal transmitted over many
frequencies through space. The signal component with the highest frequency
$f_{h}$ is received at time $t_{x}$, while simultaneously emitted components
with lower frequency arrive at $t_{x}+k$, where $k$ is the delay in seconds as
by the dispersion equation:

$k\approx 4150\times
DM\times\left(\frac{1}{f_{i}^{2}}\times\frac{1}{f_{h}^{2}}\right)$

The kernel takes samples in time across many frequency bands (channels) as
input and outputs the dedispersed samples for many different dispersion
measure $DM$ values. The kernel is parallelized such that each thread can work
on multiple samples and dispersion measures, while iterating over the
frequency bands. As input for the BAT Dedispersion kernel, we are using the
parameters from the ARTS survey on the Apertif telescope [28], which uses a
sampling rate of 24.4 KHz, 2048 DMs, and 1536 channels.

The tunable parameters of the dedispersion kernel are shown in Table VII. The
`loop_unroll_factor_channel` parameter depends on the input, as any divisor of
the number of channels can be used as a partial loop unrolling factor for the
inner loop in the kernel. When the loop unroll factor is 0, it is left to the
CUDA compiler to decide whether or not to apply loop unrolling.
`tile_stride_x` controls the stride used to vary the amount of work per
threads. When `tile_stride_x` is 0 and `tile_size_x` is larger than 1, threads
will process `tile_size_x` consecutive samples, when `tile_stride_x` is 1
threads will process `tile_size_x` samples that are each `block_size_x` apart
in the input. `tile_stride_y` works similarly but for dispersion measures in
the y-dimension.

Table VII: Tunable parameters – Dedispersion kernel in BAT. Parameter |  Values | #
---|---|---
`block_size_y` |  $\\{1,2,4,8,16n\mid 16n\in[16,512]\\}$ | 36
`block_size_y` |  $\\{4n\mid 4n\in[4,128]\\}$ | 32
`tile_size_x` |  $\\{n|n\in[1,16]\\}$ | 16
`tile_size_y` |  $\\{n|n\in[1,16]\\}$ | 16
`tile_stride_x` |  {0, 1} | 2
`tile_stride_y` |  {0, 1} | 2
`loop_unroll_factor_channel` |  {0, 1, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256, 384, 512, 768, 1536} | 21
`blocks_per_sm` |  {0, 1, 2, 3, 4} | 5

##  V Experimental Design

###  V-A Benchmarks, Hardware and Runtime Environment

We ran our benchmarks on four different systems with four different Nvidia
GPUs. These GPUs include the RTX 2080Ti, RTX 3060, RTX 3090 and RTX Titan. For
the Pnpoly, Nbody, GEMM and Convolution benchmarks we performed an exhaustive
search of the entire search space. For the Hotspot, Dedisp and Expdist
benchmarks our results are based on 10 000 random configurations from the
search space for each architecture.

##  VI Results and Discussion

###  VI-A Distribution of configurations

(a) GEMM

(b) Hotspot

(c) Dedisp

(d) Convolution

(e) Pnpoly

(f) Nbody

(g) Expdist

Figure 1: Performance distribution of configuration for all benchmarks on all
architectures

Fig 1 shows the distribution of configurations centered around the median
performing configuration. The plot extends from the worst to the best
configuration. The first thing to observe is that the distribution shapes are
significantly different between the different benchmarks, but similar in shape
across GPUs (this agrees with results observed in other benchmark sets [29]).
Most of the benchmarks have a high density of configurations around the median
and then exponential decay toward the best-performing configurations. The
Hotspot benchmark in Fig. 1(b) is an outlier among the benchmarks, with a high
density around the median configurations, but a cluster of very highly
performing configurations giving more than 10x speedup. The Nbody distribution
in Fig. 1(f) also shows a distinct high-density cluster of configurations that
perform very poorly.

###  VI-B Convergence towards optimum using Random Search

Fig. 2 shows the convergence to the optimum configuration with relative
performance (y-axis) plotted against the number of function evaluations on a
Symmetric Log scale (linear scale from 0 to 1). Results are from random
sampling 100 times from exhaustive or partial runs, with the median of the
best evaluation plotted.

We can observe that there is a significant variance in the convergence between
benchmarks, while there is less difference between GPUs. Again the Hotspot
benchmark in Fig. 2(b) is a clear outlier, with Random Search quickly
approaching a performance that is close to optimal. We stipulate that this is
the due to the size of the high performing cluster shown in Section VI-A. This
cluster is likely of a sufficient size such that random search can quickly
find a solution in this cluster, which is then close to optimal.

There are also significant differences in how quickly the other benchmarks
converge towards the optimal. Expdist in Fig. 2(g) and Nbody in Fig. 2(f)
achieve a 90% optimum performance after just 10 function evaluations. For
Dedisp (Fig. 2(c)) and PnPoly (Fig. 2(e)) it takes around 100 evaluations to
reach the same level. We can also see here how for the RTX Titan the PnPoly
benchmark shows how a single highly performing configuration can be the source
of the final jump in relative performance.

Lastly Convolution (Fig. 2(d)) and GEMM (Fig. 2(a)) require hundreds of
configurations to exceed 90%.

(a) GEMM

(b) Hotspot

(c) Dedisp

(d) Convolution

(e) PnPoly

(f) Nbody

(g) Expdist

Figure 2: Convergence towards optimum for all benchmarks on all architectures

###  VI-C Proportion of centrality

We are using the proportion of centrality metric proposed by Schoonhoven et
al. [3] to calculate the search difficulty of the benchmarks.We did not have
sufficient resources to calculate the metric for the benchmarks with the
largest search spaces, incl. Hotspot, Dedisp and Expdist. The results are
shown in Figure 3. The results indicate that local search algorithms will
generally find better performing configurations on the Convolution benchmark
compared with GEMM and Pnpoly, which are comparatively more difficult
benchmarks under this metric. This is in contrast to the results from our
Random Search results, where all three benchmarks have similar trajectories
towards the optimum.

(a) GEMM

(b) Convolution

(c) Pnpoly

Figure 3: Proportion of centrality for all benchmarks on all architectures

###  VI-D Max speedup over Median

Fig. 4 shows the speedup between the Median performance configuration of the
search space and the best possible configuration found. While most of the
benchmarks have speedups between 1.5 - 3.06x, outliers like the Hotspot
benchmark have very significant speedups from 11.12 - 11.97x.

Figure 4: Max speedup over Median configuration

###  VI-E Performance Portability

We can analyze the performance portability of configurations to measure the
degree to which configurations are specifically optimized for each
architecture. In Fig. 5 we show the relative performance compared with the
optimal configuration for each architecture as the optimal configurations are
transferred to the other architectures. The direction for this transfer is
described by the optimal configuration for the GPU labeled in each row, being
transferred to the different GPUs labeled on each of the columns.

In Fig. 5(b) we plot the portability for PnPoly. This shows that
configurations are very portable between the RTX 3060 and RTX 3090, however
configurations optimized for the RTX 3090 transfer poorly to the RTX Titan
(58.5% of optimal) and the 2080Ti (67.1%). Similarly for the Convolution
benchmark in Fig. 5(a) the optimal configuration for the RTX 3060 transfers
poorly to the RTX 2080Ti (73.3%) and RTX Titan (75.0%).

(a) Convolution

(b) Pnpoly

(c) Nbody

Figure 5: Performance portability for exhaustively searched benchmarks on all
architectures

(a) GEMM

(b) Nbody

(c) Hotspot

(d) Convolution

(e) Pnpoly

(f) Dedisp

(g) Expdist

Figure 6: Feature importance for all benchmarks on all architectures

###  VI-F Feature importance

To assess the importance of the different search parameters for the
benchmark’s objective, we train a Catboost Regression model over the dataset
and analyze this model to investigate which features it finds useful for
prediction. We use Permutation Feature Importance to then analyze the impact
of each individual parameter on the model’s predictive power.

Training the model on the datasets the majority of the benchmarks and GPUs
that use our CatBoost model is able to predict the performance of different
configurations very precisely, with an R-squared score from 0.992 and upwards
for all benchmarks except Convolution, where it ranges from 0.9268 to 0.9361.

Using this model we can then generate the feature importances in Fig. 6. We
can observe that for many of the benchmarks, especially GEMM in Fig. 6(a) and
Nbody in Fig. 6(b), many of the parameters do not appear to have any
meaningful impact on the model’s predictive performance. Although the
significance of these features may not extend to GPUs that are vastly
different from those tested in this study, our findings are generally
consistent across the various GPUs examined.

###  VI-G Feature importance impact on relevant search space

Table VIII: Search space sizes of benchmarks in BAT Benchmark | Cardinality | Constrained | Valid | Reduced | Reduce-Constrained
---|---|---|---|---|---
PnPoly | 4 092 | 4 092 | 3 734 - 3 774 | 4 092 | 3 734 - 3 774
Nbody | 9 408 | 1 568 | 1 568 | 112 | 70
Convolution | 18 432 | 9 400 | 5 220 - 5 256 | 4 700 | 4 700
GEMM | 82 944 | 17 956 | 17 956 | 17 956 | 17 956
Expdist | 9 732 096 | 540 000 | N/A | 144 | 96
Hotspot | 22 200 000 | 21 850 147 | N/A | 220 000 | 202 582
Dedisp | 123 863 040 | 107 011 905 | N/A | 3 870 720 | 3 327 135

Given our previous analysis we can reduce the search space of the benchmarks
to only include those parameters that has at least 0.05 feature importance on
any of the architectures. The results for the cardinality of the search spaces
can be seen in Table VIII. This gives an indication towards the size of the
most interesting parts of the search space for the tested GPUs. Researchers
using this benchmark can search over the full search space, but use this
information to get better insight into how their models are able to
concentrate on the most interesting parts of the search space.

###  VI-H Discussion

The results indicate the optimization parameters in our benchmarks have a
significant impact on performance. While some parameters have more impact than
others, the act of optimizing some parameters interact with other parameters.
We observe this behavior through the Permutation Feature Importance summing up
to a value much greater than 1 for many of the benchmarks. This behavior only
occurs when there are significant dependencies between features. Thus this
provides evidence towards the need for global optimization as opposed to
orthogonal search algorithms. While the importance of each parameter is
generally consistent across GPU architectures, the specific values for these
parameters need to be optimized for the target architecture. Our portability
study shows that simply transferring the optimal configuration from one
architecture to another can give as low as 58.5% of the optimal performance,
while other configurations can be ported at 99.9% of the optimal performance.
Generally this is the case for GPUs of the same family of architectures like
our RTX 3060 and RTX 3090.

##  VII Conclusion and Future work

HPC systems with GPU are becoming increasingly more complex and challenging to
hand-tune codes for. Autotuning framworks provide means that parameterizes
kernels for a range of system parameters. BAT 2.0, the new benchmarking suite
introduced in this study, provides a comprehensive framework for evaluating
the performance of optimization algorithms in modern computing systems
utilizing GPUs. The results of our analysis reveal that the optimization
parameters have a significant impact on performance and the need for global
optimization. The importance of autotuning is highlighted in the portability
study, which shows that optimal performance can only be achieved by optimizing
each application for a specific target architecture. The benchmarking suite
facilitates the study of optimization algorithms and their effectiveness in
achieving optimal performance, positioning it as a valuable tool in modern
autotuning research.

A C++-based interface towards BAT is currently under development to support
tuners like KTT [30] and CLTune [24]. Future work should thus include a
comparison between C++-based and Python-based kernel tuners.

## Acknowledgements

The authors would like to acknowledge NTNU for the PhD stipend and support of
our HPC-Lab that facilitated the development of this project. The authors
would also like to thank Asst. Prof. Fredrik Kjolstad and his compiler group,
at Stanford University, for his valuable comments and for hosting the first
author during the final period of writing this article. We would also like to
acknowledge PhD Candidate Tor Andre Haugdahl at NTNU for their valuable
feedback. The sixth author would also like to thank the Center for Geophysical
Forecasting (SFI CGF) at NTNU and RCN (NFR proj. no. 309960) for their support
during this project.

## References

  * [1] Stijn Heldens et al.  “The landscape of exascale research: A data-driven literature analysis”  In _ACM Computing Surveys (CSUR)_ 53.2 ACM New York, NY, USA, 2020, pp. 1–43
  * [2] Pieter Hijma et al.  “Optimization Techniques for GPU Programming”  In _ACM Comput. Surv._ New York, NY, USA: Association for Computing Machinery, 2022  DOI: 10.1145/3570638
  * [3] Richard Schoonhoven, Ben Werkhoven and K Joost Batenburg  “Benchmarking optimization algorithms for auto-tuning GPU kernels”  In _IEEE Transactions on Evolutionary Computation_ IEEE, 2022
  * [4] Shuai Che et al.  “Rodinia: A benchmark suite for heterogeneous computing”  In _2009 IEEE International Symposium on Workload Characterization (IISWC)_ Austin, TX, USA: IEEE, 2009, pp. 44–54  DOI: 10.1109/IISWC.2009.5306797
  * [5] Anthony Danalis et al.  “The Scalable Heterogeneous Computing (SHOC) benchmark suite”  In _Proceedings of the 3rd Workshop on General-Purpose Computation on Graphics Processing Units_ , GPGPU-3  Pittsburgh, Pennsylvania, USA: Association for Computing Machinery, 2010, pp. 63–74  DOI: 10.1145/1735688.1735702
  * [6] Scott Grauer-Gray et al.  “Auto-tuning a high-level language targeted to GPU codes”  In _2012 Innovative Parallel Computing (InPar)_ San Jose, CA, USA: IEEE, 2012, pp. 1–10  DOI: 10.1109/InPar.2012.6339595
  * [7] Filip Petrovič et al.  “A Benchmark Set of Highly-efficient CUDA and OpenCL Kernels and its Dynamic Autotuning with Kernel Tuning Toolkit” arXiv: 1910.08498  In _arXiv:1910.08498 [cs]_ , 2019  URL: http://arxiv.org/abs/1910.08498
  * [8] Ingunn Sund, Knut A. Kirkhorn, Jacob O. Tørring and Anne C. Elster  “BAT: A Benchmark suite for AutoTuners” Number: 1  In _Norsk IKT-konferanse for forskning og utdanning_ , 2021, pp. 44–57  URL: https://ojs.bibsys.no/index.php/NIK/article/view/915
  * [9] Takuya Akiba et al.  “Optuna: A Next-generation Hyperparameter Optimization Framework”  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, KDD ’19  Anchorage, AK, USA: Association for Computing Machinery, 2019, pp. 2623–2631  DOI: 10.1145/3292500.3330701
  * [10] Frank Hutter, Holger H. Hoos and Kevin Leyton-Brown  “Sequential Model-Based Optimization for General Algorithm Configuration”  In _Learning and Intelligent Optimization_ , Lecture Notes in Computer Science  Berlin, Heidelberg: Springer, 2011, pp. 507–523  DOI: 10.1007/978-3-642-25566-3_40
  * [11] Ben Werkhoven  “Kernel Tuner: A search-optimizing GPU code auto-tuner”  In _Future Generation Computer Systems_ 90, 2019, pp. 347–358  DOI: 10.1016/j.future.2018.08.004
  * [12] R. Clint Whaley, Antoine Petitet and Jack J. Dongarra  “Automated empirical optimizations of software and the ATLAS project”  In _Parallel Computing_ 27.1, New Trends in High Performance Computing, 2001, pp. 3–35  DOI: 10.1016/S0167-8191(00)00087-9
  * [13] M. Frigo and S.G. Johnson  “FFTW: an adaptive software architecture for the FFT” ISSN: 1520-6149  In _Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP ’98 (Cat. No.98CH36181)_ 3, 1998, pp. 1381–1384 vol.3  DOI: 10.1109/ICASSP.1998.681704
  * [14] Alexander Zien, Nicole Krämer, Sören Sonnenburg and Gunnar Rätsch  “The Feature Importance Ranking Measure”  In _Machine Learning and Knowledge Discovery in Databases_ , Lecture Notes in Computer Science  Berlin, Heidelberg: Springer, 2009, pp. 694–709  DOI: 10.1007/978-3-642-04174-7_45
  * [15] Liudmila Prokhorenkova et al.  “CatBoost: unbiased boosting with categorical features”  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18  Red Hook, NY, USA: Curran Associates Inc., 2018, pp. 6639–6649
  * [16] Grigori Fursin  “Collective knowledge: organizing research projects as a database of reusable components and portable workflows with common interfaces”  In _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 379.2197, 2021, pp. rsta.2020.021120200211  DOI: 10.1098/rsta.2020.0211
  * [17] Zhen Wang et al.  “FedHPO-B: A Benchmark Suite for Federated Hyperparameter Optimization” arXiv:2206.03966 [cs]  arXiv, 2022  URL: http://arxiv.org/abs/2206.03966
  * [18] Kenan Šehić, Alexandre Gramfort, Joseph Salmon and Luigi Nardi  “LassoBench: A High-Dimensional Hyperparameter Optimization Benchmark Suite for Lasso” arXiv:2111.02790 [cs]  arXiv, 2022  URL: http://arxiv.org/abs/2111.02790
  * [19] Katharina Eggensperger et al.  “HPOBench: A Collection of Reproducible Multi-Fidelity Benchmark Problems for HPO” arXiv:2109.06716 [cs]  arXiv, 2022  DOI: 10.48550/arXiv.2109.06716
  * [20] C. Nugteren  “CLBlast: A tuned OpenCL BLAS library”  In _Proceedings of the International Workshop on OpenCL_ , IWOCL ’18  ACM, 2018, pp. 5:1–5:10
  * [21] Romulo Goncalves et al.  “A spatial column-store to triangulate the Netherlands on the fly.”  In _Proceedings of the 24th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems_ , 2016, pp. 1–4
  * [22] Floris-Jan Willemsen, Rob Nieuwpoort and Ben Werkhoven  “Bayesian Optimization for auto-tuning GPU kernels”  In _2021 International Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems (PMBS)_ , 2021, pp. 106–117  DOI: 10.1109/PMBS54543.2021.00017
  * [23] B. Werkhoven, J. Maassen and H.. Bal  “Optimizing convolution operations on GPUs using adaptive tiling”  In _Future Gener. Comput. Syst._ 30 Elsevier, 2014, pp. 14–26
  * [24] Cedric Nugteren and Valeriu Codreanu  “CLTune: A Generic Auto-Tuner for OpenCL Kernels”  In _2015 IEEE 9th International Symposium on Embedded Multicore/Many-core Systems-on-Chip_ Turin, Italy: IEEE, 2015, pp. 195–202  DOI: 10.1109/MCSoC.2015.10
  * [25] Ben Werkhoven  “Kernel Tuner: A search-optimizing GPU code auto-tuner”  In _Future Generation Computer Systems_ 90, 2019, pp. 347–358  DOI: https://doi.org/10.1016/j.future.2018.08.004
  * [26] Hamidreza Heydarian et al.  “Template-free 2D particle fusion in localization microscopy”  In _Nature methods_ 15.10 Nature Publishing Group, 2018, pp. 781–784
  * [27] Alessio Sclocco, Stijn Heldens and Ben Werkhoven  “AMBER: a real-time pipeline for the detection of single pulse astronomical transients”  In _SoftwareX_ 12 Elsevier, 2020, pp. 100549
  * [28] WA Cappellen et al.  “Apertif: Phased array feeds for the Westerbork Synthesis Radio Telescope-System overview and performance characteristics”  In _Astronomy & Astrophysics_ 658 EDP Sciences, 2022, pp. A146
  * [29] Jaroslav Ol’ha, Jana Hozzová, Jan Fousek and Jiřı́ Filipovič  “Exploiting historical data: Pruning autotuning spaces and estimating the number of tuning steps”  In _Concurrency and Computation: Practice and Experience_ 32.21 Wiley Online Library, 2020, pp. e5962
  * [30] Jiří Filipovič, Filip Petrovič and Siegfried Benkner  “Autotuning of OpenCL Kernels with Global Optimizations” event-place: Portland, OR, USA  In _Proceedings of the 1st Workshop on AutotuniNg and aDaptivity AppRoaches for Energy Efficient HPC Systems_ , ANDARE ’17  New York, NY, USA: ACM, 2017, pp. 2:1–2:6  DOI: 10.1145/3152821.3152877