11institutetext: Key Laboratory of Solar Activities of National Astronomical
Observatories, Chinese Academy of Sciences,Beijing 100012, China;
22institutetext: State Key Laboratory of Space Weather, Chinese Academy of
Sciences, Beijing 100080, China

# Monte Carlo simulations of a diffusive shock with multiple scattering
angular distributions

Xin Wang NAOC fellow: . This work was funded in part by CAS-NSFC grant
10778605 and NSFC grant 10921303 and the National Basic Research Program of
the MOST (Grant No. 2011CB811401).1122<EMAIL_ADDRESS>Yihua Yan 11

###### Abstract

Aims. We independently develop a simulation code following the previous
dynamical Monte Carlo simulation of the diffusive shock acceleration under the
isotropic scattering law during the scattering process, and the same results
are obtained.

Methods. Since the same results test the validity of the dynamical Monte Carlo
method for simulating a collisionless shock, we extend the simulation toward
including an anisotropic scattering law for further developing this dynamical
Monte Carlo simulation. Under this extended anisotropic scattering law, a
Gaussian distribution function is used to describe the variation of scattering
angles in the particle’s local frame.

Results. As a result, we obtain a series of different shock structures and
evolutions in terms of the standard deviation values of the given Gaussian
scattering angular distributions.

Conclusions. We find that the total energy spectral index increases as the
standard deviation value of the scattering angular distribution increases, but
the subshock’s energy spectral index decreases as the standard deviation value
of the scattering angular distribution increases.

###### Key Words.:

acceleration of particles– shock wave–solar energetic particle

##  1 Introduction

It is well known that the diffusive acceleration model has been popular for
more than five decades since Fermi (1949) first proposed that cosmic rays
could be produced via diffusive processes. Until now, diffusive shock
acceleration (DSA) (i.e. first-order Fermi acceleration) has been extensively
applied to many physical systems, such as shocks in the solar system, in our
Galaxy, and throughout the Universe (Baring 1997; Bell 2004).

It is also well known that the nonlinear interactions in plasma usually
include such things as the turbulence of scattering wave field, cosmic ray
(CR) injection, and “back reaction” by CR pressure. These complex behaviors
have held back comprehensive understanding of the DSA and nonlinear DSA
theory. Therefore, to study the properties of the acceleration process and
dynamical behavior of the CR’s “back reaction” on the background flow,
choosing numerical simulation methods has been a primary and essential
problem. (Ostrowski 1991; Kang & Jones 1991; Malkov 1997; Ellison, Baring &
Jones 1995, 1996; Knerr, Jokipii & Ellison 1996; Berezhko et al. 1994;
Shalchi, Li & Zank  2010; Berezhko & Völk 2000). The main simulation methods
are introduced here, and a more detailed review can be seen in Kang (2001).

Monte Carlo method. In Monte Carlo simulations, the scattering processes
between individual particles with the collective background flow are simulated
around a one-dimensional parallel shock. The particle scattering process is
based on a prescribed scattering law, and collecting moments based on the
background computational grid for scattering calculation is done by particle-
in-cell (PIC) techniques. Knerr, Jokipii & Ellison (1996) successfully
developed the dynamical Monte Carlo simulations for the Earth’s bow shock with
important results in the maximum energetic particles achieving greater than
1MeV accelerated by the shock. Before the dynamical Monte Carlo simulation,
Ellison, Möbius & Paschmann (1990) had developed stationary Monte Carlo
simulations with the result that the cutoff energy accelerated by the shock
only reached 100 keV owing to the limitation on the size of bow shock. Baring,
Ellison & Jones  (1995) also used the stationary Monte Carlo method for
simulating the oblique interplanetary shocks, and the calculated results are a
good fit to the observed data. There are many works that use the Monte Carlo
method: Ellison, Baring & Jones (1996), Ellison & Double (2002), Vladimirov,
Ellison & Bykov (2006, 2008), and others.

Hybrid simulation. The particles’ motion equations are solved explicitly based
on the electromagnetic field of the background plasma. Since the proton mass
is about 2000 times that of the electron mass, the total plasma is assumed to
be a coupling of two components: one component (e.g. electrons) is treated as
a massless fluid and the other component(e.g. protons,ions) is treated as
individual particles (Leroy et al. 1982). This method also employs the PIC
techniques based on computational grids. However, the limited computational
resources imply that the extensive calculation of the electromagnetic field
uses unrealistic parameters and are unable to follow the shock for long enough
(Giacalone et al. 1993; Giacalone & Jokipii 2009).

Two-fluid model. The two-fluid model uses the diffusive-convection equation,
coupled with the gas dynamic equations, to simulate the CR’s acceleration as a
gas component and an accelerated particle component (Drury & Falle 1986; Dorfi
1990; Jones & Kang 1992). Since the CR energy density is solved instead of the
particle’s distribution function in this model, the simulation results are
based on some assumptions, such as the CR’s adiabatic index, averaged
diffusive coefficient, and injection rate. The averaged diffusive coefficient
needs to be inferred from the diffusive-advection equation. The acceleration
efficiency dependens on the assumption of the injection rate. Under the
reasonable a priori models of these free parameters, a lot of simulations have
tested the acceleration efficiency and the shock structure, and both agree
with those derived from the diffusion-advection method. And these semi-
analytic solutions that have been extensively used can be seen in recent
works: Malkov, Diamond & Völk (2000), Bednarz & Ostrowski (2001), Blasi (2002,
2004), Amato & Blasi (2005, 2006), Blasi, Amato & Caprioli (2007), Caprioli et
al. (2009a), Caprioli, Amato & Blasi (2010a).

Kinetic simulation. Within this full numerical simulation, the diffusion-
convection equation for the distribution function is solved with a momentum-
dependent diffusion coefficient and a suitable assumption of injection rate
(Kang & Jones 1991; Berezhko et al. 1994; Berezhko & Völk 2000). Unlike the
two-fluid model, the kinetic model should not assume the CR’s adiabatic index,
in addition to using the momentum-dependent diffusion coefficient instead of
the averaged diffusion coefficient. Berezhko and collaborators have studied
numerous DSA models for supernova remnants (SNRs) with the kinetic model and
conclude that about 20 % of the total SN energy transferred to the CR’s
population in the SNRs are more than the results calculated from the two-fluid
model. As for the energy spectrum, the CR spectrum in this model shows a basic
power-law spectrum in the total energy range with a concave curve at some
energy range because of the precursor structure. More detailed new studies of
the kinetic model can be referred to in these papers (Kang, Jones & Gieseler
2002; Kang & Jones 2007; Amenomori et al. 2008; Li et al. 2009; Zirakashvili &
Aharonian 2010).

In an effort to follow and extend the previous dynamical Monte Carlo
simulation (Knerr, Jokipii & Ellison 1996), we independently developed a
simulation code based on the Matlab platform using multiple scattering laws.
Our multiple scattering angular distributions consist of three Gaussian
distributions and one isotropic distribution for the scattering angles during
the scattering process. The aim of the isotropic scattering angular
distribution is to check the dynamical Monte Carlo method independently.
Besides this, we want to know how the Gaussian distributions affect the
scattering angular distribution function and the shock wave’s evolution and
propagation; even more, we expect to find the relationships between the
multiple scattering law and the shock compression ratio. To validate the
multiple scattering angular distributions, we followed the parallel-plane
collisionless shock and the particle’s acceleration using the same parameters
and data as from Earth’s bow shock, which was used in previous dynamical Monte
Carlo simulation.

In Section 2, we outline the motivations for performing these four cases with
four different scattering angular distributions. In Section 3, the specific
simulation techniques are described. In Section 4, we present the shock
simulation results and different cases with four assumptions for scattering
angle distributions. Section 5 includes a summary and the conclusions.

##  2 The model

The dynamical Monte Carlo simulation has been developed by Knerr, Jokipii &
Ellison (1996) to study Earth’s bow shocks. It gives good results for the
higher than 1MeV cutoff in energy particles and the power-law energy tail in
the energy spectra. The dynamical Monte Carlo simulation method uses the
prescribed scattering law instead of the complex electromagnetic field
calculation like in the hybrid model. In addition, the dynamical Monte Carlo
simulation need not assume the CR’s injection rate and the associated
diffusive coefficient as do the two-fluid and kinetic models. For the above
reasons, we consider that developing a simulation code by following the
previous dynamical simulation is necessary. Although the previous results
successfully agree with observed data, the authors mention that their results
show that the total compression ratio of the shock is more than 4, which
should be less than the ratio of the standard value for a nonrelativistic
shock (Pelletier 2001). The Rankine-Hugoniot (RH) jump conditions allow
deriving the relation of the compression ratio with the Mach number:
$r=(\gamma_{a}+1)/(\gamma_{a}-1+2/M^{2})$, for a nonrelativistic shock, the
adiabatic index $\gamma_{a}$ = 5/3 , if the Mach number $M\gg 1$, then the
maximum compression ratio should be 4. To validate these consistent results
from the previous model and extend this study to find what might be
responsible for the shock compression ratio, we extend the previous isotropic
scattering angular law by including an anisotropic scattering angular law.
This prescribed multiple scattering law consists of an isotropic scattering
angular distribution and an anisotropic scattering angular distribution. The
scattering angles consist of two variables of pitch angle and azimuthal angle.
Once a particle has a collision with the massive scattering centers, its pitch
angle becomes $\theta^{\prime}$=$\theta$+$\delta\theta$, and the azimuthal
angle becomes $\phi^{\prime}$=$\phi$+$\delta\phi$, where $\delta\theta$ is the
variation in the pitch angle $\theta$, and $\delta\phi$ is the variation in
the azimuthal angle $\phi$. The pitch angles $\theta$ and $\theta^{\prime}$
are both in the range $0\leq\theta,\theta^{\prime}\leq\pi$, and azimuthal
angles $\phi$ and $\phi^{\prime}$ are both in the range
$0\leq\phi,\phi^{\prime}\leq 2\pi$ on the unit sphere. The variation in the
pitch angle $\delta\theta$ and azimuthal angle $\delta\phi$ are composed of
the scattering angle, and its anisotropic character is described by the
Gaussian function $f(\delta\theta,\delta\phi)$.

Under the multiple scattering angular distribution law, four cases are
calculated with three Gaussian distributions and one isotropic random
distribution for the scattering angles. Here, the sign $\sigma$ is used to
represent the standard deviation of the Gaussian function, and the sign $\mu$
is used to represent the statistical average or expected value of Gaussian
function for the scattering angles ($\delta\theta,\delta\phi$). We catalog the
four cases as follows.

(1) Case A: the scattering angles ($\delta\theta,\delta\phi$) are distributed
with a standard deviation $\sigma=\pi$/4 and an average value $\mu=0$.

(2) Case B: the scattering angles ($\delta\theta,\delta\phi$) are distributed
with a standard deviation $\sigma=\pi$/2 and an average value $\mu=0$.

(3) Case C: the scattering angles ($\delta\theta,\delta\phi$) are distributed
with a standard deviation $\sigma=\pi$ and an average value $\mu=0$.

(4) Case D: the scattering angles ($\delta\theta,\delta\phi$) are distributed
with a standard deviation $\sigma=\infty$ and an average value $\mu=0$, with
$\delta\theta$ varying from $-\pi/2$ to $\pi/2$, and $\delta\phi$ varying from
$-\pi$ to $\pi$ isotropically.

We performed four simulations according to the four different assumptions of
the scattering angular distributions algorithm. We also assumed the scattering
time (i.e., the mean time between two scattering events) is the same constant
in the four cases as in the previous model. The idea that such a simple law
can be used to describe the entire scattering process was postulated by
Eichler (1979), based on the two-stream instability work done by Parker
(1961). Put simply, it is assumed that the turbulence generated by both
energetic particles streaming in front of the shock and by thermal particles
produces nearly elastic scattering for particles for all energies in diffusive
shocks.

##  3 Description of the method

For simulating the total properties of shocks as they evolve from formation to
a final steady state as energy increases via Fermi acceleration, we used the
dynamical Monte Carlo model which employed the PIC techniques. Because there
is no assumption of the injection rate or transparency function in PIC
techniques, the shock-heated downstream ions can freely scatter back across
the subshock into the upstream region without being thermalized, and the
superthermal particles are produced in the thermal background self-
consistently. In addition, unlike the hybrid simulation, there is no
complicated electromagnetic field calculation for individual particles,
because it is replaced by the prescribed scattering law (Ellison et.al. 1993).
To reproduce these acceleration and scattering processes, a similar simulation
box and the same parameters (see Appendix A) as the previous dynamical Monte
Carlo method are used in these new codes. As described in the previous
simulation (see Figure 8), the particles with an initial bulk velocity $U_{0}$
and a Maxwellian thermal velocity $V_{L}$ in their local frame are moving
along a parallel magnetic field $B_{0}$ in a one-dimensional box. To maintain
a continuous flux-weighted flow upstream, a new particle fluid with the same
density needs to be injected into the simulation box from the left boundary.
For a shock initialization, the reflecting wall on the right boundary is used
to reflect the incoming particles, and forms a piston shock. The model also
includes the escape of energetic particles at the upstream free escape
boundary (FEB). The FEB phenomenologically models a finite shock size or the
lack of sufficient scattering far upstream to turn particles around. Once ions
cross the FEB, they are assumed to decouple from the shock system, and are
taken as the energy losses (Jones & Ellison 1991). The size of the foreshock
region (the distance from the shock front to the FEB) thus sets a limit on the
maximum energy a particle can obtain.

As shown in Figure 8, one particle’s box frame velocity $V$ is a total
velocity, which is composed of the local thermal velocity $V_{L}$ and the bulk
fluid speed $U$ (i.e. $V=V_{L}+U$, for upstream $U=U_{0}$, for downstream
$U=0$). After one particle arrives in the downstream region, its kinetic
energy is converted into random thermal energy by dissipation processes. With
the development of these many processes, the bulk fluid speed of downstream
flow becomes zero, and the length of downstream region is extended
dynamically.

As listed in Table 2, all of the specific parameters are used in our
simulations, considering PIC techniques. The total length ($X_{max}=300$) is
divided into the number of grids ($n_{x}=600$) with a grid length ($\Delta
x=1/2$). Initial grid density of the particles ($n_{0}=650$) is set. The total
time ($T_{max}=2400$) is divided into the number of time steps ($N_{t}=72000$)
with an increment of time ($dt=1/30$). In summary, these new codes consist of
the following three substeps like the previous simulation, except for the
third substep employing the extended multiple scattering laws.

(i) Individual particles move. Particles with their velocities move along the
one-dimensional $x$ axis:

$X^{t}_{k}=X^{t-1}_{k}+(V_{x})^{t}_{k}dt,t\in[1,t_{max}],k\in[1,k_{max}],$ (1)

where

$(V_{x})_{k}=(V_{Lx})_{k}+(U_{k})_{x},$ (2)

$(U_{k})_{x}=\frac{1}{n_{k}}\displaystyle\sum_{i=1}^{n_{k}}(V_{x})_{i}.$ (3)

here $t_{max}$=72000, $k_{max}$=600, and “k” represents the index of the
computational grid, $(U_{k})_{x}$ represents the bulk fluid speed of the
computational grid along to the $\hat{x}$ direction, and the value of $U_{k}$
is obtained from substep (ii). Since we are simulating a diffusive shock based
on a one-dimensional parallel magnetic field, the fluid quantities only vary
in the $\hat{x}$ direction.

(ii) Mass collection. Summation of particle masses and velocities are
calculated at the center of each computational grid:

$P_{k}=\displaystyle\sum_{i=1}^{n_{k}}m_{p}(V_{x})_{i},k=1,2,...k_{max}$ (4)
$U_{k}=\frac{1}{n_{k}}\displaystyle\sum_{i=1}^{n_{k}}(V)_{i},k=1,2,...k_{max},$
(5)

where $n_{k}$ is the number density of particles in the “k” grid, representing
the mass of the computational grid. Here, $P_{k}$ is the total momentum of the
protons in the “k” grid, $m_{p}$ is the mass of an individual proton, and
$U_{k}$ is the average bulk fluid speed of the grid (i.e. the velocity of the
scattering center). The collected grid-based mass and momentum densities will
directly decide the velocity of the scattering center $U_{k}$. The particle’s
total velocity $V$ in the box frame is decided by Equation (2). Once the value
of $U_{k}$ becomes zero, the shock front is decided by the position of the
corresponding grid, and it means that the shock is formed and the length of
the downstream region is extended dynamically. Similarly, if the value of
$U_{k}$ is between $U_{0}$ and zero, it means that the foreshock region or
precursor (i.e. between the FEB and the shock front) is formed by the “back
pressure.” The FEB and the shock front both dynamically move away from the
reflective wall with a shock velocity $v_{sh}$.

(iii) Applying multiple scattering laws. A certain fraction of the particles
are chosen to scatter the background scattering center with their
corresponding scattering angles according to the prescribed scattering angular
distributions. The average number of scattering events occurring in an
increment of time $dt$ depends on the scattering time scale $\tau$, and the
scattering rate is presented by

$R_{s}=dt/\tau,$ (6)

where $R_{s}$ is the probability of the scattering events occurring in an
increment of time. The candidates with their local velocities and scattering
angles scatter off the grid-based scattering centers. These individual
particles do not change their routes until they are selected to scatter once
again. So the particle’s mean free path is proportional to the local thermal
velocities in the local frame with

$\lambda\propto V_{L},$ (7)

for simplicity, we take its formula as

$\lambda=V_{L}\cdot\tau.$ (8)

For the individual protons, the grid-based scattering center can be seen as a
sum of individual momenta. So these scattering processes can be taken as the
elastic collisions. In an increment of time, once all of the candidates
complete these elastic collisions, the momentum of the grid-based scattering
center is changed. In turn, the momentum of the grid-based scattering center
will affect the momenta of the individual particles in their corresponding
grid in the next increment time. One complete time step consists of the above
three substeps. The total simulation temporally evolves forward by repeating
this time step sequence. To calculate the scattering processes accurately and
produce an exponential mean free path distribution, the time step should be
less than the scattering time (i.e. $dt<\tau$).

The presented multiple scattering law simulations are developed on the Matlab
platform. Any one of the four cases can occupy the CPU time for about seven
weeks on a 3. 4GHZ (MF) CPU per core. To speed up the running programs, the
parallel algorithm should be used on a high performance computer (HPC).

##  4 Results & discussions

We present all of the shock profiles for the shock simulations of the four
cases in Figure 1, and we present all aspects of simulation results including
the density and velocity profiles, compression ratios, analysis of the heating
and acceleration, and energy spectrum, as well as the correlations between the
shock compressions with the energy spectral index. For the convenience of
comparison and discussion, we list the specific calculated items in Table 1.
Here, $\sigma$ is the standard deviation of the scattering angular
distribution function $f(\delta\theta,\delta\phi)$, and $\mu$ is the average
value of scattering angles $(\delta\theta,\delta\phi)$. The subshock
compression ratios $r_{sub}$ are calculated from the velocity profiles in the
same shock frame reference. The compression ratios $r_{u}$ and $r_{\rho}$ are
calculated from velocity profiles and density profiles, respectively. The
total energy spectral index $\Gamma_{tot}$ and the subshock’s energy spectral
index $\Gamma_{sub}$ are calculated from compression ratios $r_{u}$ and
$r_{sub}$, respectively. The last two rows are shown as scaled values.

Table 1:  Results of calculation with an initial energy of $E_{0}$=1.3105keV. Items | Case A | Case B | Case C | Case D
---|---|---|---|---
$f(\delta\theta,\delta\phi)$ | Gaussian | Gaussian | Gaussian | Isotropy
$\sigma$ | $\pi/4$ | $\pi/2$ | $\pi$ | $\infty$
$\mu$ | 0 | 0 | 0 | 0
$v_{sub}$ | 0\. 1118 | 0\. 1465 | 0\. 1733 | 0\. 2159
$v_{sh}$ | -0. 0433 | -0. 0535 | -0. 0617 | -0. 0733
$X_{sh}$ | 196 | 171\. 5 | 152 | 124
$FEB$ | 106 | 81\. 5 | 62 | 34
$\rho_{2}$ | 5000 | 4225 | 3759 | 3303
$\rho_{1}$ | 650 | 650 | 650 | 650
$r_{sub}$ | 2\. 5421 | 3\. 0207 | 3\. 3903 | 3\. 9444
$r_{u}$ | 7\. 9231 | 6\. 6031 | 5\. 8649 | 5\. 0909
$r_{\rho}$ | 7\. 6936 | 6\. 5001 | 5\. 7836 | 5\. 0815
$\Gamma_{tot}$ | 0\. 7167 | 0\. 7677 | 0\. 8083 | 0\. 8667
$\Gamma_{sub}$ | 1\. 4727 | 1\. 2423 | 1\. 1275 | 1\. 0094
$V_{Lmax}$ | 10\. 8251 | 16\. 0001 | 17\. 7824 | 20\. 5286
$E_{cutoff}$ | 1\. 10MeV | 2\. 41MeV | 2\. 98MeV | 4\. 01MeV
$|V_{sh}|$ | 58\. 16km/s | 71\. 87km/s | 82\. 88km/s | 98\. 46km/s

As shown in Figure 1, the present isotropic Case D largely appears similar to
the results from previous dynamical simulations by Knerr, Jokipii & Ellison
(1996). In addition, all aspects of the shock wave structure, density and
velocity profiles, compression ratio and energy spectra present in isotropic
Case D also give similar results to the previous outcome. The specific results
of the present isotropic Case D are shown in the fifth column of Table 1:
$v_{sh}$=-0.0733, $X_{sh}$=124, $FEB$=34, $r_{sub}$=3.9444, $r_{u}$=5.0909,
$r_{\rho}\sim r_{u}$ with a difference of 0.18%, $\Gamma_{sub}$=1.0094,
$\Gamma_{tot}$=0.8667, $|V_{sh}|$=98.46km/s, and $E_{cutoff}$=4.01MeV. As a
comparison, the corresponding results of the previous simulation are also
given: $v_{sh}$=-0.0720, $X_{sh}$=127.5, $FEB$=37.5, $r_{sub}$=3.20,
$r_{u}$=5.20, $r_{\rho}\sim r_{u}$ with a difference of 2.3% ,
$\Gamma_{sub}$=1.1818, $\Gamma_{tot}$=0.8571, $|V_{sh}|$=96.9 km/s, and
$E_{cutoff}\sim$4.00MeV. Among comparable results, a slightly larger
differences in the values of the $r_{sub}$ and $\Gamma_{sub}$ between the two
isotropic simulations would be distributed between different sizes of the
subshock region, decided differently. As seen from the comparison of the
results coming from two independent simulation codes, the present simulation
code successfully produced good agreement in the results with those in
previous dynamical Monte Carlo simulation. Therefore, the present simulation
code is based on the Matlab platform without using a supercomputer that can
independently validate the previous dynamical simulation method using a
completely different code for that supercomputer. Next, we offer a series of
discussions about the different cases considering the specific aspects of the
simulation results for diffusive shock.

Figure 1: Four cases of density profiles for the entire simulation box vs the
time. The dashed line represents the position of the FEB in each plot.

Figure 2: The velocity and density profiles (in the box frame) vs the position
at the end of the simulation ($T_{max}$=2400) in the four cases. The vertical
dashed line in each plot represents the position of the FEB, from case A to
case D.

###  4.1 Shock Profiles

Figure 1 shows time sequences of the density profiles of four cases. In each
plot, a shock forms and moves away from the reflective wall, and the dashed
line represents the FEB position with the time in each case parallel to the
shock front position. We can see that both the shock position and the FEB
position are moving with a virtually constant velocity from the beginning of
the simulation to the end of the simulation (i.e. $T_{max}$=2400) in each
case. Simultaneously, as far as the positions of the FEB are concerned, we can
see that the FEB position at the end of the simulation is significantly
different in four different cases. As for the average density fluctuation in
the downstream region, there are also apparent changes in different cases, and
case A has the slowest shock propagation speed among these four cases. Case D
has the lowest average density profile in the downstream region among these
four cases. Because from Cases A to D the only difference is the prescribed
scattering angular distribution, we conclude that these differences of the
results for shock propagation speed and density profiles are decided by the
standard deviation value $\sigma$ of the scattering angular distribution.

Figure 2 shows four cases of density and velocity profiles at the end of the
simulations. From Cases A to D, the position of the FEB approaches zero as the
value of the standard deviation $\sigma$ increases. The effects of the
accelerated particles are clearly seen in the upstream smoothing of the
velocity profiles in each case. In the simulations, when high-energy particles
cross the shock front and diffuse upstream, they contribute negatively to the
velocity profile. This reduces the grid-based velocity in the zones upstream
of the shock, which in turn affects particles that are scattering in that
region. In fact, the accelerated particles slow and heat the incoming flow and
smooth the shock transition by their “back reaction.” As is obvious to see
from the velocity and density plots, different scattering angular
distributions produce different effects on the shock wave evolution. For the
examples presented here, we consider that a difference of approximately 40.93%
of the shock velocity is contributed by the scattering angle distribution.

###  4.2 Compression ratios

Here, we compare the compression ratios calculated from the velocity profiles
with those from the density relationships. First, the value of the total
compression ratio can simply be calculated from the formula

$r_{u}=u_{1}/u_{2}$ (9)

where $u_{1}=u_{0}+|v_{sh}|$ , $u_{2}=|v_{sh}|$, and $u_{1}(u_{2})$ is the
upstream (downstream) velocity in the shock frame. The shock velocity at the
end of the simulation ($T_{max}$) can be derived from the formula

$v_{sh}=(X_{max}-X_{sh})/T_{max}$ (10)

where $X_{max}$=300, $T_{max}$=2400, and $X_{sh}$ is the position of the shock
at the end of the simulation (see Table 1). The specific calculated results
are shown in Table 1.

But in terms of the specific shock structure as seen in Figure 3, an accurate
subshock compression ratio calculation should be more complicated. In any one
of the cases in Figure 3 (plotted in the box reference frame), we show the
specific aspects of a shock modified by an energetic particle population whose
mean-free-path is an increasing function of momentum. The shock structure in
each plot consists of three main parts: precursor, subshock, and downstream.
The smooth precursor is on the largest length scale between the FEB and near
shock position $X_{sh}$, where the fluid speed gradually decreases from value
$U_{0}$ to $v_{sub}$. The size of this precursor is almost the mean-free-path
length of the maximum energy accelerated particles. One of the smallest scales
is the subshock region with a sharp deflection of the fluid speed decreasing
from $v_{sub}$ to $v_{box}=0$. The downstream region changes after the fluid
speed becomes $v_{box}=0$ by microphysical dissipation processes. The gas
subshock is just an ordinary discontinuous classical shock embedded in the
comparably larger scale energetic particle shock (Berezhko & Ellison  1999).
The value of $v_{sub}$ is determined by a sharp deflection of smooth curves in
velocity profiles near the shock front, and the value of the subshock velocity
increases from cases A, B, and C to Case D (i.e.
$(v_{sub})_{A}<(v_{sub})_{B}<(v_{sub})_{C}<(v_{sub})_{D}$). All of the
velocity profiles are based on the box frame. That value of the box frame’s
velocity is zero ($v_{box}=0$) in all cases. The subshock compression ratio
$r_{sub}$ is calculated from the formula $r_{sub}=(v_{sub}+|v_{sh}|)/|v_{sh}|$
. For the sake of the comparison of the values of $r_{sub}$ in different
cases, the subshock compression ratios are calculated in the same shock frame
reference, and the calculated results of $r_{sub}$ are shown in Table 1.

Figure 3: Velocity profiles in the shock region at the end of the simulation
($T_{max}$=2400) in the four cases; the vertical solid and dotted lines
indicate the shock front and FEB in each plot, respectively; the horizontal
solid, dotted, dot-dashed and dashed lines show the values of the shock
velocity $v_{sh}$, velocity of box frame $v_{box}$, subshock velocity
$v_{sub}$ and initial bulk velocity $U_{0}$, respectively. Two vertical bars
in each plot represent the two deflections of velocity, the upper bar
represents the part of the shock precursor and the lower one represents the
subshock. Figure 4: Density profiles in an entire simulation box at the end of
the simulation ($T_{max}$=2400) in the four cases. The vertical solid line is
located in the position of the shock front, the upper horizontal dot-dashed
line represents the value of the downstream density $\rho_{2}$, and the lower
horizontal dashed line indicates the value of the upstream density $\rho_{1}$
in each plot.

We then have the calculations of the compression ratio from the density
relationships between the upstream and downstream flows:

$r_{\rho}=\rho_{2}/\rho_{1}$ (11)

where $\rho_{1}=n_{0}$ is the upstream density, and $\rho_{2}$ is the
downstream density. This value is presented by which is the average value of
$n_{k}$ over the downstream region.

$\rho_{2}=\frac{1}{(k_{max}-k_{sh})}\displaystyle\sum_{k=k_{sh}}^{k_{max}}(n_{k})$
(12)

where $n_{k}$ is the number density of particles in the “k” grid,
$k_{sh}=x_{sh}/dx$ is the grid index of the shock at the end of the simulation
($T_{max}=2400$), and $k_{max}=600$ is the grid index of the $X_{max}$.

Figure 4 shows the complete density plots of the four cases at the end of the
simulation. The value of the upstream density $\rho_{1}$ is the same constant
value, which is equal to the initial density $n_{0}$ in each case. The value
of the downstream density $\rho_{2}$ decreases from cases A, B, and C to Case
D (i.e. $(\rho_{2})_{A}>(\rho_{2})_{B}>(\rho_{2})_{C}>(\rho_{2})_{D}$).
Similarly, the detailed calculation results of the compression ratios
$r_{\rho}$ are listed in Table 1. As listed in Table 1, the values of the
subshock compression ratios, $r_{sub}$=2.5421, $r_{sub}$=3.0207,
$r_{sub}$=3.3903, and $r_{sub}$=3.9444, corresponding to cases A, B, C, and D,
respectively, are all lower than the standard value of $r=4$. Unfortunately,
the values of total compression ratio $r_{u}$ and $r_{\rho}$ in each case are
both higher than the standard value of $r=4$. But Knerr, Jokipii & Ellison
(1996) consider that, if energy is lost from the system (e.g. by particles
escaping via FEB), it is possible to produce a shock with a total compression
ratio that is higher than the standard value predicted by the Rankine-Hugoniot
(RH) conditions. We have examined the mass and energy losses via the FEB in
each case. The results definitely show that the case with more energy losses
would produce a higher total compression ratio than those in the case with
less energy loss. Consequently, we consider that the energy loss rates would
be affected by the prescribed scattering law. In any case, the energy losses
are always an important and interesting problem in the nonlinear diffusive
shock acceleration theory, so we will perform more precise research focusing
on these problems in later papers. In addition, although the values of $r_{u}$
are correspondingly slightly higher than the value of $r_{\rho}$ in each case,
all these differences are less than 3%, and the specific difference in each
case is 2.9%,1.5%,1.3% and 0.18% corresponding to cases A, B, C and D. As seen
from Figs. 3 and 4, the value of the total compression ratio $r_{tot}$,
determined from the velocity profiles, is more consistent with the density
profiles in each case (i.e. the total compression ratios $r_{tot}$ in all
cases are satisfied by the Rankine-Hugoniot (RH) conditions :
$u_{1}/u_{2}=\rho_{2}/\rho_{1}$). Therefore, it is not difficult for us to
conclude that the total compression ratios $r_{u}$ and $r_{\rho}$ decrease as
the value of the standard deviation $\sigma$ of the scattering angular
distribution increases, but the subshock compression ratio $r_{sub}$ increases
as the value of the standard deviation $\sigma$ of the scattering angular
distribution increases.

###  4.3 Heating & acceleration

Figure 5: The scatter plots of the particle’s thermal velocities in the local
frame vs its position at the end of the simulations ($T_{max}$=2400), and the
vertical dashed line and solid line in each plot indicate the approximate
position of the FEB and shock front, respectively. Only the ratio of 1/12 of
the total number of particles are plotted.

Here we contrast between two important aspects of the heating and acceleration
processes in the diffusive shock acceleration. Figure 5 shows the particles’
scatter plots in four cases at the end of the simulations ($T_{max}$=2400). In
each case, particles with local velocity scatter into the simulation box’s
position. A large amount of particles that do not get injected into the Fermi
acceleration mechanism and that have lower thermal velocities stay in the
downstream region, and a few particles with higher energy via multiple shock
encounters can move far away from the shock front and even escape from the
FEB. From Cases A to D, more and more particles are injected into the Fermi
acceleration mechanism, and they gain greater and greater maximum energy.
Obviously, the maximum thermal velocity in Case D would be several times that
of the ones in Case A. We supposed that this difference is mainly contributed
by the scattering angular distribution function $f(\delta\theta,\delta\phi)$.
In short, the majority of particles, which flow toward the shock, cross the
shock only once, and (after scattering) remain fairly stationary in the
downstream region, which would consist of the “heated” elements, and a few
high-energy particles represent the “power-law” part of the simulated
particles flows. Actually, the “back pressure” from the accelerated particles
via their back reaction reduces the incoming fluid speed, leading to a
smoothed precursor heating. Therefore, an anisotropic scattering angular
distribution in the shocks dominate the “gas heating” process in the simulated
plasma, and isotropic scattering angular distribution in the shocks play an
important role in nonlinear acceleration of the energetic particles via the
Fermi mechanism, and they dominate the “precursor heating”.

###  4.4 Energy spectra

Figure 6: The final energy spectrum in the four cases and the initial energy
spectrum plots. The thick solid line with a narrow peak at $E_{0}=$1.3105keV
represents the initial Maxwell energy distribution. The solid, dashed, dot-
dashed, and dotted extended curves indicate the simulated particles’ energy
spectral distribution, averaged over the entire downstream region, at the end
of the simulations ($T_{max}$=2400), corresponding to Cases A, B, C, and D,
respectively. Most particles cross the shock only once, producing the large
broad peak centered at E${}_{A}\sim$ 0.05keV, E${}_{B}\sim$0\. 1keV,
E${}_{C}\sim$ 0.15keV, and E${}_{D}\sim$0.20keV in Cases A, B, C, and D,
respectively. However, some particles gain enough energy via the Fermi
acceleration mechanism to produce the “power-law” tail in the energy spectrum
with the cutoff at $E_{A}$=1.10 MeV, $E_{B}$=2.41 MeV, $E_{C}$=2.98 MeV, and
$E_{D}$=4.01 MeV corresponding to Cases A, B, C, and D. All these energy
spectra are plotted in the same shock frame.

Spectra calculated in the shock frame from the initial and final particle
(proton) energy distributions in all cases are shown in Figure 6. The energy
units in this plot are derived from the scaling parameters presented in Table
2 (see Appendix A). Initially, all particles move toward the wall with a
certain thermal spread in energy. A narrow peak at E=1.3105keV represents the
initial Maxwell energy distribution. The four extended curves indicate the
simulated particle energy spectral distribution, averaged over the entire
downstream region, at the end of the simulations, corresponding to the four
cases, respectively. The majority of the particles cross the shock only once,
producing an expanded energy spectrum with a central peak at E${}_{A}\sim$
0.05keV, E${}_{B}\sim$0.1keV, E${}_{C}\sim$ 0.15keV, and E${}_{D}\sim$0.20keV
in Cases A, B, C, and D, respectively. However, as is shown in Figure 6, the
minority of the particles gain enough energy via the Fermi acceleration
mechanism to produce the “power-law” tail in the energy spectrum with the
cutoff at $E_{A}$=1.10 MeV, $E_{B}$=2.41 MeV, $E_{C}$=2.98 MeV and
$E_{D}$=4.01 MeV corresponding to Cases A, B, C and D, respectively. For more
details about the calculated results, see Table 1. It is evident from Figure 6
that the values of the central peak of the extended energy spectra in the four
cases are far from the initial energy peak in their respective order. This
means the values of the central peak in each case increase as the value of the
standard deviation of the scattering angular distribution increases, and each
extended curve shows a harder power-law slope in its high-energy tail as the
expand energy range increases, respectively. Therefore, we can see that the
case of applying an anisotropic scattering angular distribution function will
produce a slightly softer energy spectrum, and the case of applying an
isotropic scattering angular distribution will produce a slightly harder
energy spectrum.

###  4.5 Spectral index & compression ratios

Usually, we could predict the power-law energy spectral index from diffusive
shock acceleration theory:

$dJ/dE\propto E^{-\Gamma}$ (13)

where $dJ/dE$ is the energy flux and $\Gamma$ is the energy spectral index,
and

$\Gamma=(r+2)/(2\times(r-1)).$ (14)

According to Equation 14, we substituted the values of the compression ratio
$r$ with two group values of $r_{tot}$ and $r_{sub}$ obtained in each case.
Then, the two group energy spectral indices $\Gamma_{tot}$ and $\Gamma_{sub}$
in each case are calculated. Two groups’ spectral index values are listed in
Table 1 as $\Gamma_{A}$= 0.7167, $\Gamma_{B}$ =0.7677, $\Gamma_{C}$ =0.8083,
and $\Gamma_{D}$=0.8667 in the total group $\Gamma_{tot}$, and $\Gamma_{A}$=
1.4727, $\Gamma_{B}$ =1.2423, $\Gamma_{C}$ =1.1275, and $\Gamma_{D}$=1.0094 in
subshock group $\Gamma_{sub}$, corresponding to the cases A, B, C, and D,
respectively. As shown in Figure 7, from Cases A to D, all of the values of
the subshock’s energy spectral index $\Gamma_{sub}>1$ and show that a slightly
harder power-law slope in the respective order. From Cases A to D, all of the
values of the total energy spectral index $\Gamma_{tot}<1$ show a slightly
decreasingly deviation from the power-law slope in the respective energy
spectrum.

Figure 7: The correlation of the deviation value of the Gaussian distribution
vs the energy spectral index. The triangles represent the total energy
spectral index in each case. The circles indicate the subshock’s energy
spectral index in each case. From Cases A to D, all of the values of the
subshock’s energy spectral index $\Gamma_{sub}>1$ show a slightly harder
power-law slope in the respective order. From Cases A to D, all of the values
of the total energy spectral index $\Gamma_{tot}<1$ show a slightly decreasing
deviation from the power-law slope in the respective energy spectrum.

##  5 Summary and conclusions

We followed the previous dynamical Monte Carlo simulation by using a new code
based on the Matlab platform independently, and presented the same results as
the outcome from the previous simulation. In addition, we successfully
extended the simulation include the multiple scattering angular distributions
using these new codes to study the diffusive shock acceleration mechanism
further.

In conclusion, the comparison of the calculated results come from different
extensive cases, we find that the total energy spectral index increases as the
standard deviation value of the scattering angular distribution increases, but
the subshock’s energy spectral index decreases as the standard deviation value
of the scattering angular distribution increases. In these multiple scattering
angular distribution simulation cases, the prescribed scattering law dominates
the shock structure and plays an important role in balancing whether the
particles have more heating or more acceleration. In other words, the cases
with anisotropic scattering distribution give the overall velocity-deflection
precursor sizes, which are larger than the isotropic case, and give a
relatively greater “heating” effect or less “acceleration” effect on
background flows than does the isotropic case.

As a result, the shock compression ratio and the energy spectral index are
both modified naturally by the prescribed scattering law. Specifically, the
cases of applying an anisotropic scattering distribution function will produce
a slightly softer subshock’s energy spectral index, and the case of applying
an isotropic scattering angular distribution will produce a slightly harder
subshock’s energy spectral index. Simultaneously, from the isotropic case to
the anisotropic case, the total energy spectrum shows an increasing deviation
from the “power-law” distribution.

In addition, although we find no case producing the total compression ratio
which should be less than the standard value 4 according to the Rankine-
Hugoniot (RH) jump conditions, the fact is clear that the prescribed
scattering angular distribution function would have an effect on the total
compression ratio. If there is a suitably prescribed scattering law that leads
to much less energy loss, it is possible to constrain the total compression
ratio to be less than 4.

###### Acknowledgements.

The authors would like to thank Doctors G. Li, Hongbo Hu, Siming Liu, Xueshang
Feng, and Gang Qin for many useful and interesting discussions concerning this
work. In addition, we also appreciate Profs. Qijun Fu and Shujuan Wang, as
well as other members of the solar radio group at NAOC.

## References

  * Amenomori et al. (2008) Amenomori, M.; Bi, X. J.; Chen, D. & coworkers, 2008 Astrophys. J., 678, 1165\.
  * Amato & Blasi (2005) Amato, E. & Blasi, P., 2005, M.N.R.A.S.Lett., 364, 76
  * Amato & Blasi (2006) Amato, E. & Blasi, P., 2006, M.N.R.A.S., 371, 1251
  * Bednarz & Ostrowski (2001) Bednarz, J. & Ostrowski, M., 2001, M.N.R.A.S., 310, L13.
  * Baring, Ellison & Jones (1995) Baring, M. G., Ellison, D. C. & Jones, F. C., 1995, Adv. Space Res., 15, 397
  * Baring (1997) Baring, M. G., 1997, Very High Energy Phenomena in the Universe; Morion Workshop, Edited by Y. Giraud-Heraud and J. Tran Thanh Van, 1997., p.97
  * Bell (2004) Bell, A.R., 2004, M.N.R.A.S., 353, 550
  * Berezhko et al. (1994) Berezhko E.G., Yelshin V.K. & Ksenofontov L.T. 1994, Astropart. Phys., 2, 215
  * Berezhko & Ellison (1999) Berezhko E.G. & Ellison, D. C. 1999, Astrophys. J., 526, 385
  * Berezhko & Völk (2000) Berezhko E.G. & Völk H.J. 2000, Astron. Astrophys., 357, 283
  * Blasi (2002) Blasi, P., 2002, A.Ph., 16, 429
  * Blasi (2004) Blasi, P., 2004, A.Ph., 21, 45
  * Blasi, Amato & Caprioli (2007) Blasi, P., Amato, E. & Caprioli, D., 2007, M.N.R.A.S., 375, 1471
  * Caprioli et al. (2009a) Caprioli, D., Blasi, P., Amato, E. & Vietri, M., 2009a, M.N.R.A.S., 395, 895
  * Caprioli, Amato & Blasi (2010a) Caprioli, D., Amato, E. & Blasi, P., 2010a, A.Ph., 33, 160
  * Drury & Falle (1986) Drury, L. O’C. & Falle, S. A. E. G. 1986, M.N.R.A.S., 223, 353
  * Dorfi (1990) Dorfi, E.A., 1990, Astron. Astrophys., 234, 419
  * Ellison, Möbius & Paschmann (1990) Ellison, D. C., Möbius, E. & Paschmann, G. 1990, Astrophys. J., 352, 376
  * Ellison et.al. (1993) Ellison, D. C., Giacalone, J. , Burgess, D., & Schwartz, S. J., 1993, J. Geophys. Res., 98, 21058
  * Ellison, Baring & Jones (1995) Ellison, D. C., Baring, M. G. & Jones, F. C., 1995, Astrophys. J., 453, 873
  * Ellison, Baring & Jones (1996) Ellison, D. C., Baring, M. G. & Jones, F. C., 1996, Astrophys. J., 473, 1029
  * Ellison & Double (2002) Ellison, D. C. & Double, G. P., 2002, A.Ph., 18, 213
  * Eichler (1979) Eichler, D., 1979, Astrophys. J., 229, 419
  * Fermi (1949) Fermi, E., 1949, Phys. Rev. Lett.75, 1169
  * Giacalone et al. (1993) Giacalone, J., Burgess, D, Schwartz, S. J. & Ellison, D. C., 1993, Astrophys. J., 402, 550
  * Giacalone & Jokipii (2009) Giacalone, J. & Jokipii, J. R. 2009, Astrophys. J., 701, 1865.
  * Pelletier (2001) Pelletier, G. 2001, Lecture Notes in Physics, , 576, 58
  * Jones & Ellison (1991) Jones, F. C. & Ellison, D. C. 1991, Space Sci. Rev., 58, 259.
  * Jones & Kang (1992) Jones, T. W. & Kang, H. 1992, Astrophys. J., 396, 575
  * Kang & Jones (1991) Kang H. & Jones T.W. 1991, M.N.R.A.S., 249, 439
  * Kang (2001) Kang, H., 2001, Astrophys. J., 550, 737
  * Kang, Jones & Gieseler (2002) Kang, H., Jones, T. W. & Gieseler, U.D.J., 2002, Astrophys. J.579, 337
  * Kang & Jones (2007) Kang, H., Jones, T. W., 2007, A.Ph., 28, 232
  * Knerr, Jokipii & Ellison (1996) Knerr, J. M., Jokipii, J. R. & Ellison, D. C. 1996, Astrophys. J., 458, 641
  * Leroy et al. (1982) Leroy, M. M., Winske, D., Goodrich, C. C., Wu, C. S., & Papadopoulos, K., 1982, J. Geophys. Res., 87, 5081
  * Li et al. (2009) Li, G., Webb, G., Shalchi, A. & Zank, G. P., 2009, 18th Annual International Astrophysics Conference. AIP Conference Proceedings, 1183, 57
  * Malkov (1997) Malkov, M. A., 1997, Astrophys. J., 485, 638
  * Malkov, Diamond & Völk (2000) Malkov, M. A., Diamond P. H. & Völk, H. J., 2000, Astrophys. J.Letters, 533, 171
  * Parker (1961) Parker, E. N., 1961, Nucl. Energy, Vol.C2, 146
  * Ostrowski (1991) Ostrowski, M. , 1991, M.N.R.A.S., 249, 551
  * Shalchi, Li & Zank (2010) Shalchi, A. , Li, G. & Zank, G. P., 2010, Astrophysics and Space Science, 325, 99.
  * Vladimirov, Ellison & Bykov (2006) Vladimirov, A., Ellison, D. C. & Bykov, A., 2006, Astrophys. J., 652, 1246
  * Vladimirov, Ellison & Bykov (2008) Vladimirov, A., Ellison, D. C. & Bykov, A., 2008, Astrophys. J., 688, 1084
  * Zirakashvili & Aharonian (2010) Zirakashvili, V. N. & Aharonian, F. A., 2010, Astrophys. J., 708, 965

##  Appendix A Simulation box & parameters

With respect to the validity and consistency of verifying the previous dynamic
Monte Carlo simulation method, the present simulation program uses the same
simulation box and identical parameters as the previous dynamical Monte Carlo
simulation (Knerr, Jokipii & Ellison 1996). The schematic diagram of the
simulation box is shown in Figure 8 and all the simulation parameters are
listed in Table 2.

Figure 8:  Shock is produced by supersonic flow toward the stationary
reflective wall to the right. Continuous inflow of new particles occurs at the
left boundary. Inflow velocity is $U_{0}$, a particle’s total velocity is $V$,
the local frame particle velocity is $V_{L}$, and thermal velocity
$<V_{L}>=v_{0}$. The two circles represent one typical particle in the
upstream region and one in the downstream region, respectively. The vertical
solid line represents the shock front, the vertical dashed line represents the
FEB, the velocity of the shock is $V_{shock}$, size of the foreshock region is
$Xfeb=90$, the upstream flow velocity $U=U_{0}$, and the downstream flow
velocity is $U=0$. The magnetic field $B_{0}$ and inflow velocity $U_{0}$ are
both normal to the shock front (see Knerr, Jokipii & Ellison 1996). Table 2:
The parameters of the simulated cases

Inflow velocity |  $u_{0}$=0. 3 | 403km/s
---|---|---
Thermal speed |  $\upsilon_{0}$=0. 02 | 26\. 9km/s
Scattering time |  $\tau$=0. 833 | 0\. 13s
Box size |  $X_{max}$=300 | $10R_{e}$
Total time |  $T_{max}$=2400 | 6\. 3minutes
Time step size |  $dt$=1/30 | 0\. 0053s
Number of zones |  $nx$=600 | . . .
Initial particles per cell |  $n_{0}$=650 | . . .
FEB distance |  $X_{feb}$=90 | $3R_{e}$
111Scaling used a box size = 10$R_{e}$ (where $R_{e}$ represents the Earth’s
radius) and the box frame inflow velocity $u_{0}$= 403km/s. This implies the
following scale factors for distance, velocity, and time:
$X_{scale}$=10$R_{e}$/300, $v_{scale}$=403km $s^{-1}$/0.3, and
$t_{scale}$=$x_{scale}$/$v_{scale}$. Here, the Mach number M=11.6.
Dimensionless or normalized numbers are used in the text to describe our
simulations, except for specifically highlighted examples (see Knerr, Jokipii
& Ellison 1996).


scaled_IEEEeqnarray

(1)

# Integration of a systolic array based hardware accelerator into a DNN
operator auto-tuning framework
††thanks: The research leading to these results was supported by the
Competence Center Karlsruhe for AI Systems Engineering (CC-KING
https://www.ai-engineering.eu) sponsored by the Ministry of Economic Affairs,
Labour and Tourism Baden-Württemberg. It was conducted within the project
KIsSME (Artificial Intelligence for selective near-real-time recordings of
scenario and manoeuvre data in testing highly automated vehicles, Grant No.
19A20026C), funded by the German Federal Ministry for Economic Affairs and
Climate Action.

Federico Nicolás Peccia  FZI Research Center for Information Technology
Karlsruhe, Germany
<EMAIL_ADDRESS>Oliver Bringmann  University of Tübingen
Tübingen, Germany
<EMAIL_ADDRESS>
###### Abstract

The deployment of neural networks on heterogeneous SoCs coupled with custom
accelerators is a challenging task because of the lack of end-to-end software
tools provided for these systems. Moreover, the already available low-level
schedules and mapping strategies provided by the accelerator developers for
typical tensor operations are not necessarily the best possible ones for each
particular use case. This is why frameworks which automatically test the
performance of the generated code on a specific hardware configuration are of
special interest. In this work, the integration between the code generation
framework TVM and the systolic array-based accelerator Gemmini is presented. A
generic schedule to offload the GEneral Matrix Multiply (GEMM) tensor
operation onto Gemmini is detailed, and its suitability is tested by executing
the AutoTVM tuning process on it. Our generated code achieves a peak
throughput of 46 giga-operations per second (GOPs) under a 100 MHz clock on a
Xilinx ZCU102 FPGA, outperforming previous work. Furthermore, the code
generated by this integration was able to surpass the default hand-tuned
schedules provided by the Gemmini developers in real-world workloads.

###### Index Terms:

RISC-V, FPGA, TVM, Gemmini, Accelerator, Code Generation

##  I Introduction

Heterogeneous SoCs generators like Chipyard [1] or HEROv2 [2] are quickly
being adopted for an increasing amount of use cases thanks to their great
adaptability. This kind of generator exposes an enormous amount of possible
configurations to the user, thus enabling the generation of SoCs tailor-made
for specific applications. Thanks to the open nature of these projects,
hardware accelerators like Gemmini [3], VTA [4] or NVDLA [5] are being
developed to offload specific workloads from the CPU, allowing the deployment
of complex algorithms onto edge platforms.

Although the developers of these accelerators normally provide hand-tuned
kernels to offload commonly used tensor operations into the accelerator, the
high configurability of these SoC generators becomes a problem: there is no
guarantee that these default schedules will provide the best throughput across
all possible SoC configurations.

This is why automatic code generation and evaluation tools are becoming
increasingly popular [6, 7, 8]. These present enormous advantages, as they can
easily generate different mapping options for each tensor operator and test
them to automate the process of finding the best schedule parameters given a
particular SoC configuration (a process known as auto-tuning). By measuring on
a physical hardware platform, the impact of other system components on the
accelerator’s operation is also taken into account.

To demonstrate the advantages of this auto-tuning process for hardware
accelerators, this work presents the integration of the Gemmini accelerator
into the TVM deployment framework. The paper is organized as follows: first, a
scheduling search space for a GEMM operation on a generic systolic array
accelerator is proposed in Section II. Then, Section III presents the
integration of the Gemmini accelerator into the TVM framework. 111As of the
date of submission of this paper, the merging of this integration into the
main TVM branch is still in progress. Finally, Section IV compares the auto-
tuning of different GEMM workloads against a reference implementation [9], and
demonstrates that our scheduling definition achieves improvements in terms of
giga-operations per second (GOPs) for all workloads. The measurements are
expanded by presenting the autotuning results for operators extracted from the
Baidu DeepBench Workloads, outperforming the default handcrafted kernels
provided by the Gemmini developers.

(a)

1:config.accel()

2:for $i_{o}=0$ to $(M/tile_{m_{1}})-1$ do

3: for $j_{o}=0$ to $(N/tile_{n_{1}})-1$ do

4: move.in($D^{\prime}$)

5: for $k_{o}=0$ to $(K/tile_{k_{1}})-1$ do

6: move.in($A^{\prime}$)

7: move.in($B^{\prime}$)

8: for $i_{i}=0$ to $(tile_{m_{1}}/tile_{m_{2}})-1$ do

9: for $j_{i}=0$ to $(tile_{n_{1}}/tile_{n_{2}})-1$ do

10: for $k_{i}=0$ to $(tile_{k_{1}}/tile_{k_{2}})-1$ do

11: gemm.DIMxDIM($C^{\prime}$,$A^{\prime}$,$B^{\prime}$)

12: end for

13: end for

14: end for

15: end for

16: for $i_{i}=0$ to $(tile_{m_{1}}/tile_{m_{2}})-1$ do

17: for $j_{i}=0$ to $(tile_{n_{1}}/tile_{n_{2}})-1$ do

18: move.out($C^{\prime}$)

19: end for

20: end for

21: end for

22:end for

(b)

(c)

Figure 1: (1(a)) shows the proposed GEMM schedule parameters for an
accelerator based on a $DIM\times DIM$ systolic array able to execute a
generic GEMM with form $C=A\times B+D$. (1(b)) shows an example generated
pseudocode for the operation, and (1(c)) shows a graphical representation of
the move of data in and out of the accelerator.

##  II Scheduling a GEMM operation on a systolic array

In 1982, Kung [10] presented the advantages of systolic array architectures
for the execution of matrix operations. These have been especially attractive
to accelerate neural network operators because of their high data reuse.
Several accelerators were built using systolic arrays (or similar processing
element distributions) as its core [11, 12, 13, 14, 15, 16].

To correctly schedule a tensor operation on a hardware accelerator, one should
be aware of the kind of instructions provided by the accelerator. These can be
classified into two distinct categories:

  * •

Low level: these instructions allow the programmer fine control of the
behaviour of the accelerator. These are typically instructions that enable
moves of data in and out of the accelerator’s internal memory, and others that
execute/dispatch the most basic supported computation. An example of this kind
of instructions can be found on the Angel-Eye [17] or the VTA [4]
accelerators.

  * •

Layer wide: these instructions take care of the burden of managing the
individual instructions, by exposing to the programmer a higher-level
interface to execute an entire tensor operation, like the NVDLA [5].

For a systolic array-based accelerator, we are interested in the first kind of
instructions, because they allow the programmer or the tuning software fine-
grained control over the generated schedules. To correctly schedule them, the
code generator framework should take these 4 factors into account:

  1. 1.

Configuration of the hardware: these instructions should be generated only
when the configuration of the accelerator changes, and not on each new
operator, to avoid unnecessary reconfigurations.

  2. 2.

Move of data into the accelerator’s SRAM: if the input matrices fit entirely
on the accelerator’s internal memory, one could choose to first move the
entire matrices in, and only then start to generate the compute instructions.
But perhaps interleaving move and compute instructions can actually avoid idle
times and thus generate faster compute schedules. The proposed schedule should
be able to achieve this load balancing between different kinds of
instructions.

  3. 3.

Computation: the idle time of the systolic array should be minimized. There
should always be data available in the accelerator’s SRAM for the systolic
array to use for its computation.

  4. 4.

Move of data out into the external DRAM: two different options could be chosen
when generating these instructions: either patches of results are moved out as
soon as they are ready, or multiple patches are stored in the accumulator
before bulk moving all of them out.

Figure 2: Integration workflow example for a neural network formed by a fully
connected layer and a softmax layer

To limit the explosion of the schedule parameter search space, certain
hardware limitations and basic assumptions can be taken into account during
this stage:

  * •

Schedules which do not respect the maximum limitation of columns and rows for
each move should automatically be dropped.

  * •

The schedule cannot generate configurations where data would overflow the
accelerator’s SRAM: hardware information should be used as a limit for the
generated move instructions source and destination addresses.

  * •

Each generated GEMM should try to utilize the systolic array to its full
extent, to prevent idle processing elements.

The configuration parameters of the proposed schedule presented in Fig. 1
along with the accelerator parameters can effectively generate code that
considers all the previous mentioned points:

  * •

$tile_{m,n,k}$: two-level tiling for each corresponding computation axis. The
first level moves data into the accelerator’s memory, and the second
partitions the moved data into smaller GEMMs to fit the systolic array, trying
to maximize its usage.

  * •

$parallel\\_accumulations$: denotes how many output patches are accumulated
simultaneously in the accumulator, before moving them out. Modifies the
position of the move.out instruction in the generated code.

  * •

$apply\\_double\\_buffer$: represents if double buffering should be applied on
the move.in of data, weights, both or none. Takes into account the bank size
of the scratchpad, so that double-buffered data is written to different banks,
thus preventing bank access conflict.

  * •

$exchange\\_axis$: allows to reorder the computation of axis M and N, useful
when M and N are not equal.

  * •

$WS/OS$: if supported by the accelerator, configures the systolic array to
work on output stationary mode or weight stationary mode [18].

  * •

$mvout\\_big\\_block$: if supported by the accelerator, enables the generation
of move.out instructions which move more than $DIM\times DIM$ size patches.
Useful to explore the trade-off between burst of smaller move.out instructions
versus bigger instructions.

##  III Integrating Gemmini into TVM

To validate the aforementioned schedule search space, the TVM [7] framework
was chosen to take advantage of its AutoTVM module. This auto-tuning process
goes over the parameter search space, measures each configuration on physical
hardware, and exports the best schedule configuration. To avoid an infeasible
amount of measurements, XGB model-based tuners are available [19].

The Gemmini [3] systolic array-based accelerator was selected as a test
platform because of its open-source nature. Gemmini was developed by the UC
Berkeley and is part of the Chipyard ecosystem. It works in a tightly-coupled
manner with a RISC-V CPU, using the Rocket Chip Coprocessor (RoCC) interface
to control the accelerator with help from custom instructions.

Gemmini uses a systolic array of $DIM\times DIM$ multiply-and-accumulate (MAC)
processing elements to perform matrix multiplications. The data is consumed
from a scratchpad made up of banked SRAMs and is stored in a series of banked
SRAMs equipped with adder units known as the accumulator. A DMA engine
connected to the System Bus (directly to the L2 cache) is used to get data in
and out of the accelerator’s SRAMs. Gemmini is also able to apply scaling
factors during the move in and out of data, and also other common operations
like ReLu or max pooling.

Gemmini’s RoCC instructions can be grouped into the following categories:

  * •

Configuration: these instructions configure the input, output and execution
pipelines of the accelerator.

  * •

Move: these instructions move specific amounts of rows and columns of data
into the SRAM or out into the DRAM.

  * •

Execution: dispatch the actual execution of the $DIM\times DIM$ GEMM to the
systolic array.

  * •

Flush and fence: general maintenance instructions.

  * •

Loop: ”CISC” type instructions for commonly used operations. They take away
the burden of manually scheduling the intrinsic instructions by generating
them directly on the hardware using FSMs. Gemmini developers also provide some
handcrafted layer-wise C functions, which internally call these loop
instructions.

Gemmini uses a decoupled access-execute architecture, with dedicated
controllers to manage the move.in, execute and move.out instructions
independently. A ROB is included to detect hazards between instructions and to
issue them to their respective controller.

Fig. 2 shows how the developed integration between TVM and Gemmini works.
First, the TensorFlow Lite quantized model is imported. Here, the model is
transformed into the Relay IR dataflow graph representation of TVM. Then,
pattern matching is used to replace subgraphs of operations with custom
operators supported by Gemmini. These operators consist of a computation
definition and its schedule: a set of parametrized loop transformations. For
the GEMM operator, two different schedules were developed: one that generates
the calls to the intrinsic instructions parametrized as presented in Section
II, and one that replaces the entire loops with a call to the default
handcrafted layer-wise C function.

During the schedule transformation, several pragmas are used to tag specific
loops and then replace them with the Gemmini intrinsic instructions in the
following compilation pass. The tensorization feature of TVM was used to
insert the two instructions (preload and compute) needed by Gemmini to compute
a GEMM. TVM’s C code generator was used to create the source code file that
executes the operators, and the standard header file provided by Gemmini was
included in the generated file by TVM to correctly reference Gemmini’s
instruction macros.

1632641282565121024$0$$20$$40$$60$Workload size GOPs [1] no L2[1] L2Ours no
L2Ours L2 (a) Best schedules found by AutoTVM in our implementation against
previous work.

1632641282565121024$0$$20$$40$$60$Workload size GOPs CISC no L2CISC L2Ours no
L2Ours L2 (b) Best schedules found by AutoTVM in our implementation against
the default ”CISC instruction” based schedules.

Figure 3: Results across different GEMM workloads. For each workload,
$M=N=K=workload\>size$

###  III-A Quantization management

TensorFlow Lite’s quantization scheme [20] uses symmetric quantization for the
weights and biases of each layer, and asymmetric for the input and output of
each layer 222see
https://www.tensorflow.org/lite/performance/quantization_spec?hl=en. Because
the multiplication of two quantized matrices with different zero points is not
straightforward, correction terms have to be subtracted to get the correct
result as seen in Eq. (2). Finally, to transform $Q^{\prime}_{C}$ into the
output quantization regime of the layer, a requantization operator is inserted
by TVM described by Eq. (3).

An easy approach would be to accelerate only the matrix multiplication using
the systolic array (the term 1 of Eq. (2)) and then execute the correction
terms subtraction and the re-quantization operator on the CPU. But this is
slow and doesn’t exploit the benefits of TVM’s compilation framework, so a
transformation of the GEMM operation is done during the compilation pass,
allowing TVM to fold the remaining constants into the bias of the layer. The
final proposed solution for a quantized matrix multiplication with bias
addition is described in Eq. (4) and (5). The scaling factor $s_{d}/s_{c}$ is
inserted as output scaling factor in the configuration of the move.out
instruction of Gemmini, thus achieving the acceleration of the entire
operator, without extra tensor operations executed on the CPU.

$Q^{\prime}_{C_{(m,n)}}=\sum_{0}^{k}\left(Q_{A_{(m,k)}}-zp_{a}\right)*Q_{B_{(k,n)}}=\\\
\sum_{0}^{k}Q_{A_{(m,k)}}*Q_{B_{(k,n)}}-\sum_{0}^{k}zp_{a}*Q_{B_{(k,n)}}$ (2)

$\displaystyle Q_{C}=zp_{c}+\frac{s_{d}}{s_{c}}*Q^{\prime}_{C}$ (3)
$\displaystyle
Q_{C_{(m,n)}}=\frac{s_{d}}{s_{c}}\left[\left(\sum_{0}^{k}Q_{A_{(m,k)}}*Q_{B_{(k,n)}}\right)+Q^{\prime}_{D_{(m,n)}}\right]$
(4) $\displaystyle
Q^{\prime}_{D_{(m,n)}}=Q_{D_{(m,n)}}-\sum_{0}^{k}zp_{a}*Q_{B_{(k,n)}}+\frac{s_{d}}{s_{c}}*zp_{c}$
(5)

##  IV Experiments

Table I: Selected Baidu DeepBench workloads Id | M | N | K
---|---|---|---
15 | 64 | 1 | 1216
49 | 128 | 1 | 1024
63 | 512 | 1 | 512
73 | 512 | 2 | 512
84 | 1024 | 4 | 512
1549637384$0$$2$$4$$6$$8$Id GOPs Default no L2Default L2Ours no L2Ours L2
Figure 4: Best schedules found by AutoTVM for the Baidu DeepBench dataset
using our implementation.

In order to be able to compare our results against the ones reported by [9],
our test setup was also built using the Rocket Core [21] together with a
$16\times 16$ Gemmini accelerator, with a 256 KB scratchpad (4 SRAM banks with
4096 rows each) and a 64 KB accumulator (1 SRAM bank with 1024 rows). The
inputs and weights are represented using 8 bits and the accumulated values
using 32 bits. In the implemented design which uses the L2 cache, the SiFive
inclusive L2 cache was used. Like the original paper, the resulting hardware
was also implemented for a Xilinx Zynq UltraScale+ ZCU102 FPGA running at 100
MHz. The only difference between our setup and the one from [9] is the size of
the L2 cache because they did not report the selected size. In our work, the
default size used in the Chipyard project was selected: 512 KB.

Fig. 3 presents the best schedule performance found by AutoTVM, and compares
it with the previous work and the default Gemmini schedules. An XGB tuner [19]
with early stopping equal to 500 iterations was used to traverse the search
space. The amount of operations for a generic GEMM of form
$C_{[M,N]}=A_{[M,K]}\times B_{[K,N]}+D_{[M,N]}$ was defined as $OP=2\times
M\times N\times K+M\times N$.

Although Fig. 3(b) shows that our AutoTVM schedules are better than the CISC-
based schedules for almost all workloads, it fails to find better schedules
for workloads 256 and 512 when the L2 cache is activated. There are two
possible explanations for this behaviour. The first one would be that the XGB
tuner is stuck in a local optimum, but this was verified to not be the case by
executing an exhaustive grid-search of all the schedule parameter search
space: the XGB found schedules are indeed the best possible schedules our
implementation can generate for that workloads. The second possible
explanation lies in a load balancing feature of the Gemmini’s CISC
instructions. The hardware FSMs monitor the proportion of each instruction
type in the ROB and can pause the generation of each kind of instruction, to
maximize the overlap between move and execute operations. This behaviour can
not yet be replicated using our TVM integration, and further work needs to be
done to analyse how to implement a similar feature.

To show the effectiveness of this schedule search space on real-world
examples, a set of dense layer workflows taken from the Baidu DeepBench
dataset [22] were selected (Table I). Fig. 4 presents the result of the
AutoTVM tuning process executed on these workloads, using the same tuner
parameters and FPGA bitstreams as in the previous measurements.

##  V Conclusions

This work proposed a schedule parameter space for a GEMM tensor operation. The
proposed schedule was configured into the TVM deep learning compiler, and
integrated with the Gemmini systolic array hardware accelerator. We
demonstrate that this schedule parameter space allows the autotuning process
to find faster schedules than previous work, and also faster schedules than
the hardcoded schedules provided by the expert developers of Gemmini on almost
all tested workloads.

Future works will also add other operators to the integration, like
convolutions and depthwise convolutions. The results of the acceleration of
entire neural networks will also be presented.

This work did not try to improve how much time the autotuning process takes.
The search space of scheduling parameters was traversed using an XGB tuner
approach, which can be improved to converge faster as shown in [23]. In future
works, optimization techniques to speed up this process should be
investigated, to be able to find the best scheduling parameters using the
minimum amount of actual hardware measurements.

## References

  * [1] A. Amid, D. Biancolin, A. Gonzalez, D. Grubb, S. Karandikar, H. Liew, A. Magyar, H. Mao, A. Ou, N. Pemberton, P. Rigge, C. Schmidt, J. Wright, J. Zhao, Y. S. Shao, K. Asanović, and B. Nikolić, “Chipyard: Integrated design, simulation, and implementation framework for custom socs,” _IEEE Micro_ , vol. 40, no. 4, pp. 10–21, 2020.
  * [2] A. Kurth, B. Forsberg, and L. Benini, “Herov2: Full-stack open-source research platform for heterogeneous computing,” _CoRR_ , vol. abs/2201.03861, 2022\. [Online]. Available: https://arxiv.org/abs/2201.03861
  * [3] H. Genc, S. Kim, A. Amid, A. Haj-Ali, V. Iyer, P. Prakash, J. Zhao, D. Grubb, H. Liew, H. Mao, A. Ou, C. Schmidt, S. Steffl, J. Wright, I. Stoica, J. Ragan-Kelley, K. Asanovic, B. Nikolic, and Y. S. Shao, “Gemmini: Enabling systematic deep-learning architecture evaluation via full-stack integration,” in _Proceedings of the 58th Annual Design Automation Conference (DAC)_ , 2021.
  * [4] T. Moreau, T. Chen, L. Vega, J. Roesch, E. Yan, L. Zheng, J. Fromm, Z. Jiang, L. Ceze, C. Guestrin, and A. Krishnamurthy, “A hardware-software blueprint for flexible deep learning specialization,” 7 2018. [Online]. Available: http://arxiv.org/abs/1807.04188
  * [5] G. Zhou, J. Zhou, and H. Lin, “Research on nvidia deep learning accelerator,” in _2018 12th IEEE International Conference on Anti-counterfeiting, Security, and Identification (ASID)_ , 2018, pp. 192–195.
  * [6] L. Zheng, C. Jia, M. Sun, Z. Wu, C. H. Yu, A. Haj-Ali, Y. Wang, J. Yang, D. Zhuo, K. Sen, J. E. Gonzalez, and I. Stoica, “Ansor : Generating high-performance tensor programs for deep learning,” _CoRR_ , vol. abs/2006.06762, 2020. [Online]. Available: https://arxiv.org/abs/2006.06762
  * [7] T. Chen, T. Moreau, Z. Jiang, H. Shen, E. Q. Yan, L. Wang, Y. Hu, L. Ceze, C. Guestrin, and A. Krishnamurthy, “TVM: end-to-end optimization stack for deep learning,” _CoRR_ , vol. abs/1802.04799, 2018. [Online]. Available: http://arxiv.org/abs/1802.04799
  * [8] S. Zheng, Y. Liang, S. Wang, R. Chen, and K. Sheng, “Flextensor: An automatic schedule exploration and optimization framework for tensor computation on heterogeneous system,” in _Proceedings of the Twenty-Fifth International Conference on Architectural Support for Programming Languages and Operating Systems_ , ser. ASPLOS ’20. New York, NY, USA: Association for Computing Machinery, 2020, p. 859–873. [Online]. Available: https://doi.org/10.1145/3373376.3378508
  * [9] P. Xu and Y. Liang, “Automatic code generation for rocket chip rocc accelerators,” _CARRV_ , 2020.
  * [10] Kung, “Why systolic architectures?” _Computer_ , vol. 15, pp. 37–46, 1 1982\. [Online]. Available: http://ieeexplore.ieee.org/document/1653825/
  * [11] Y. H. Chen, T. Krishna, J. S. Emer, and V. Sze, “Eyeriss: An energy-efficient reconfigurable accelerator for deep convolutional neural networks,” _IEEE Journal of Solid-State Circuits_ , vol. 52, 2017.
  * [12] Y. H. Chen, T. J. Yang, J. S. Emer, and V. Sze, “Eyeriss v2: A flexible accelerator for emerging deep neural networks on mobile devices,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , vol. 9, pp. 292–308, 6 2019.
  * [13] Z.-G. Liu, P. N. Whatmough, and M. Mattina, “Systolic tensor array: An efficient structured-sparse gemm accelerator for mobile cnn inference,” 5 2020\. [Online]. Available: http://arxiv.org/abs/2005.08098
  * [14] X. Wei, C. H. Yu, P. Zhang, Y. Chen, Y. Wang, H. Hu, Y. Liang, and J. Cong, “Automated systolic array architecture synthesis for high throughput cnn inference on fpgas,” vol. Part 128280, 2017.
  * [15] R. Xu, S. Ma, Y. Wang, X. Chen, and Y. Guo, “Configurable multi-directional systolic array architecture for convolutional neural networks,” _ACM Transactions on Architecture and Code Optimization_ , vol. 18, 2021.
  * [16] N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, A. Borchers, R. Boyle, P. luc Cantin, C. Chao, C. Clark, J. Coriell, M. Daley, M. Dau, J. Dean, B. Gelb, T. V. Ghaemmaghami, R. Gottipati, W. Gulland, R. Hagmann, C. R. Ho, D. Hogberg, J. Hu, R. Hundt, D. Hurt, J. Ibarz, A. Jaffey, A. Jaworski, A. Kaplan, H. Khaitan, D. Killebrew, A. Koch, N. Kumar, S. Lacy, J. Laudon, J. Law, D. Le, C. Leary, Z. Liu, K. Lucke, A. Lundin, G. MacKean, A. Maggiore, M. Mahony, K. Miller, R. Nagarajan, R. Narayanaswami, R. Ni, K. Nix, T. Norrie, M. Omernick, N. Penukonda, A. Phelps, J. Ross, M. Ross, A. Salek, E. Samadiani, C. Severn, G. Sizikov, M. Snelham, J. Souter, D. Steinberg, A. Swing, M. Tan, G. Thorson, B. Tian, H. Toma, E. Tuttle, V. Vasudevan, R. Walter, W. Wang, E. Wilcox, and D. H. Yoon, “In-datacenter performance analysis of a tensor processing unit,” _ACM SIGARCH Computer Architecture News_ , vol. 45, 2017\.
  * [17] K. Guo, L. Sui, J. Qiu, J. Yu, J. Wang, S. Yao, S. Han, Y. Wang, and H. Yang, “Angel-eye: A complete design flow for mapping cnn onto embedded fpga,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 37, 2018.
  * [18] A. Samajdar, Y. Zhu, P. N. Whatmough, M. Mattina, and T. Krishna, “Scale-sim: Systolic CNN accelerator,” _CoRR_ , vol. abs/1811.02883, 2018. [Online]. Available: http://arxiv.org/abs/1811.02883
  * [19] T. Chen, L. Zheng, E. Yan, Z. Jiang, T. Moreau, L. Ceze, C. Guestrin, and A. Krishnamurthy, “Learning to optimize tensor programs,” vol. 2018-December, 2018.
  * [20] B. Jacob, S. Kligys, B. Chen, M. Zhu, M. Tang, A. G. Howard, H. Adam, and D. Kalenichenko, “Quantization and training of neural networks for efficient integer-arithmetic-only inference,” _CoRR_ , vol. abs/1712.05877, 2017. [Online]. Available: http://arxiv.org/abs/1712.05877
  * [21] K. Asanović, R. Avizienis, J. Bachrach, S. Beamer, D. Biancolin, C. Celio, H. Cook, D. Dabbelt, J. Hauser, A. Izraelevitz, S. Karandikar, B. Keller, D. Kim, J. Koenig, Y. Lee, E. Love, M. Maas, A. Magyar, H. Mao, M. Moreto, A. Ou, D. A. Patterson, B. Richards, C. Schmidt, S. Twigg, H. Vo, and A. Waterman, “The rocket chip generator,” EECS Department, University of California, Berkeley, Tech. Rep. UCB/EECS-2016-17, Apr 2016. [Online]. Available: http://www2.eecs.berkeley.edu/Pubs/TechRpts/2016/EECS-2016-17.html
  * [22] S. Narang and G. Diamos, “Baidu deepbench,” _GitHub Repository_ , 2017.
  * [23] D. Rieber, M. Reiber, O. Bringmann, and H. Fröning, “Hw-aware initialization of dnn auto-tuning to improve exploration time and robustness,” 5 2022. [Online]. Available: http://arxiv.org/abs/2205.15568