# Lagrangian Numerical Methods for Ocean Biogeochemical Simulations

Francesco Paparella1,2
Marina Popolizio2

###### Abstract

We propose two closely–related Lagrangian numerical methods for the simulation
of physical processes involving advection, reaction and diffusion. The methods
are intended to be used in settings where the flow is nearly incompressible
and the Péclet numbers are so high that resolving all the scales of motion is
unfeasible. This is commonplace in ocean flows. Our methods consist in
augmenting the method of characteristics, which is suitable for
advection–reaction problems, with couplings among nearby particles, producing
fluxes that mimic diffusion, or unresolved small-scale transport. The methods
conserve mass, obey the maximum principle, and allow to tune the strength of
the diffusive terms down to zero, while avoiding unwanted numerical
dissipation effects.

1Division of Sciences and Mathematics
New York University Abu Dhabi

2Dipartimento di Matematica e Fisica
Università del Salento Lecce

Keywords: Ocean biogeochemistry; lagrangian methods; advection reaction
diffusion; unresolved flows.

##  1 Introduction

Biogeochemical problems in oceanography are usually expressed in terms of
coupled advection-reaction-diffusion equations involving scalar fields,
sometimes in large number, representing chemical species, biological species,
or functional groups (see, e.g., [1]). These fields are advected by the ocean
currents, are subject to diffusion, and interact nonlinearly with each other.

A generic, abstract form of oceanographical biogeochemical equations is the
following

$\displaystyle\frac{\partial c_{1}}{\partial t}+\boldsymbol{u}\cdot\nabla
c_{1}=$ $\displaystyle D_{1}\nabla^{2}c_{1}+f_{1}(c_{1},\ldots,c_{n})$
$\displaystyle\vdots$ (1) $\displaystyle\frac{\partial c_{n}}{\partial
t}+\boldsymbol{u}\cdot\nabla c_{n}=$ $\displaystyle
D_{n}\nabla^{2}c_{n}+f_{n}(c_{1},\ldots,c_{n})$

where $c_{1},\ldots,c_{n}$ are the scalar fields, $\boldsymbol{u}$ is the
water velocity field in the region of interest, which is assumed to be known,
$D_{1},\ldots,D_{n}$ are the diffusion coefficients, and the functions
$f_{1},\ldots,f_{n}$ specify the local interactions among the scalar fields.

The relative importance of the transport and diffusion terms is quantified by
the Péclet numbers

$Pe_{l}=\frac{UL}{D_{l}}$

where $U$ and $L$ are, respectively, a characteristic speed and a
characteristic length associated to the velocity field $\boldsymbol{u}$. The
relative importance of the transport and reaction terms is quantified by the
Damköhler numbers

$Da_{l}=\frac{L}{U\tau_{l}}$

where $\tau_{l}$ is a characteristic time scale associated with the reaction
described by $f_{l}$.

The Damköhler number for phytoplankton may range from negligibly small up to
$O(10)$ [2]. While large values of the Damköhler number may amplify the
patchiness of a reacting scalar as compared to a non–reacting one [3, 4] and
make the problem stiff, the true source of numerical difficulties in
biogeochemical applications lies in the enormous size of the Péclet number.

If one takes the diffusivities to be the molecular ones (or computed from the
mean square displacement of trajectories of individual plankton cells) then
the Péclet numbers may easily exceed $10^{10}$. Such a large value is
reflected in the fact that ocean tracers (temperature, salinity, etc.) show
structures from the scale of ocean basins down to submillimetric scales. Even
accounting for a continuing rapid pace of improvement in computer
technologies, it is quite obvious that, in the foreseeable future, no
numerical code will be able to resolve such a wide interval of scales.

In the absence of reaction terms, a reasonable way to deal with unresolved
small scales is to parameterize the advective fluxes due to the unresolved
scales with diffusion operators (often in a more complicated form than simple
Laplacians). To this end there is an impressive array of techniques, ranging
from explicitly adding new terms to the equations (e.g. in turbulence
closures), to using flux or slope limiters (e.g. in finite volume methods), to
advection and interpolation (e.g. in semi-lagrangian methods) or dealiasing
and filters (e.g. in pseudo–spectral methods). A review of numerical methods
used for geophysical flows is given in [5]. In all these cases, however, the
strength of the diffusive terms is determined not just by the physical
parameters of the problem, but also by the size of the mesh. In fact, all
these techniques may be viewed as different ways to average out the subgrid
scales. Thus, in the presence of unresolved small scales, the values of the
scalar fields at each grid node must be understood not as a pointwise
evaluation of a function, but as an average over a spatial region having an
extension comparable with the size of a computational mesh.

Early studies already showed that changing the strength of the diffusive
fluxes representing the unresolved scales may have a dramatic impact on the
reaction terms [2, 6] and warned that a “mean field” approach might be
inappropriate for modeling plankton dynamics. Later studies, conducted using
realistic ocean models, showed strong fluctuations in plankton productivity
depending on the advection scheme used and, most importantly, on the
resolution [7, 8, 9]. The most recent assessment of the importance of the
unresolved structures is found in [10].

As a first step to understand these results we need to observe that, for the
full set of equations (1), one faces the overwhelming difficulty that an
averaging operator $\overline{}$ does not commute with nonlinear reaction
terms:
$f_{l}(\overline{c_{1}},\ldots,\overline{c_{n}})\neq\overline{f_{l}(c_{1},\ldots,c_{n})}$.
Because reactions are formally pointwise one would need to compute
$\overline{f_{l}(c_{1},\ldots,c_{n})}$, but all that current grid–based codes
can do is to compute $f_{l}(\overline{c_{1}},\ldots,\overline{c_{n}})$. The
wide chasm of unresolved scales means that the mesh-averaged values
$\overline{c_{1}},\ldots,\overline{c_{n}}$ may be substantially different from
their pointwise counterpart $c_{1},\ldots,c_{n}$. As we shall see in the
following, the bias produced by this effect may have either sign, depending,
among other things, on the initial conditions.

In the absence of any diffusive effect, that is, setting $D_{1,\ldots,n}=0$ in
(1), it is arguably better to avoid any discretization involving an Eulerian
grid, and use a straightforward implementation of the method of
characteristics. This leads to the following simple Lagrangian numerical
scheme: we uniformly seed the domain $\Omega$ with $M$ particles, having
position $\boldsymbol{x}_{i}$, $i=1,\ldots,M$, and then numerically solve

$\begin{cases}\dot{\boldsymbol{x}}_{i}&=\boldsymbol{u}(\boldsymbol{x}_{i},t)\\\
\dot{c}_{1;i}&=f_{1}(c_{1;i},\ldots,c_{n;i})\\\ &\vdots\\\
\dot{c}_{n;i}&=f_{n}(c_{1;i},\ldots,c_{n;i})\end{cases}$ (2)

with one among many viable ODE solvers. Here and in the following we use the
shorthand notation $c_{l;i}=c_{l}(\boldsymbol{x}_{i},t)$ for the scalars
sampled at the location of each particle. It is important to appreciate that,
even when the number of particles is too small to fully sample the small-scale
structures present in the full solution of the PDEs, the values $c_{l;i}$
remain unaffected by the sparsity of the sampling, and are only affected by
inaccuracies in the solution of the ODEs (2), due, e.g., to an imperfect
knowledge of the velocity field $\boldsymbol{u}$. This scheme is thus immune
from the averaging problem discussed above. If, as is the case in
oceanographic applications, the velocity field $\boldsymbol{u}$ is
divergenceless, or nearly so, then an initially uniform sampling will remain
uniform, or nearly so, at all future times. In this context the lack of a
structured grid is just a nuisance: diagnostic and data analysis tasks may be
performed after resampling the numerical solutions of (2) on a regular grid of
choice, using, e.g., the methods discussed in [11, §5.3].

Unfortunately, the method of characteristics is not directly applicable to
biogeochemical problems: the complete absence of diffusive effects in (2)
would lead to paradoxical effects. For instance, if a water mass containing
some phytoplankton but poor of nutrients were brought close to water masses
devoid of phytoplankton but nutrient–rich, fluxes associated to small–scale
motions would seed some plankton in the nutrient–rich water masses, leading,
if the conditions are right, to a bloom. With the scheme (2) a particle full
of phytoplankton could be brought arbitrarily close to a particle full of
nutrients and yet there would be no exchanges between the two: the plankton
would wither, and the nutrients would remain unused.

In this paper we show how to augment the simple Lagrangian scheme (2) with
couplings among nearby particles designed to mimic diffusive effects or, more
generally, fluxes due to small-scale, unresolved transport processes. In order
to be acceptable, such a coupling must possess the following three properties

  1. 1.

respect mass conservation;

  2. 2.

obey the maximum principle;

  3. 3.

allow to recover the scheme (2) in the limit $D_{l}\to 0$.

The importance of mass conservation is fairly obvious. Even for models using
non–conserving reaction terms, there is no reason to introduce uncontrollable
numerical sources and sinks of scalars. Schemes that do not obey the maximum
principle may create maxima and minima unbounded by the maxima and minima of
the initial conditions. In particular, scalar fields that should be non-
negative (e.g. the concentration of a chemical species) may locally develop
negative values, which, in turn, yield meaningless results with most reaction
models. Being able to recover the scheme (2) means that one is free to tune
the strength of the diffusive effects on the basis of modeling considerations
alone, and not because of numerical requirements. We propose two distinct
couplers that satisfy all these three properties. Of the two methods that we
propose, the first is based on an integral formulation, the second is an
heuristic recipe based on physical considerations. The two methods are
distinct in the way used to enforce mass conservation. In both cases, however,
the maximum principle is a direct consequence of the fact that the
concentration of each particle after a diffusive step is determined as an
average involving the concentrations of nearby particles. Free parameters,
appearing in both methods, can be use to tune the strength of the diffusive
effects to extremely low values, or to zero, thereby maintaining the particles
uncoupled.

Particle–based methods are not a novelty. Smoothed particle hydrodynamics
(SPH) has proved to be very suitable for highly compressible astrophysical
problems, but flexible enough to be applied in many other settings [12],
including heat conduction [13]. However, we felt that achieving all three of
the above properties might be not straightforward with an SPH–inspired
approach, therefore our methods are not based upon the differentiation of a
smooth kernel. Other particle–based methods, closer to the spirit of the
present work, have been proposed for diffusion and advection–diffusion
equations [14, 15], but did not gain a large popularity.

Few are the instances in which Lagrangian methods have been applied to
geophysical problems. Nearly all numerical ocean models use grid–based
methods, with the notable exception of the so–called “slippery sack” model
[16]. This was initially a purely adiabatic, Lagrangian scheme, which was
later augmented with a diffusive coupling between nearby particles [17]. The
Lagrangian scheme (2) has been successfully applied to explain some
incongruences between ecological models and observations [18]. When augmented
with a diffusive coupling it has been used to explain the Fourier spectrum of
a plankton concentration field [19]. We are not aware of other applications of
Lagrangian schemes to ocean biogeochemistry. There exists more work on
Lagrangian methods for modeling the atmosphere. In particular, a method based
on contour advection and surgery has been highly successful in reproducing the
observed distribution of stratospheric ozone [20, 21]. For a recent survey on
Lagrangian methods in atmospheric sciences see [22].

The enormous potential of diffusively–coupled Lagrangian methods in
biogeochemistry is illustrated by a simple example, inspired by the results
obtained with a much more realistic model in [23]. In (1) we set $n=2$ and
choose a two–dimensional, incompressible velocity field
$\boldsymbol{u}=(-\psi_{y},\psi_{x})$ defined through the streamfunction
$\psi(x,y)=\sin(x)\sin(y)$ on the doubly–periodic domain
$(x,y)\in[0,2\pi)\times[0,2\pi$). The reaction terms are

$f_{1}(c_{1},c_{2})=-r\,c_{1}c_{2},\quad f_{2}(c_{1},c_{2})=+r\,c_{1}c_{2},$
(3)

with $r=0.2$. We may see the scalar field $c_{2}$ as the spatial density of a
consumer that grows at the expense of a resource whose density is $c_{1}$. The
initial conditions are:

$c_{1}(x,y,0)=\cos^{2}(x/2),\quad c_{2}(x,y,0)=10^{-4}.$ (4)

We compute six solutions of this problem for progressively smaller
diffusivities and correspondingly higher resolutions. The six meshes have
$128\cdot 2^{k}$ points in each direction, and the diffusivities are
$D_{1}=D_{2}=10^{-3}\cdot 2^{-2k}$, $k=0,\ldots,5$. At each resolution, using
substantially lower diffusivities would lead to severe oscillations and
numerical instabilities. The solid lines in Figure 1A show the time evolution
of the spatial average of $c_{2}$ (that is, the mean consumer density). The
dots show the same quantity computed by using the Lagrangian scheme (2)
augmented with one of the two diffusive couplers that will be presented in the
following (namely, that of section 2.2). The six Lagrangian solutions all use
just $128^{2}$ particles, and they differ only in the strength of the
diffusive coupling.

Figure 1: Spatial average of the field $c_{2}$ as a function of time. The
solid lines are results obtained with a pseudo–spectral code, with
progressively higher resolution and correspondingly lower diffusivity (see
text). The dots are results obtained with a Lagrangian code using the coupler
of sec. 2.2 with $128^{2}$ particles, where the strength of the diffusive
coupling between particles is set as to match that of the pseudo–spectral
computations. Panel A: calculations starting from the initial condition (4).
Panel B: calculations starting from the initial condition (5). Figure 2: Field
$c_{2}$ at time $t=100$. One quarter of the whole domain is shown. A)
pseudo–spectral scheme on a $128\times 128$ points grid. B) pseudo–spectral
scheme on a $4096\times 4096$ points grid. C) Lagrangian scheme with $128^{2}$
particles and a diffusion matching that of A). D) Lagrangian scheme with
$128^{2}$ particles and a diffusion matching that of D). All cases use the
initial conditions (5).

In this particular example, because of the quadratic nonlinearity, the same
amount of resource $c_{1}$ is consumed faster if it is spatially concentrated
than if it is spread out on a larger surface but at lower concentrations. Thus
smaller diffusivities, which better preserve the concentration peaks of the
resource, yield a faster growth of the spatially averaged field $c_{2}$. In
other words, they yield a higher productivity of the consumer.

One might then be lead to hope that, just as unresolved turbulence can be
usefully approximated by effective diffusion terms, in the same way effective
reaction terms should be sought, representing the large–scale effects of the
small–scale chemistry, with parameters tuned as a function of the resolution
of the model. Here we give an example showing that this hope is unlikely to be
fulfilled: we just change the initial conditions (4) with

$c_{1}(x,y,0)=\left(\sin\left(\frac{x}{2}\right)\sin\left(\frac{y}{2}\right)\right)^{4},\quad
c_{2}(x,y,0)=\left(\cos\left(\frac{x}{2}\right)\cos\left(\frac{y}{2}\right)\right),^{4}$
(5)

and repeat the same calculations described above. Because the resource and the
consumer are now initially segregated into two nearly non–overlapping blobs,
larger diffusivities bring in contact the resource and the consumer more
quickly. As a result, we obtain the opposite effect as before: the growth of
the spatially averaged consumer is fastest at the lowest resolution, and
declines as the resolution is increased (Figure 1B). Thus, hypothetical
effective reaction terms intended to reproduce at low resolution the results
obtained at highest resolution with the chemistry (3) should achieve the no
small feat of adjusting the productivity that they yield not just to the
resolution, but to the initial conditions, too.

The diffusively–coupled Lagrangian scheme, having a diffusivity tunable
independently of the resolution, is not affected by these problems, and
reproduces fairly well with just $128^{2}$ particles the results of the
pseudo–spectral code using the same strengths of the diffusive coupler as
those used for Figure 1A.

The four panels of Figure 2 show the field $c_{2}$ at time $t=100$ as computed
by the pseudo–spectral scheme with 128 and 4096 grid points (panels A, B), and
by the Lagrangian scheme (panels C, D) with diffusivities matching those of
the pseudo–spectral calculations. The Lagrangian solutions are visualized by
plotting partially overlapping colored squares centered at the particles’
positions, rather than by resampling the solution on a regular grid. This
choice makes evident that the Lagrangian solution in panel D), reproduces the
same range of fluctuations as the solution in the panel B), even though it
obviously cannot resolve the fine structures created by the advective
dynamics.

The rest of the paper is organized as follows: in the following section we
describe the diffusive couplers; in section 3 we compare the results obtained
through our Lagrangian methods against known exact solutions or numerical
solutions obtained with a pseudo–spectral code at much higher resolution; in
section 4 we briefly discuss how to efficiently implement the methods; finally
some concluding remarks are offered in section 5.

##  2 Diffusive couplers

We are not going to attempt a discretization of the Laplacian operator:
evaluating the second derivatives of a field on a set of randomly distributed
points and then devising a numerical scheme that satisfies mass conservation
and the maximum principle would be quite challenging. Of the two methods that
we propose, the first the discrete counterpart of a convolution with the heat
kernel; the second represents diffusive processes as exchanges of mass among
nearby particles. Both methods have free parameters, which determine the
strength of the diffusive effects. More precisely, they determine the rate at
which the variance of a scalar field is dissipated. In section 3.1 we give an
objective, quantitative way to attach an effective diffusivity to a given set
of parameters.

We feel that the first coupler has a more mathematically elegant formulation.
However, it requires an iterative procedure to converge, which may make it
slow. The second coupler is little more than a recipe to destroy variance, but
its computational cost scales linearly with the number of particles.

The problem of diffusively coupling the Lagrangian particles in (2) is greatly
simplified if one takes a fractional step approach [24, §17.2]. The reaction
and advection terms are solved by integrating the ODEs (2) from time $t$ to
time $t+\tau$, then a separate diffusive step, which solves the heat equation,
is performed. Here, for notational simplicity, we illustrate the methods for
the case of one scalar field. The generalization of the methods to the $n$
scalar fields of the full PDEs (1) is straightforward.

###  2.1 First coupler

In place of a discretized form of the heat equation, we seek a discretized
form of its solution; the latter, for a scalar field $c$, is given by the
following convolution integral

$c(\boldsymbol{x},t+\tau)=\int_{\Omega}k(\boldsymbol{x},\boldsymbol{y},\tau)c(\boldsymbol{y},t)\,\mathrm{d}\boldsymbol{y}$
(6)

where the kernel $k$ is the fundamental solution of the heat equation in the
domain $\Omega$ subject to the desired boundary conditions. In
$\mathbb{R}^{d}$ the kernel is

$k(\boldsymbol{x},\boldsymbol{y},\tau)=\left(\frac{1}{4\pi
D\tau}\right)^{\frac{d}{2}}\exp\left(-\frac{\left\|\boldsymbol{x}-\boldsymbol{y}\right\|^{2}}{4D\tau}\right)$
(7)

where $D$ is the diffusion coefficient of the heat equation.

Given $M$ points $\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{M}$ in $\Omega$,
let $W_{ij;\tau}$ be the elements of a matrix representing a discrete
counterpart of the convolution (6) evaluated at the points
$\boldsymbol{x}_{i}$, $\boldsymbol{x}_{j}$ and across a time interval $\tau$.
By analogy with the properties of the kernel (7), we shall assume $W$ to be a
non–negative, symmetric matrix. The simplest discretization of the convolution
(6) is given by

$c_{i}(t+\tau)=\sum_{j=1}^{M}W_{ij;\tau}c_{j}$ (8)

where we use the shorthands $c_{i}=c(\boldsymbol{x}_{i},t)$ and
$c_{i}(t+\tau)=c(\boldsymbol{x}_{i},t+\tau)$. If

$\sum_{j=1}^{M}W_{ij;\tau}=1,$ (9)

that is, each column of $W$ sums to $1$, then the expression (8) is just a
weighted average of all the concentration values $\\{c_{i}\\}$. Therefore, it
satisfies the maximum principle in the form:

$\min_{i=1,\ldots,M}\\{c_{i}\\}\leq
c_{i}(t+\tau)\leq\max_{i=1,\ldots,M}\\{c_{i}\\}.$ (10)

If each row of $W$ sums to $1$, i.e.

$\sum_{i=1}^{M}W_{ij;\tau}=1$ (11)

then the expression (8) satisfies the conservation of mass in the form

$\sum_{i=1}^{M}c_{i}(t+\tau)=\sum_{j=1}^{M}\left(\sum_{i=1}^{M}W_{ij;\tau}\right)c_{j}=\sum_{j=1}^{M}c_{j}.$
(12)

Thus, if the discrete kernel $W$ is a _doubly–stochastic_ matrix [25], i.e. it
satisfies both (9) and (11), then the discrete model (8) obeys both the
maximum principle and the conservation of mass.

Let us now describe how to construct such a discrete kernel $W$. Initially we
define a crude discretization of the exact kernel (7) as follows

$K_{ij;\tau}=\begin{cases}\exp\left(-\frac{\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|^{2}}{4\mathcal{D}\tau}\right),&\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|<m\sqrt{2\mathcal{D}\tau}\\\
0,&\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|\geq
m\sqrt{2\mathcal{D}\tau}\end{cases}$ (13)

where the nominal diffusivity $\mathcal{D}$ must be intended as a free
parameter. The kernel $K$ has a cut–off determined by $m$, also a free
parameter, to avoid computing the negligible contribution of pairs of
particles too far away from each other. Because $K$ is not, in general, a
doubly–stochastic matrix, we need to find a doubly–stochastic surrogate of
$K$.

The problem of rescaling a given matrix into a doubly–stochastic one is named
_balancing_ , and dates back to the 1930s. Since then, a large number of
applications has been solved by resorting to the balance of matrices (see,
e.g., [26] for a rich list of examples).

We say that a matrix $K$ can be balanced if there exist two diagonal matrices,
$\text{diag}(\boldsymbol{a})$ and $\text{diag}(\boldsymbol{b})$, such that

$W=\text{diag}(\boldsymbol{a})K\,\text{diag}(\boldsymbol{b})$ (14)

is doubly–stochastic. The fundamental theorem addressing this problem for
non–negative matrices is due to Sinkhorn and Knopp [25]. Starting from any
vector $\boldsymbol{a}_{0}$ with positive elements, they propose the following
iteration:

$\boldsymbol{b}_{k+1}=\left(K^{T}\boldsymbol{a}_{k}\right)^{-1};\quad\boldsymbol{a}_{k+1}=\left(K\boldsymbol{b}_{k}\right)^{-1}$
(15)

where the reciprocal is intended to be applied element–wise. Their theorem
then states that the process converges to a doubly–stochastic matrix of the
form (14) with $\boldsymbol{a}=\lim_{k\to\infty}\boldsymbol{a}_{k}$,
$\boldsymbol{b}=\lim_{k\to\infty}\boldsymbol{b}_{k}$, if $K$ has _total
support_. A matrix $K$ is said to have total support if every positive entry
in $K$ can be permuted into a positive diagonal with a column permutation.
Under the conditions of the theorem the balancing is unique: $K$ can be turned
into one and only one doubly–stochastic matrix by means of multiplication by
diagonal matrices (which are themselves unique up to a scalar factor).

Our crude discretization of the Gaussian kernel, the matrix (13), has total
support, because it is symmetric and has a positive main diagonal. Therefore,
if $K_{ij}$ is a non–zero element, then the column permutation that swaps
column $i$ with column $j$ brings to the main diagonal $K_{ij}$, $K_{ji}$, and
no other element; the main diagonal thus remains positive. We can then define
the discrete convolution kernel $W$ that appears in (8) as the balancing of
$K$. For our purposes it is important to note that $K$ and $W$ have the same
pattern of zeros, therefore the particle pairs coupled by $W$ are all and only
those coupled by $K$.

###  2.2 Second coupler

A way to represent small-scale irreversible mixing processes is suggested by
physical intuition, along the following heuristic argument, similar to those
used in [17, 19]. When two fluid particles happen to be close enough, they
will exchange some portion of their mass, and, thus, of their advected
scalars. Let $q_{ij}\geq 0$ be the mass fraction exchanged between the $i-$th
and the $j-$th particle, which are assumed to have the same mass. This
fraction may be a function of the distance
$\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|$ and may be assumed to
be zero when the distance exceeds some fixed threshold. Thus the concentration
of the scalar $c$ after a diffusion step at the position of the $i-$th
particle will be

$c_{i}(t+\tau)=c_{i}-\sum_{j=1}^{M}q_{ij}c_{i}+\sum_{j=1}^{M}q_{ij}c_{j}$ (16)

where the first sum represents the losses to other particles, and the second
sum represents the gains from other particles. The above expression can be re-
arranged as

$c_{i}(t+\tau)=\left(1-\sum_{j=1}^{M}q_{ij}\right)c_{i}+\left(\sum_{j=1}^{M}q_{ij}\right)\overline{c}_{i}$
(17)

where the overline denotes the weighted average
$\overline{c}_{i}=\sum_{j=1}^{M}q_{ij}c_{j}/\sum_{j=1}^{M}q_{ij}$. If

$0\leq\sum_{j=1}^{M}q_{ij}\leq 1$ (18)

equation (17) shows that $c_{i}(t+\tau)$ is a linear interpolation between
$c_{i}$ and $\overline{c}_{i}$, and therefore the maximum principle is
satisfied.

In addition, it is straightforward to verify that
$\sum_{i}c_{i}(t+\tau)=\begin{subarray}{c}\sum_{i}\end{subarray}c_{i}$, and
therefore the expression (16) conserves mass.

As exchange fraction we shall use

$q_{ij}=\begin{cases}\frac{p}{\left(4\pi\mathcal{D}\tau\right)^{\frac{d}{2}}}\,\exp\left(-\frac{\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|^{2}}{4\mathcal{D}\tau}\right),&\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|<m\sqrt{2\mathcal{D}\tau}\\\
0,&\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|\geq
m\sqrt{2\mathcal{D}\tau}\end{cases}$ (19)

where $p$, $\mathcal{D}$ and $m$ are free parameters and $d$ is the
dimensionality of the space. This particular choice is loosely suggested by
the fact that if the scalar field carried by the $i-$th particle at time $t$
were represented by a delta function, a diffusion process having diffusivity
$\mathcal{D}$, after a time $\tau$ would spread out the scalar over the whole
domain with a resulting concentration proportional to
$\exp\left(-\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|^{2}/(4\mathcal{D}\tau)\right)$.
The cut–off for large distances is also physically motivated: the small-scale,
unresolved advective motions that this diffusion process is supposed to
represent, cannot occur at an arbitrarily large speed; therefore, in a finite
time $\tau$ only particles closer than some threshold length may exchange
mass.

Special care must be taken in choosing $p$ small enough as to enforce the
condition (18). A useful rule of thumb is:

$\frac{p}{(4\pi\mathcal{D}\tau)^{\nicefrac{{d}}{{2}}}}<\frac{1}{N(m\sqrt{2\mathcal{D}\tau})},$
(20)

where $N(h)$ is the average number of particles that fall into a sphere of
radius $h$.

###  2.3 Boundary conditions

So far we have discussed the diffusive couplers as if the computational domain
were unbounded. When the domain is limited, any condition enforced along its
boundaries is reflected in the kernel $k$ appearing in the convolution
solution (6), which ceases to be a simple Gaussian function.

In the case of periodic boundary conditions, the kernel is an infinite sum of
Gaussians, one for each of the periodic images. For example, on the segment
$[0,2\pi)$ the kernel is

$k(x,y,\tau)=\sum_{n\in\mathbb{Z}}\frac{1}{\sqrt{4\pi
D\tau}}\exp\left(-\frac{\left(x-y+2n\pi\right)^{2}}{4D\tau}\right).$ (21)

If $m\sqrt{2D\tau}<\pi$, and we accept to approximate to zero the exponential
when its argument is larger than or equal to $m$ (as we do in (13) and in
(19)), then only one term gives a non–zero contribution in the sum. This shows
that the expressions (13) and (19) remain valid for periodic boundary
conditions, provided that the norms
$\left\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\right\|$ which appear in those
expressions are considered as the minimum distance in the periodic domain
between the particle $i$ and the particle $j$.

Another common boundary condition prescribes that the flux of tracers across
any portion of the boundary has to be zero. When no particle is seeded outside
of the domain, this condition is automatically enforced by both the diffusive
couplers presented here. There is, however, a pitfall that needs to be brought
to light. This is most easily illustrated in a one–dimensional domain. Let us
consider the half–line $[0,\infty)$. If we impose no–flux (a.k.a Neumann)
boundary conditions at $x=0$, then the heat kernel is

$k(x,y,\tau)=\frac{1}{\sqrt{4\pi
D\tau}}\left[\exp\left(-\frac{\left(x-y\right)^{2}}{4D\tau}\right)+\exp\left(-\frac{\left(x+y\right)^{2}}{4D\tau}\right)\right].$
(22)

This can be deduced by imposing an even symmetry to the initial condition
which extends the problem to the whole line, and then restricting the solution
back to the half–line. The even symmetry enforces the boundary condition. This
implies that the points at $x>0$ do exchange fluxes across the boundary with
their mirror images at $x<0$, but do so as to keep equal to zero the net flux
at $x=0$. If these virtual fluxes across the boundary are not taken into
account, then, in proximity of the boundaries, the diffusivity of the scalar
field is underestimated, even though the no–flux boundary condition is still
correctly enforced. A solution to this problem might consist in using ghost
particles strategically placed outside the domain so as to represent an
even–symmetric field across it. In more than one dimension, this would be
relatively straightforward only for straight boundaries, and would quickly
escalate to a challenging problem for boundaries of arbitrary shape. However,
the contribution of the mirror images is important only within a distance of
$O(\sqrt{2D\tau})$ from the boundary. In high–Péclet number, under–resolved
simulations, this distance would be comparable to or smaller than the
inter–particle distance. We thus feel that attempting to fix this issue may
not be worth the effort. In the following when we mention “no–flux boundary
condition” we refer to the straightforward case in which no ghost particles
are used.

In the test cases we have not used the Dirichlet boundary condition. However
we anticipate no difficulties in implementing this condition by distributing
particles along the boundary and fixing their concentrations to a prescribed
value. The same considerations about mirror images and ghost particles,
subject to the appropriate symmetry, apply to this case as well.

##  3 Results

###  3.1 Advection and diffusion

A first test for the diffusive couplers introduced in the previous section is
to compare their performance for advection–diffusion problems in cases in
which small–scale structures are progressively formed and eventually become
under–resolved. An analytically–solvable, well–known, but non trivial test
case is the following [27]:

$\frac{\partial c}{\partial t}+y\frac{\partial c}{\partial x}=D\nabla^{2}c$
(23)

with initial condition

$c(x,y,0)=\cos(x).$ (24)

In a domain vertically unbounded and horizontally periodic with period of
$2\pi$, the problem (23,24) has the exact solution

$c(x,y,t)=e^{-D\left(t+\frac{t^{3}}{3}\right)}\cos\left(x-yt\right)$ (25)

which develops arbitrarily high wavenumbers in the $y-$direction as times
progresses due to the tipping over of the tracer streaks operated by the
shearing flow (Figure 3).

Figure 3: Numerical solution of (23,24) using the first coupler (§2.1). The
parameters of the discretized kernel (13) are $d=2$, $m=8$,
$\sqrt{2\mathcal{D}\tau}=\pi/512$, $\tau=0.1$. The second coupler, with the
parameters of Figure 4, produces visually indistinguishable results.

Multiplying (23) by $c$, averaging, and using (25) after an integration by
parts, one finds the following explicit expression for the rate of dissipation
of scalar variance

$-\frac{d}{dt}\left\langle\frac{c^{2}}{2}\right\rangle=D\left\langle\left|\nabla
c\right|^{2}\right\rangle=\frac{D}{2}\left(1+t^{2}\right)e^{-2D\left(t+\frac{t^{3}}{3}\right)}.$
(26)

Where the angular brackets denote a spatial average over one horizontal period
and an arbitrary vertical length.

In Figure 4 this expression is compared with the results obtained using the
two couplers discussed in sec. (2).

Figure 4: Rate of dissipation of scalar variance for the problem (23,24). Blue
curve: results from the numerical simulation of Figure 3. Green curve: results
using the second coupler (§2.2), with parameters $p=1.38\cdot 10^{-5}$, $m=4$,
$\sqrt{2\mathcal{D}\tau}=\pi/256$, $\tau=0.1$ in the exchange fraction (19).
Red curve: expression (26) with $D=3.23\cdot 10^{-6}$.

The numerical computations use the domain $[0,2\pi)\times[-\pi,3\pi]$,
periodic in $x$ and with no–flux boundary conditions in $y$. The number of
particles is $128\times 256$. The averages are computed in the central part of
the domain, shown in Figure 3. The left–hand side of (26) is then computed
from the particles’ concentrations. The value of the diffusivity $D$ in the
right–hand side of (26) is least–squares fitted to the numerical results. The
fit extends from the beginning of the simulation up to the time of maximum
dissipation. The value of the parameter $p$ in the second coupler is tuned in
order to match the fitted value of $D=3.23\ldots\cdot 10^{-6}$ obtained with
the first coupler with at least two significant digits.

The match with the exact dissipation rate becomes inaccurate at later times,
because when the stripes become under–resolved the tracer variance is aliased
to lower wave numbers, and thus it is not damped as quickly as it should have
been. The advantages of the Lagrangian approach, however, should become clear
by observing that with a pseudo–spectral code at a comparable resolution, the
lowest diffusivity must be $D\approx 10^{-3}$ in order to avoid significant
spurious oscillations. With that diffusivity the dissipation rate peaks at
time $t\approx 10$ instead than $t\approx 70$, by which time the streaks have
all but disappeared. With the Lagrangian approach we are able to obtain a
qualitatively correct shape of the dissipation curve corresponding to a
diffusivity that is three orders of magnitude smaller.

In fact, for each choice of the parameters, we can define the _effective
diffusivity_ of the method as the value $D$ in the right-hand side of (26)
that best fits the growing part of the numerical dissipation curve. This
value, in general, does not coincide with the _nominal diffusivity_
$\mathcal{D},$ which appears in (13) and (19) and depends on the parameters as
we shall discuss below.

Using the first coupler, in the discrete kernel (13) we set the cut–off radius
$m\sqrt{2\mathcal{D}\tau}=h$ to be
$h=\nicefrac{{\pi}}{{8}},\,\nicefrac{{\pi}}{{16}},\,\nicefrac{{\pi}}{{32}},\,\nicefrac{{\pi}}{{64}}$.
For each of these values we consider $m=3,4,6,8,12,16$. Fixing the value of
the time step (we use $\tau=0.1$) the nominal diffusivity is then determined
as

$\mathcal{D}=\frac{h^{2}}{2\tau m^{2}}.$ (27)

Figure 5 shows the effective diffusivity as a function of the nominal
diffusivity for the above values of $h$ and $m$.

Figure 5: Effective diffusivity $D$ as a function of nominal diffusivity
$\mathcal{D}$ for the first coupler (§2.1). Different symbols correspond to
different values of the cut–off radius $h$. Different points with the same
symbol correspond to different values of $m$. The nominal diffusivity is then
given by (27). The black dashed line is the identity $D=\mathcal{D}$.

Points that have the same $\nicefrac{{h}}{{m}}$ ratio yield nearly the same
effective diffusivity. In other words, for fixed $\mathcal{D}$, the effective
diffusivity is fairly insensitive to the cut-off radius $h$, even when this is
so small that only very few particles are involved: when
$h=\nicefrac{{\pi}}{{64}}$ only $\pi$ particles, on average, fall within a
disc of radius $h$.

At high nominal diffusivities, the effective diffusivity nearly coincides with
the nominal one: $D(\mathcal{D})\approx\mathcal{D}$. At low nominal
diffusivities the effective diffusivity appears to be proportional to the
square of the nominal one: $D(\mathcal{D})\propto\mathcal{D}^{2}$. Further
tests suggest that the constant of proportionality scales as the square root
of the particle density, and that the switch between the two regimes occurs
when the standard deviation $\sqrt{2\mathcal{D}\tau}$ of the discrete kernel
(13) is of the same order of magnitude as the average distance between nearest
particles. We did not further investigate the reasons of this change of slope
and postpone an in-depth examination of the issue to a further work.

Figure 6 shows

Figure 6: Effective diffusivity $D$ as a function of nominal diffusivity
$\mathcal{D}$ for the second coupler (§2.2). Symbols have the same meaning as
in Figure 5. Blue markers refer to computations with $p=10^{-4}$, red to
$p=10^{-5}$, green to $p=10^{-6}$. The black dashed lines are the functions
$D=10^{n}\mathcal{D}$, with $n=-5,-4,\ldots,0$.

the effective diffusivity obtained with the second coupler as a function of
the coupler’s parameters appearing in the exchange fraction (19). The markers
relative to $h=\nicefrac{{\pi}}{{64}}$, $m=12,16$ are absent, because with
those parameters the condition (18) does not hold: the method violates mass
conservation and blows up.

The cut-off radius is determined as specified above for the first coupler, and
the expression (27) for the nominal diffusivity still holds. As in that case,
the effective diffusivity is fairly insensitive to the cut-off radius $h$ when
the ratio $\nicefrac{{h}}{{m}}$ is kept fixed. In contrast with the first
coupler, the effective diffusivity appears to be roughly proportional to the
nominal one across the whole range of diffusivities that we have tested. The
effective diffusivity also appears to be roughly proportional to the parameter
$p$.

The effective diffusivity of the second coupler also depends on the density of
the particles. If, keeping all other parameters the same, we double the
average number of particles that fall within a disk of radius $h$, we find,
from (16) and (19), that the average mass exchanged on a time step by each
particle with its neighbors doubles. Thus the effective diffusivity is
proportional to the particle density.

###  3.2 Reaction and diffusion

The methods described in the present work are designed for cases in which the
Péclet numbers are extremely high. However, it cannot be excluded that some
geophysical flows may, occasionally, be characterized by less extreme Péclet
numbers. It is thus of interest to verify what may be the performance of the
methods when the advection terms are not dominant over the diffusion ones. In
the limit of zero Péclet numbers, the equations (1) reduce to
reaction–diffusion equations. Even though we are not proposing our methods for
this class of problems, we found informative to use one of them as a test
case.

Here we will consider the well-known Fisher–Kolmogorov–Petrovskii–Piskunov
equation, namely

$\frac{\partial c}{\partial t}=D\nabla^{2}c+c(1-c).$ (28)

For non-negative $c$, this equation describes the propagation of fronts
joining a stable ($c=1$) and an unstable ($c=0$) region (e.g. [28]§13.2).
There exist solutions with fronts propagating at any speed $V\geq 2\sqrt{D}$.
However, for a very large class of initial conditions, in particular those
whose derivative has compact support, the propagation speed is the minimal
one[29]: $V=2\sqrt{D}$.

When the function $c$ assumes negative values the solution generally blows–up
to minus infinity in a finite time. It is thus important to avoid numerical
solution methods that generate spurious oscillations. In particular, this may
be a problem when the diffusion coefficient is small, because the thickness of
the front is also proportional to $\sqrt{D}$. Thus, low diffusivities imply
high gradients in the traveling front.

We produce a numerical approximation of (28) by uniformly random seeding
$128^{2}$ particles in the square $[0,2\pi]\times[0,2\pi]$. We use no–flux
boundary condition. Initially, all particles have a concentration of zero,
except those having a coordinate $x<0.2$, whose concentration is set to one.
We then advance the solution with time steps of length $\tau=0.1$ by
alternating one of the diffusive couplers of sec. 2 and the evaluation of the
exact solution of the equation $\dot{c}=c(1-c)$.

Figure 7: Speed of propagation of the front in the solution of equation (28)
as a function of the effective diffusivity. Symbols have the same meaning as
in Figure 5. Blue markers refer to computations using the first coupler
(sec.2.1) and the the green ones to computations using the second coupler
(2.2) with $p=10^{-5}$. The black dashed line is the theoretical speed
$V=2\sqrt{D}$.

In Figure 7 we plot the propagation speed of the front as a function of the
effective diffusivity of the method, evaluated as detailed in the previous
subsection. The first coupler gives the best results, while the second coupler
overestimates the speeds by about a factor 2.5. With both couplers the front
propagation speed appears to be proportional to the square root of the
diffusivity, as in the exact solution, except at very low diffusivities, where
the front speed declines somewhat faster than the exact scaling. This
excessive slow–down is in qualitative agreement with what was found in a
stochastically forced version of equation (28). The primary effect of the
random forcing was that of damping the leading tail of the propagating front,
thus slowing it down [30]. We speculate that the random arrangement of the
particles may play the role of the stochastic forcing.

The front is well–resolved only at the lowest diffusivities. When $D\approx
10^{-3}$ the thickness of the front becomes comparable with the interparticle
distance. Thus most of the results of Figure 7 refer to runs in which the
front is poorly resolved or not resolved at all. When the front is not
resolved, the separation between the region where $c=1$ and $c=0$ appears as a
jagged line, with meanders of characteristic size determined by the
interparticle distance.

We could not run this test case with a cut–off radius $h=\pi/64$, because this
length results to be smaller than the percolation threshold: due to the random
inhomogeneities in the distribution of the particles, after a short transient,
no particle with concentration zero is found at a distance less than $h$ from
a particle with concentration higher than zero, thus the front stops
propagating. In figure 7, we used $h=\pi/48$, instead. This elucidates the
disadvantage of not having a velocity field stirring the particles: although
Poissonian random gaps in the distribution of particles exist even in the
presence of a stirring velocity field, they open and close as time progresses,
rather than remaining static, and are thus far less damaging, as the results
of the other tests should clearly illustrate.

###  3.3 Advection, reaction and diffusion at different Damköhler numbers

We now return to the simple resource–consumer model (3) to test the
performance of the Lagrangian couplers when the Damköhler number is changed.
Here we do so by letting the reaction rate assume the values
$r=0.04,\,0.2,\,1,\,5$, while keeping in all cases the same velocity field,
which is defined by the following streamfunction

$\psi(x,y,t)=\left[(n\bmod 2)\sin(x+\phi_{n})-\left(1-(n\bmod
2)\right)\sin(y+\phi_{n})\right]$ (29)

where $n=\left\lfloor t\right\rfloor$ (the largest integer smaller than $t$),
“$\bmod$” denotes the remainder of the integer division, and $\phi_{n}$ is a
uniformly random phase chosen in $[0,2\pi)$. This is an example of a “random
renewing flow” (see e.g. [31]§11.1) which is very effective at mixing an
advected scalar field. The characteristic spatial scale of the flow is
constant, but an advected field is subject to a continuous process of
stretching and folding that produces a cascade of progressively smaller
scales.

Our benchmarks are numerical solutions of the problem (1) with the chemistry
(3) and the velocity field induced by (29), solved on a uniform grid with
$4096^{2}$ nodes, on the doubly–periodic domain $[0,2\pi)\times[0,2\pi)$, with
a Fourier–Galerkin pseudo–spectral code, and a diffusion coefficient
$D=\nicefrac{{0.003}}{{32^{2}}}\approx 2.9\cdot 10^{-6}$. A slightly larger
diffusivity was used than in the computations of Figure 1 at the same
resolution, because at higher reaction rates the solution develops higher
gradients in the concentration fields. We thus have tuned $D$ so as to obtain
a solution free of spurious oscillations at $r=5$, and we have kept that value
for all the reaction rates. We use both the uniform consumer initial condition
(4) and the non overlapping blobs initial condition (5).

Figure 8: Mean (panels A,B) and standard deviation (panels C,D) of the
consumer concentration field as a function of time using the chemistry (3) and
the stirring field (29). Panels A,C refer to the initial conditions (4);
panels B,D to the initial conditions (5). Different colors denote different
reaction rates, as specified in the legend of panel D. Solid lines refer to
results obtained with a pseudo–spectral code on a grid with $4096^{2}$ nodes.
Dotted and dashed lines refer to the Lagrangian method with $128^{2}$
particles and respectively, the coupler of section 2.1 and of section 2.2.

Against the benchmark we compare the results obtained using the Lagrangian
method with the couplers of section 2. For the first coupler we use a cut–off
radius $h=\pi/64$ and $m=5.8$. For the second coupler we use $h=\pi/32$,
$m=4$, $p=10^{-5}$. In both cases $128^{2}$ particles were used, the time step
is $\tau=0.1$ and the ODEs (2) are solved with the standard fourth–order
Runge–Kutta scheme. The results are summarized in Figure 8.

Both Lagrangian methods reproduce very well the time evolution of the mean of
the chemical fields, and reasonably well their standard deviation. If measured
with the criterion of section 3.1, the parameters above yield an effective
diffusivity slightly higher ( $D\approx 1.1\cdot 10^{-5}$) than the
diffusivity of the pseudo–spectral code (the criterion suggests $m\approx 8$
for the first coupler and $m\approx 7.5$ for the second). If the effective
diffusivity matches that of the pseudo–spectral code, in the later stages of
the simulation the standard deviation remains too high and decays at a rate
clearly slower than in the pseudo–spectral benchmark. This occurs because, as
stirring cascades the chemical tracers to unresolved small scales, the
variance relative to those scales is aliased back to larger scales, where it
is damped at an incorrect, lower rate. Using an ad-hoc higher effective
diffusivity initially gives a slight underestimation of the standard deviation
and, later on, a slight overestimation, while producing what we consider to be
an acceptable approximation of a dynamics that requires a resolution 32 times
higher to be fully resolved.

##  4 Implementation details

An efficient implementation of the diffusive couplers of section 2 requires a
fast algorithm for finding all the particles falling within a distance $h$
from any given particle. This _fixed–radius near neighbors search_ is a
classic problem in computational geometry. For arbitrary distribution of
points, it can be solved by arranging the points in tree data structures such
as quad–trees or Kd–trees (see e.g. [32] chap. 5, 14). The use of trees leads
to algorithms with a computational cost of $O\left(M\log(M)\right)$, where $M$
is the number of points. When, as in our case, the particles are uniformly
distributed, it is more convenient to use a lattice and hashing method, which
has a computational cost of $O(M)$ [33].

The computational domain is overlaid with a regular lattice with square meshes
of size $h$. To each mesh is assigned a unique index. For simplicity we use
row–order indexing, although the Z–order indexing might improve cache
efficiency. The particles are kept in a list, sorted according to the index of
the mesh that contains each particle, which is easily computed from the
particle’s position. The sorting is performed by means of the _counting_
algorithm (e.g. [34]§8.2), which does not use pairwise comparisons, and has a
complexity of $O(M)$. A hash table associates each mesh index with the first
particle in the sorted list having that index. Thus, accessing all particles
within the same mesh is an $O(1)$ operation, because each mesh, on average,
contains the same number of particles, due to their uniform distribution. To
find all the particles within a distance $h$ from a given one, one needs to
compute the distance of the given particle from all the particles in the same
mesh and in some of the adjacent meshes (three in 2D or four in 3D). After
each time step, the particle list is sorted again, and the hash table is
updated. If the size of the mesh $h$ is decreased as the number of particles
$M$ increases in such a way as to maintain constant the average number of
particles in each mesh, then the fixed–radius neighbor search problem is
solved in $O(M)$ time. We did not attempt yet to produce a parallel version of
our prototype code. However we don’t expect to face unusual difficulties or
harsh performance penalties by pursuing a straightforward domain partitioning
strategy, in which each processor takes care of a contiguous block of meshes.

In the case of the coupler of section 2.2, the computation of the exchange
fraction (19) only increases the prefactor in the asymptotic scaling of the
fixed–radius neighbor search. The overall algorithm is thus $O(M)$.

In the case of the coupler of section 2.1, an analysis of the computational
cost is more complicated, because it needs to take into account the cost of
balancing the discrete kernel (13). The analysis of balancing algorithms is
still an open problem, and we settled for the venerable Sinkhorn–Knopp
algorithm only because it is extremely simple to implement. An assessment of
the performances and of the relative merits of balancing algorithms, in
particular on distributed–memory parallel architectures, is beyond the scope
of this paper, and might become the subject of a future work.

##  5 Discussion and conclusions

In this paper we have investigated the viability of Lagrangian numerical
methods to approximate the solution of advection–reaction–diffusion equations
in cases where it is impossible to resolve all the scales of motion, as is
commonplace for biogeochemical problems.

The methods consist in alternating a purely Lagrangian step that solves the
advection–reaction part of the equations with the method of characteristics,
with a diffusive step that couples the particles moving along the
characteristic lines of the advection–reaction problem. Two such couplers have
been proposed. One amounts to a discrete version of the convolution with a
Gaussian kernel, the other prescribes the exchanges between nearby particles
of small portions of the mass carried by each. In both cases the resulting
scheme conserves mass, respects the maximum principle and allows to tune the
diffusivity down to zero, where the couplers have no effect, and the method of
characteristics is recovered.

We have carried several tests comparing the methods against exact solutions of
advection–diffusion and reaction–diffusion problems, and against fully
resolved numerical solutions of advection–reaction–diffusion problems obtained
using a pseudo–spectral method run at significantly higher resolution than
that of the Lagrangian code. In all cases the results have been fairly good,
except in the case of the reaction–diffusion test, where the lack of an
advection term that stirs the particles hampers the performance of the method.
However, even in this unfavorable case, the methods are able to recover in a
roughly correct way the main features of the solution and their scaling as a
function of the diffusivity.

Of course, when it is impossible to resolve all the spatial scales present in
the solution, no method should be considered as completely satisfactory, and
it is very likely that special cases could be found where it would perform far
from well. For example, we don’t expect our Lagrangian method to perform
brilliantly in reproducing the propagation of chemical fronts stirred by
steady cellular flows. The speed of those fronts critically depends on an
accurate description of the tails of the tracer distribution in proximity of
the hyperbolic stagnation points at the cell boundaries [35]. When the spatial
structures are severely under resolved those tails are not reproducible, and
the resulting speed is then unlikely to be correct. On the other hand,
chemical fronts of that kind are not present in the oceans, and stagnation
points, albeit present, are not steady; typical ocean mixing processes involve
shearing, or stretching and folding dynamics, and in those cases our approach
seems to be satisfactory.

This paper does not suggest that our Lagrangian methods are competitive with,
or even comparable to, a fully resolved numerical solution obtained with an
Eulerian method, but rather that, by allowing to control the diffusivity
independently of the resolution, the Lagrangian methods offer, when resolution
can’t be further increased, a much better compromise than equally unresolved
Eulerian methods. In this respect, diffusive couplers like those presented
here could be seen more as a subgrid–scale parameterization of sorts, rather
than as a discretization of a diffusion operator such as the Laplacian that
appears in (1).

While we believe that the present work is a successful proof–of–concept, some
additional steps will be required in order to incorporate it into a realistic
ocean circulation model. The current prototype implementation needs to be
extended to three spatial dimensions and to distributed–memory parallel
architectures. In the present form the couplers only represent homogeneous and
isotropic diffusive processes. In ocean models, anisotropy is necessary, at
least in the vertical direction, and the possibility to allow for
spatially–dependent diffusivities is desirable. Finally, the existing Eulerian
parameterizations for the sources and sinks of tracers, due to interactions
with the bottom, with the air, and through river run–off must be adapted to
the Lagrangian framework. These goals will probably be easier to achieve by
modifying the coupler of section 2.2 where subgrid–scale fluxes are
represented explicitly and locally as exchanges of mass between particles.
They may be more demanding with the coupler of section 2.1, which requires the
balancing of a matrix, a process that involves all particles simultaneously,
even when the discretized kernel that couples them has a cut–off at a finite
distance.

### Acknowledgements

Part of this work was funded by Università del Salento through “Progetto
5xmille per la Ricerca”. We have benefited from discussions with Zouhair
Lachkar, Olivier Paulis and Marcello Vichi. We are indebted to Clare Eayrs,
Marina Levi and Shafer Smith who read and commented earlier drafts of the
paper.

## References

  * [1] Vichi, M., Pinardi, N., and S. Masina. “A generalized model of pelagic biogeochemistry for the global ocean ecosystem. Part I: theory.” Journal of Marine Systems, 64 (2007) 89–109.
  * [2] Pasquero, C. "Differential eddy diffusion of biogeochemical tracers." Geophysical research letters, 32 (2005), L17603.
  * [3] Mahadevan, A., and J. W. Campbell. "Biogeochemical patchiness at the sea surface." Geophysical Research Letters, 29 (2002), 32-1–32-4.
  * [4] Richards, K. J., and S. J. Brentnall. "The impact of diffusion and stirring on the dynamics of interacting populations." Journal of theoretical biology, 238 (2006), 340–347.
  * [5] Durran D. R. Numerical Methods for Fluid Dynamics: With Applications to Geophysics. Springer, 2nd edition, Berlin, 2010.
  * [6] Brentnall, S. J., et al. "Plankton patchiness and its effect on larger-scale productivity." Journal of Plankton Research, 25 (2003), 121-140.
  * [7] Lévy, M., A. Estublier, and G. Madec. "Choice of an advection scheme for biogeochemical models." Geophysical Research Letters, 28 (2001), 3725–3728.
  * [8] Lévy, M., et al. "Large-scale impacts of submesoscale dynamics on phytoplankton: Local and remote effects." Ocean Modelling, 43 (2012), 77–93.
  * [9] Lévy, M., and A. P. Martin. "The influence of mesoscale and submesoscale heterogeneity on ocean biogeochemical reactions." Global Biogeochemical Cycles, 27 (2013), 1139–1150.
  * [10] Martin, A. P., et al. "An observational assessment of the influence of mesoscale and submesoscale heterogeneity on ocean biogeochemical reactions." Global Biogeochemical Cycles, 29 (2015) 1421–1438.
  * [11] Hockney, R. W., and J. W. Eastwood. Computer simulation using particles. Adam Hilger, Bristol, 1988.
  * [12] Monaghan, J. J. "Smoothed particle hydrodynamics." Annual review of astronomy and astrophysics, 30 (1992), 543–574.
  * [13] Cleary, P. W., and J. J. Monaghan. "Conduction modelling using smoothed particle hydrodynamics." Journal of Computational Physics, 148 (1999), 227–264.
  * [14] Degond, P., and S. Mas-Gallic. "The weighted particle method for convection-diffusion equations. I. The case of an isotropic viscosity." Mathematics of computation, 53 (1989) 485–507.
  * [15] Degond, P., and F. J. Mustieles. "A deterministic approximation of diffusion equations using particles." SIAM Journal on Scientific and Statistical Computing, 11 (1990), 293–310.
  * [16] Haertel, P. T., and D. A. Randall. "Could a pile of slippery sacks behave like an ocean?" Monthly weather review, 130 (2002): 2975–2988.
  * [17] Haertel, P. T., L. Van Roekel, and T. G. Jensen. "Constructing an idealized model of the North Atlantic Ocean using slippery sacks." Ocean Modelling, 27 (2009), 143–159.
  * [18] Koszalka, I., et al. "Plankton cycles disguised by turbulent advection." Theoretical population biology, 72 (2007), 1–6.
  * [19] Bracco, Annalisa, Sophie Clayton, and Claudia Pasquero. "Horizontal advection, diffusion, and plankton spectra at the sea surface." Journal of Geophysical Research: Oceans, 114 (2009) C02001.
  * [20] Edouard, S., et al. "The effect of small-scale inhomogeneities on ozone depletion in the Arctic." Nature, 384 (1996), 444–447.
  * [21] Mariotti, A., et al. "The evolution of the ozone “collar” in the Antarctic lower stratosphere during early August 1994." Journal of the atmospheric sciences, 57 (2000), 402–414.
  * [22] Lin, J., et al. eds. Lagrangian modeling of the atmosphere. Geophysical Monograph Series Vol. 200. John Wiley & Sons, 2013.
  * [23] Martin, A. P., et al. "Patchy productivity in the open ocean." Global Biogeochemical Cycles 16.2 (2002).
  * [24] LeVeque, R. J. Finite volume methods for hyperbolic problems. Cambridge Texts in Applied Mathematics Vol. 31. Cambridge University Press, Cambridge (2002).
  * [25] Sinkhorn, Richard, and Paul Knopp. "Concerning nonnegative matrices and doubly stochastic matrices." Pacific Journal of Mathematics, 21 (1967) 343–348.
  * [26] Knight, P. A., and D. Ruiz. "A fast algorithm for matrix balancing." IMA Journal of Numerical Analysis, 33 (2013), 1029–1047.
  * [27] Rhines, P. B., and W. R. Young. "How rapidly is a passive scalar mixed within closed streamlines?." Journal of Fluid Mechanics, 133 (1983), 133–145.
  * [28] Murray J. D. (2007) Mathematical Biology I. An Introduction. Springer, Berlin, 3rd ed.
  * [29] van Saarloos W. “Front propagation into unstable states”. Physics Reports 386 (2003), 29–222.
  * [30] Doering C. R., Mueller C., and P. Smereka. “Interacting particles, the stochastic Fisher–Kolmogorov–Petrovsky–Piscounov equation, and duality.” Physica A 325 (2003) 243–259.
  * [31] Childress S. and A. D. Gilbert (1995) Stretch, Twist, Fold: The Fast Dynamo. Lecture Notes in Physics Monographs, Springer, Berlin.
  * [32] de Berg, M., van Krefeld M., Overmars M., and O. Schwarzkopf (2008) Computational Geometry. Algorithms and Applications. Springer, Berlin, third edition.
  * [33] Bentley J. L., Stanat D. F., and E. H. Williams Jr. “The complexity of finding fixed-radius near neighbors.” Information Processing Letters 6 (1977) 209-212.
  * [34] Cormen, T. H., Leiserson C. E., Rivest R. L., and C. Stein (2009) Introduction to Algorithms. MIT press, Cambridge.
  * [35] Tzella A. and J. Vanneste. “FKPP Fronts in Cellular Flows: The Large-Péclet Regime.” SIAM Journal on Applied Mathematics 75 (2015) 1789–1816.


# MRFI: An Open Source Multi-Resolution Fault Injection Framework for Neural
Network Processing

Haitong Huang12, Cheng Liu1, Xinghua Xue12, Ying Wang1, Huawei Li1, Xiaowei
Li1  1 Institute of Computing Technology, Chinese Academy of Sciences
2 University of Chinese Academy of Sciences  {huanghaitong21s,
<EMAIL_ADDRESS>
###### Abstract

To ensure resilient neural network processing on even unreliable hardware,
comprehensive reliability analysis against various hardware faults is
generally required before the deep neural network models are deployed, and
efficient error injection tools are highly demanded. However, most existing
fault injection tools remain rather limited to basic fault injection to
neurons and fail to provide fine-grained vulnerability analysis capability. In
addition, many of the fault injection tools still need to change the neural
network models and make the fault injection closely coupled with normal neural
network processing, which further complicates the use of the fault injection
tools and slows down the fault simulation. In this work, we propose MRFI, a
highly configurable multi-resolution fault injection tool for deep neural
networks. It enables users to modify an independent fault configuration file
rather than neural network models for the fault injection and vulnerability
analysis. Particularly, it integrates extensive fault analysis functionalities
from different perspectives and enables multi-resolution investigation of the
vulnerability of neural networks. In addition, it does not modify the major
neural network computing framework of PyTorch. Hence, it allows parallel
processing on GPUs naturally and exhibits fast fault simulation according to
our experiments.

##  I Introduction

Deep neural network models are increasingly used in safety-critical
applications like autonomous driving and avionics, but computational errors
caused by imperfect hardware such as process variations, defects and
impurities in materials, and noise can lead to unexpected inference and cause
severe consequences in these safety-critical applications. To guarantee
resilient deep learning, it is crucial to thoroughly evaluate the reliability
of neural network models before their deployment. While it is prohibitively
expensive to perform fault injection to realistic chips, simulation based
error injection tools such as PyTorchFI [1], Ares [2], TensorFI [3], TorchFI
[4], and MindFI [5] have been extensively developed. These tools typically
provide basic functionalities such as error injection and reliability analysis
to general neural network processing and significantly alleviate the design
efforts required for neural network reliability analysis in presence of
various hardware errors.

Nevertheless, it has become evident that existing fault injection tools are
insufficient to address complex neural network reliability analysis for the
following aspects. First, they generally utilize hardware fault models such as
bit flip model to simulate the influence of hardware errors on neural network
processing, but they do not model any specific hardware architectures and the
correlation between the neuron-wise error injection and realistic hardware-
based fault injection remains unclear. The accuracy of the fault simulation
needs to be calibrated such that it can be utilized as a common metric to
reliability study of neural network processing. Second, there is a lack of
standard error injection metrics and neural network evaluation metrics that
can be utilized to align the fault simulation results and distinct fault
simulation results may be obtained from these different fault injection tools,
which dramatically affects the reliability study of neural network processing.
For instance, many prior fault injection tools fail to clarify the model
quantization during reliability evaluation [6] and biased evaluation results
are obtained because of the quantization setups. Third, many of the fault
injection tools are closely coupled with the neural network frameworks such as
PyTorch. The error model and error injection may vary with the different
neural network layers [7], which also complicates the fault simulation. In
addition, fine-grained selective fault-tolerant design [8] has also been
explored and it poses great challenges to the fault simulation and evaluation
of existing fault injection tools. These unique fault simulation and
evaluation requirements further make the fault injection frameworks difficult
to scale and port to GPUs or other parallel computing systems.

To address the above problems, we present MRFI, a highly configurable open-
sourced fault injection tool for reliability investigation and evaluation of
neural network processing. On the one hand, it provides unified error models,
neural network quantization setups, and error injection mechanisms to ensure
consistent metrics of error injection and reliability evaluation. Moreover, it
is compared to typical fault injection platforms with more detailed
architectural details and calibrated for more fair comparison. On the other
hand, it provides a multi-resolution vulnerability investigation of neural
network processing from distinct angles and granularities, which can be
utilized to guide selective fault-tolerant neural network design conveniently.
In addition, MRFI utilizes a tree structure for flexible error injection
configurations and reliability evaluation of different neural network modules
without modifying the neural network processing engines, which will not affect
the underlying parallel computing engines on GPUs and can achieve the state-
of-the-art fault simulation.

The contributions of this paper can be summarized as follows:

  * •

We present an open-sourced multi-resolution fault injection framework MRFI and
have it calibrated with architecture-specific fault injection frameworks.

  * •

MRFI provides unified error injection and reliability evaluation with flexible
simulation configurations. Particularly, it allows fine-grained vulnerability
analysis from different angles and can be utilized for selective protection
against various hardware errors.

  * •

MRFI is built on top of PyTorch and does not affect the computing engine of
PyTorch. Hence, the fault injection and reliability evaluation can be scaled
to parallel computing systems like GPUs conveniently and achieve competitive
simulation speed.

TABLE I: Comparison of different fault injection frameworks including MRFI,
PyTorchFI[1], and Ares[2]

| MRFI | PFI | Ares
---|---|---|---
activaion injection
permanent fault | ✓ | ✓ | ×
transient fault | ✓ | ✓ | ✓
weight injection
permanent fault | ✓ | ✓ | ✓
transient fault | ✓ | × | ×
number of faults
fixed number fault | ✓ | ✓ | ×
random error rate | ✓ | × | ✓
fault injection on quantized data
quantized fault | ✓ | ✓ | ✓
full quantization | ✓ | × | ✓
quantization method selection | ✓ | × | ×
quantization parameter settings | ✓ | limited | limited
fault injection on floating point data
fixed value inject | ✓ | ✓ | ×
random perturbation | ✓ | × | ✓
random bit flip | ✓ | × | ×
internal observers
clean activations | ✓ | × | ×
fault impact | ✓ | × | ×
functional flexibility
custom error model | ✓ | ✓ | ×
multi-resolution evaluation | B,N,C,L1 | N,L | L
fine-grained simulation configuration
by python program | ✓ | limited | ×
by configuration file | ✓ | × | ×

  * 1

B: bit, N: neuron, C: channel, L: layer

##  II Related Work

Research on the reliability and fault tolerance of deep learning revolves
around error simulation, fault evaluation, fault-tolerant model design, and
fault-tolerant hardware design. Error simulation methods can be divided into
hardware-level simulation and model-level simulation. Hardware-level
simulation is generally more accurate but has limited generality and
performance due to the need to simulate the underlying logic of specific
hardware architectures. Notable works in this area include NVIDIA SASSIFI for
fault injection simulation on GPUs[9] and fault simulation on FPGAs[8] [10]
[11] [12]. These error simulation methods often require tens to thousands of
times the normal execution time.

On the other hand, model-level evaluation and simulation have better
performance, making them widely studied in previous works. Some notable works
in this area include [1] [3] [2] [13] [5] [4]. Some works have shown that
accurate hardware-level reliability results can be obtained through model-
level simulation. For example, in Are[2], the authors obtained consistent
measurements with silicon measurement by injecting errors at the model level.
Yi He et al. [7] mapped hardware unit error rates to model neurons to obtain
accurate model-level evaluation results. Xinghua Xue et al. [14] improved the
accuracy of computational evaluation by focusing on operations. Some works aim
to improve the performance of error simulation. For example, PyTorchFI[1]
optimized the implementation of fault injection simulation, reducing the
additional time overhead of single-point fault injection to a negligible
level. The closest state-of-the-art tools to our work are PyTorchFI and Ares,
which offer more comprehensive functionality compared to other works. Table I
shows a summary of their features.

For fault evaluation, recent works have proposed improved statistical metrics
and mathematical modeling analysis to further reduce the number of error
simulations required. In HarDNN[15], the authors explored various surrogate
test metrics to replace simple model accuracy measurements, resulting in a
10-times reduction in simulation count under the same simulation accuracy.
BinFI [16] models error propagation and introduces binary search to reduce the
number of simulations. [17] achieved a two-order-of-magnitude speedup in fault
evaluation by using statistical modeling. Previous fault injection frameworks
did not consider these optimizations, whereas our framework’s flexibility
allows researchers to conveniently perform these evaluations.

After fault evaluation, safety-critical systems often require fault-tolerant
design. Fault-tolerant design can be applied at both the software and hardware
levels, such as adjusting models for fault tolerance or introducing additional
hardware redundancy computations [18] [19] [20] [21] [22]. These solutions
include Selective Hardening [8], Algorithm-Based Fault Tolerance [23],
adjustment of quantization parameter settings [6], adjustment of activation
layer values to limit numerical magnitude [24], fault-tolerant retraining
[25], and more. Evaluating these solutions requires a more flexible fault
injection framework.

##  III MRFI Framework

Figure 1: The MRFI tool integrates multiple modules. The configuration part
can be used for simple and detailed experiments. The execution module provides
a variety of predefined functions.

###  III-A MRFI Overview

The overall structure of MRFI is illustrated in Figure 1. The framework as a
whole can be divided into two parts: the configuration part and the execution
modules, which are further divided into several functional modules.

MRFI provides two main approaches to configure different experimental
requirements: EasyConfig for simple and preliminary experiments, and a
detailed error injection configuration tree for more complex needs. For simple
experiments, EasyConfig allows users to specify the layers and fault injection
method through a config file. This configuration applies uniformly to the
specified set of layers, enabling users to conduct fault injection experiments
by loading the corresponding error injection configuration, similar to working
with a typical PyTorch model. To address complex and precise error injection
requirements without modifying the experimental model, MRFI employs a detailed
error injection configuration tree. This tree records the error injection
configurations for various modules and layers of a model. Users can easily
adjust the parameters within the configuration tree to meet their specific
needs. MRFI’s error injection module then dynamically performs error injection
based on the configured settings.

The execution modules of MRFI offer multiple commonly used injection methods
that encompass the majority of error injection requirements for experiments.
For enhanced flexibility, users can also customize the execution module
according to their specific needs.

The structured configuration approach of MRFI offers several advantages over
previous manual error injection implementations. In the traditional manual
approach, experimenters typically need to modify both the neural network model
code and the error injection parameters to accommodate different error
injection requirements. This process is hard especially for modern large-scale
models with numerous layers. Additionally, manually adjusting error injection
parameters for multiple fine-grained settings in dozens of model layers is
cumbersome, making experimental management challenging. MRFI overcomes these
issues by decoupling the error injection configuration from the model
definition and organizing it into separate modules, facilitating error
injection experiments with specific configurations.

###  III-B Integration with PyTorch

MRFI is built upon PyTorch, the widely adopted neural network training
framework. Error injection in MRFI is achieved through PyTorch hook, allowing
for seamless integration with the model. PyTorch serves as the foundation for
MRFI due to its extensive user base and well-established software ecosystem.
Additionally, we have found that PyTorch’s dynamic graph approach aligns well
with error injection scenarios that necessitate dynamic modification and
configuration. The hook function in PyTorch enables dynamic observation and
modification of intermediate data within network modules, making it highly
suitable for fine-grained error injection experiment configuration. An
alternative implementation approach involves adding or replacing layers in the
model with fault injection operations. However, this method requires prior
knowledge of the original layer implementations, which is impractical given
the diverse range of user-defined modules and layers. It is worth noting that
error injection at the model level only requires access to the model’s weights
and feature map information for modification. By leveraging hooks, rewriting
of the underlying computation implementation is not necessary. Our execution
module also mainly uses PyTorch tensors, which allows the error injection
process to be naturally used on the CPU and GPU to maximize performance.

###  III-C Major MRFI Components

As depicted in Figure 1, MRFI incorporates three key components for fault
injection: quantizer, selector, and error model. Additionally, an observer
feature is available for monitoring numerical values and their error impact
within the model.

####  III-C1 Quantizer

In MRFI, we recognize the significance of quantization methods in evaluating
fault tolerance. Quantization is commonly employed during inference for fixed-
point or integer computations, differing from the floating-point storage and
computation used during model training. Previous studies highlight the varying
fault-tolerance characteristics of different numerical representations and the
influence of quantization choices on model accuracy and fault tolerance. To
facilitate accurate fault simulation for quantized models, MRFI offers
multiple configurable quantization options. These simulated quantizations
encompass the quantization and de-quantization steps, converting tensors from
floating-point to integer format and vice versa. The quantizer serves as a
wrapper for the fault injection component, enabling the provision of integer
tensors to the selector and error model for simulating integer computation and
error occurrence on real hardware. It is worth noting that the quantizer is
optional, and users can omit it if they solely require simulations with
floating-point models.

####  III-C2 Selector

The selector determines the locations within the tensor where errors are
injected. By defining the selector, users can specify the manner in which
errors manifest. For example, a fixed position selector can be employed to
evaluate the impact of permanent errors. Conversely, a random position
selector serves as a random engine to simulate the occurrence of soft errors.
Soft errors are often associated with a bit error rate, and the random
position selector efficiently generates erroneous positions based on this
rate. Users also have the flexibility to customize filtering masks at the
selector level or implement custom selectors to simulate fine-grained
selective fault-tolerant protection.

####  III-C3 Error model

The error model dictates how errors are injected into the numerical values at
the selected error positions. The error model may rely on the original
numerical values, such as random bit flips for integers or floating-point
numbers, fixed bit flips, or the addition of perturbations following a normal
distribution. Alternatively, the error model may be independent of the
original numerical values, involving stuck-at 0 faults, faults with specific
fixed numerical values, or random uniform integers or floating-point numbers.
MRFI incorporates these common error models, while also allowing users to
define custom error models conveniently.

####  III-C4 Observer

MRFI provides internal observers to enhance the visibility of numerical values
and their error impact within the model. Observers that do not require fault
injection can be employed to monitor the distribution of internal activation
values. For example, an observer can track the maximum and minimum values to
determine the dynamic range of activations, aiding in the determination of
model quantization parameters. Another type of observer enables the comparison
between fault injection and non-fault injection scenarios. Before initiating
fault injection, this observer conducts a golden run comparison analysis. It
can utilize a counter to tally the number of numerical values affected by
fault injection or employ an RMSE observer to gauge the extent of error
propagation.

##  IV Experiment Results

In this section, we present the results of various error injection experiments
conducted using MRFI. Our focus is on demonstrating experimental approaches
not covered by previous error injection tools, highlighting the necessity of
utilizing a more fine-grained error injection framework.

For our experiments, we primarily employ VGG and ResNet models evaluated on
the ImageNet dataset. In subsection IV-B, we also utilize a smaller dataset,
cifar-10, to visualize error propagation.

###  IV-A Basic BER-Accuracy fault inject experiment

(a) VGG16 weights (b) ResNet18 weights (c) VGG16 activations (d) ResNet18
activations

Figure 2: Relation between bit error rates and classification error rate of
different data type.

We begin by testing and evaluating MRFI’s performance under simple
configurations. Specifically, we conduct random bit flipping error injection
on the quantization weight of the VGG16 model using the fixed-point
quantization settings of 2.12 and 3.13. Additionally, MRFI supports different
data types and quantization methods. For comparison, we include the original
32-bit floating point (f32) and layer-wise range quantization (Qlayer).

As illustrated in Figure 2(a)(b), our results of fixed-point quantization with
settings 2.12 and 3.13 align with those obtained from Ares[2], exhibiting a
significant increase in bit error rates at $BER=10^{-6}$. In contrast,
floating-point types exhibit lower error tolerance, while layered quantization
demonstrates superior fault tolerance, corroborating the conclusions from
SNR[6]. In Figure 2(c)(d), the injection of activation values exhibits a
similar trend, but not entirely the same.

While the above experiments could be conducted separately using different
fault injection tools in the past, the advantage of MRFI lies in its
integration of diverse functionalities within a unified framework.

###  IV-B Observing error propagation

To gain deeper insights into error propagation within the network, MRFI
provides a straightforward operation to apply error observers to all network
modules. These observers enhance transparency and offer diverse perspectives
for in-depth research on error impact.

Firstly, we employ a simple equation comparison observer to assess the number
of neurons affected by a single injection. We perform one-bit flipping on the
activation values and weights of the first convolution layer (conv1) of VGG11,
and then calculate the average number of errors across 10,000 input images.
The golden run accuracy and accuracy after activation injection both yield
$76.49\%$, while the accuracy after weight injection slightly decreases to
$76.34\%$. As shown in Figure 3, the number of affected neurons increases
rapidly along the layers, with nearly all neurons in the fully-connected
layers displaying differences from their original values. This indicates the
network’s resilience to errors, while also highlighting that the number of
affected neurons is not a good indicator of the severity of error effects.

Figure 3: Number of affected neurons of one bit flip on first layer in VGG-11.

(a) MAE (b) RMSE

Figure 4: Evaluate error propagation by better metrics. Figure 5: Visualize of
LeNet internal layers and error propagation.

Next, we evaluate error propagation using more robust metrics. Inspired by
[15][17], MRFI provides mean absolute error (MAE) and rooted mean square error
(RMSE) observers. In Figure 4, both MAE and RMSE reveal that weight injection
has a significantly greater influence than activation value injection when
performing single-point injection on the first convolution layer. This finding
aligns with the accuracy results. Additionally, although the number of
affected neurons rapidly increases during the propagation process, these error
measures remain relatively stable.

MRFI also includes observers for recording activation and error values
directly, facilitating the visualization of internal feature maps in neural
networks. Figure 5 showcases the raw neuron values input by LeNet on a single
image from the cifar-10 dataset, along with the errors caused by a single
activation injection in the first layer. We observe that errors occur in fixed
positions across all channels in the convolution layer’s feature map, while
nearly all neurons in the fully connected layers experience irregular
disturbances.

###  IV-C Quantization Settings Resilience Study

Recent studies have highlighted the close relationship between quantization
methods and neural network fault tolerance. Therefore, our framework
integrates quantization settings into its configuration. In subsection IV-A,
we compared the differences between two quantization methods using a simple
global quantization setting. With MRFI, we can now quantize each submodule of
the network separately to study these differences at a fine-grained level.

Figure 6 illustrates the differences in RMSE metrics for various error
injections on different layers and blocks of ResNet18. A larger RMSE value
indicates greater sensitivity to errors. It is evident that global
quantization settings exhibit significant variations in fault tolerance and
perform poorly at each layer, resulting in lower fault tolerance than layered
quantization. Conversely, layerwise quantization demonstrates stable error
resilience. Additionally, we observe that weights are more sensitive to errors
during global quantization, whereas weights and activation values exhibit
similar sensitivities during layerwise quantization.

The quantization parameters of the model also impact fault tolerance. In MRFI,
we employ layerwise quantization and set different quantization scaling
factors for each layer of the model. We evaluate accuracy and RMSE in ResNet
as an example. As shown in Figure 7, larger scale factors generally lead to
worse fault tolerance.

In practical deployment scenarios, models are typically quantized after
training, and the quantization parameters depend on both model accuracy and
target hardware requirements. The experimental differences outlined above
emphasize the need for a thorough assessment of the quantization methods’
impact when deploying security-critical neural network models.

Figure 6: The difference in error sensitivity of quantization methods at
different network modules.

(a) Accuracy (b) RMSE

Figure 7: Fault tolerance under different quantization scale parameter.

###  IV-D Evaluate Fault-Tolerant Resilience at Different Level

MRFI’s flexible configuration capabilities enable error injection analysis at
various levels, providing robust support for selective fault-tolerant design.
In this section, we demonstrate the differences in fault tolerance based on
bit, channel and pixel.

####  IV-D1 Bit sensitivity evaluation

We inject random errors into each bit of the float16 floating-point type and
fixed-point 3.13 quantization type of the VGG-11 model, observing the model’s
accuracy at different injection rates. Figure 8 illustrates the model accuracy
for the high-order bits, while the low-order bits demonstrate greater
insensitivity to error injection. In float16, the highest bit represents the
sign, followed by 5 exponent bits, and finally 10 mantissa bits. From Figure
8, we observe that high-order exponent bits are highly sensitive to errors,
whereas error injection into the mantissa bits does not significantly affect
the model’s accuracy. The sign bit exhibits a sensitivity level lower than
that of high exponent bits but higher than that of low exponent and mantissa
bits. For fixed-point 3.13 quantization, a regular pattern emerges in the
fault tolerance of each bit, with lower bits consistently displaying better
fault tolerance than higher bits.

####  IV-D2 Channel sensitivity evaluation

Figure 9 showcases the differences in sensitivity indicators for each channel
in the activation values of the first three layers of VGG-11. The three
convolutional layers contain 64, 128, and 256 feature map channels,
respectively. By employing layerwise quantization parameters, the differences
in fault tolerance among different layers are small. After sorting the RMSE
indicators for each channel in ascending order, we observe that the variance
in channel sensitivity within a layer surpasses the variance in sensitivity
between layers. This suggests that selective fault-tolerant design at the
channel level can yield enhanced fault tolerance outcomes.

####  IV-D3 Pixel sensitivity evaluation

It is also possible to evaluate the pixel sensitivity differences of feature
maps at the spatial level. We use MRFI to randomly inject errors into each
channel of different pixels with different activation values in the ResNet two
layers, resulting in the output RMSE as shown in Figure 10. We can see that
pixels in the center position are usually more sensitive to errors, which may
be related to the dataset features of ImageNet, as most targets are located at
the center of the image.

(a) float16 (b) fixed 3.13 quantization

Figure 8: Bit sensitivity evaluation of VGG-11 among all layers. Figure 9:
Channel sensitivity evaluation on three convolution layers of VGG-11. Figure
10: Pixel sensitivity evaluation on two convolution layers of ResNet18 using
RMSE metric.

###  IV-E Performance

In addition to accuracy, the performance of error injection and simulation in
neural networks is a crucial factor for conducting experiments. MRFI has
implemented optimizations in error injection, resulting in comparable time
costs between error injection and the golden run. Below, we will refer the
state of the art tools PyTorchFI and Ares to conduct two similar experiments
and measure the time consumption. PyTorchFI has been optimized for single
point error injection, making its overhead negligible compared to forward
propagation. However, PyTorchFI does not provide random error generation based
on injection rate, so we compared it with the random sampling method used by
Ares as the baseline.

####  IV-E1 Single point injection

Table II presents the evaluation time for single point injection errors on
GPUs and CPUs. It demonstrates that the execution time of the MRFI framework
can be disregarded relative to the model’s computation. This indicates that
although our tool is designed for complex injection configurations, it has
also achieved target performance close to PyTorchFI.

####  IV-E2 Injection according to bit error rate

For random multi-point injection, Ares conducts acceptance-rejection sampling
on each bit, while MRFI uses selectors to choose injection error positions
based on the expected number calculated using Equation 1. This reduces the
number of random number generations required. However, a discrepancy arises
when the injection error rate is low, as $n_{inject}$ is a discrete integer.
To address this issue, MRFI adopts the Poisson sampling strategy described in
Equation 2, ensuring consistent results in theory. Figure 11 compares the
performance and accuracy between direct probability sampling and optimization
methods. While using a selector can significantly reduce the cost of sampling
calculations, it may introduce jumps at low injection rates. On the other
hand, Poisson sampling maintains consistent results with the baseline sampling
methodology.

$n_{inject}=\lfloor N_{neurons}\times BER\rceil$ (1)

$N_{inject}\sim Poisson(N_{neurons}\times BER)$ (2)

TABLE II: Time cost comparison on single point injection | Acc. | GPU time | CPU time
---|---|---|---
golden run | 68.87% | 86.83±0.06s | 2603±2s
1p fault inject | 68.33% | 87.94±0.04s | 2607±3s

(a) Time cost (b) Accuracy

Figure 11: Comparison of performance and accuracy of different random sampling
methods for random multi-point injection.

##  V Conclusion

In this paper, we have introduced MRFI, an open-source multi-resolution fault
injection framework for neural network processing. MRFI enables refined error
assessments, considering the impact of model quantization and supporting fine-
grained error injection configurations. To simplify configuration, MRFI
decompose complex configuration requirements into multiple modules and offers
a user-friendly configuration interface. Our experiments have shown that MRFI
can evaluate errors from multiple aspects with practical needs.

## References

  * [1] A. Mahmoud, N. Aggarwal, A. Nobbe, J. R. S. Vicarte, S. V. Adve, C. W. Fletcher, I. Frosio, and S. K. S. Hari, “PyTorchFI: A Runtime Perturbation Tool for DNNs,” _Proceedings - 50th Annual IEEE/IFIP International Conference on Dependable Systems and Networks, DSN-W 2020_ , pp. 25–31, 2020.
  * [2] B. Reagen, U. Gupta, L. Pentecost, P. Whatmough, S. K. Lee, N. Mulholland, D. Brooks, and G. Y. Wei, “Ares: A framework for quantifying the resilience of deep neural networks,” in _Proceedings - Design Automation Conference_ , vol. Part F1377, 2018.
  * [3] Z. Chen, N. Narayanan, B. Fang, G. Li, K. Pattabiraman, and N. DeBardeleben, “Tensorfi: A flexible fault injection framework for tensorflow applications,” _Proceedings - International Symposium on Software Reliability Engineering, ISSRE_ , vol. 2020-Octob, pp. 426–435, 2020.
  * [4] B. Goldstein, S. Srinivasan, N. K. Mellempudi, D. Das, L. Santiago, V. C. Ferreira, N. Solon, S. Kundu, and F. M. G. França, “Reliability evaluation of compressed deep learning models,” in _2020 IEEE 11th Latin American Symposium on Circuits Systems (LASCAS)_ , 2020.
  * [5] Y. Zheng, Z. Feng, Z. Hu, and K. Pei, “Mindfi: A fault injection tool for reliability assessment of mindspore applicacions,” in _2021 IEEE International Symposium on Software Reliability Engineering Workshops (ISSREW)_ , 2021, pp. 235–238.
  * [6] E. Ozen and A. Orailoglu, “SNR: Squeezing Numerical Range Defuses Bit Error Vulnerability Surface in Deep Neural Networks,” _ACM Trans. Embed. Comput. Syst._ , vol. 20, no. 5s, 2021. [Online]. Available: https://doi.org/10.1145/3477007
  * [7] Y. He, P. Balaprakash, and Y. Li, “Fidelity: Efficient resilience analysis framework for deep learning accelerators,” _Proceedings of the Annual International Symposium on Microarchitecture, MICRO_ , vol. 2020-Octob, pp. 270–281, 2020.
  * [8] F. Libano, B. Wilson, J. Anderson, M. J. Wirthlin, C. Cazzaniga, C. Frost, and P. Rech, “Selective hardening for neural networks in fpgas,” _IEEE Transactions on Nuclear Science_ , vol. 66, no. 1, pp. 216–222, 2018.
  * [9] S. K. S. Hari, T. Tsai, M. Stephenson, S. W. Keckler, and J. Emer, “SASSIFI: An architecture-level fault injection tool for GPU application resilience evaluation,” _ISPASS 2017 - IEEE International Symposium on Performance Analysis of Systems and Software_ , vol. 1, no. 1, pp. 249–258, 2017\.
  * [10] G. Gambardella, J. Kappauf, M. Blott, C. Doehring, M. Kumm, P. Zipf, and K. Vissers, “Efficient error-tolerant quantized neural network accelerators,” in _2019 IEEE International Symposium on Defect and Fault Tolerance in VLSI and Nanotechnology Systems (DFT)_. IEEE, 2019, pp. 1–6.
  * [11] D. Xu, Z. Zhu, C. Liu, Y. Wang, H. Li, L. Zhang, and K.-T. Cheng, “Persistent fault analysis of neural networks on fpga-based acceleration system,” in _2020 IEEE 31st International Conference on Application-specific Systems, Architectures and Processors (ASAP)_. IEEE, 2020, pp. 85–92.
  * [12] D. Xu, Z. Zhu, C. Liu, Y. Wang, S. Zhao, L. Zhang, H. Liang, H. Li, and K.-T. Cheng, “Reliability evaluation and analysis of fpga-based neural network acceleration system,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 29, no. 3, pp. 472–484, 2021.
  * [13] G. Li, S. K. S. Hari, M. B. Sullivan, T. Tsai, K. Pattabiraman, J. S. Emer, and S. W. Keckler, “Understanding error propagation in deep learning neural network (DNN) accelerators and applications,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC 2017, Denver, CO, USA, November 12 - 17, 2017_ , B. Mohr and P. Raghavan, Eds. ACM, 2017, pp. 8:1–8:12. [Online]. Available: https://doi.org/10.1145/3126908.3126964
  * [14] X. Xue, H. Huang, C. Liu, T. Luo, L. Zhang, and Y. Wang, “Winograd convolution: A perspective from fault tolerance,” in _Proceedings of the 59th ACM/IEEE Design Automation Conference_ , 2022, pp. 853–858.
  * [15] A. Mahmoud, S. K. S. Hari, C. W. Fletcher, S. V. Adve, C. Sakr, N. R. Shanbhag, P. Molchanov, M. B. Sullivan, T. Tsai, and S. W. Keckler, “Hardnn: Feature map vulnerability evaluation in cnns,” _ArXiv_ , vol. abs/2002.09786, 2020\.
  * [16] Z. Chen, G. Li, K. Pattabiraman, and N. DeBardeleben, “Binfi: An efficient fault injector for safety-critical machine learning systems,” in _Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis_ , ser. SC ’19. New York, NY, USA: Association for Computing Machinery, 2019. [Online]. Available: https://doi.org/10.1145/3295500.3356177
  * [17] H. Huang, X. Xue, C. Liu, Y. Wang, T. Luo, L. Cheng, H. Li, and X. Li, “Statistical modeling of soft error influence on neural networks,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , pp. 1–1, 2023.
  * [18] C. Liu, Z. Gao, S. Liu, X. Ning, H. Li, and X. Li, “Special session: Fault-tolerant deep learning: A hierarchical perspective,” in _2022 IEEE 40th VLSI Test Symposium (VTS)_. IEEE, 2022, pp. 1–12.
  * [19] W. Li, X. Ning, G. Ge, X. Chen, Y. Wang, and H. Yang, “Ftt-nas: Discovering fault-tolerant neural architecture,” in _2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC)_. IEEE, 2020, pp. 211–216.
  * [20] X. Ning, G. Ge, W. Li, Z. Zhu, Y. Zheng, X. Chen, Z. Gao, Y. Wang, and H. Yang, “Ftt-nas: Discovering fault-tolerant convolutional neural architecture,” _ACM Transactions on Design Automation of Electronic Systems (TODAES)_ , vol. 26, no. 6, pp. 1–24, 2021.
  * [21] J. J. Zhang, T. Gu, K. Basu, and S. Garg, “Analyzing and mitigating the impact of permanent faults on a systolic array based neural network accelerator,” in _2018 IEEE 36th VLSI Test Symposium (VTS)_. IEEE, 2018, pp. 1–6.
  * [22] B. Reagen, P. Whatmough, R. Adolf, S. Rama, H. Lee, S. K. Lee, J. M. Hernández-Lobato, G.-Y. Wei, and D. Brooks, “Minerva: Enabling low-power, highly-accurate deep neural network accelerators,” in _2016 ACM/IEEE 43rd Annual International Symposium on Computer Architecture (ISCA)_. IEEE, 2016, pp. 267–278.
  * [23] S. K. S. Hari, M. B. Sullivan, T. Tsai, and S. W. Keckler, “Making convolutions resilient via algorithm-based error detection techniques,” _IEEE Transactions on Dependable and Secure Computing_ , vol. 19, no. 4, pp. 2546–2558, 2022.
  * [24] B. Ghavami, M. Sadati, Z. Fang, and L. Shannon, “Fitact: Error resilient deep neural networks via fine-grained post-trainable activation functions,” in _2022 Design, Automation & Test in Europe Conference & Exhibition (DATE)_, 2022, pp. 1239–1244.
  * [25] D. Xu, M. He, C. Liu, Y. Wang, L. Cheng, H. Li, X. Li, and K.-T. Cheng, “R2f: A remote retraining framework for aiot processors with computing errors,” _IEEE Transactions on Very Large Scale Integration (VLSI) Systems_ , vol. 29, no. 11, pp. 1955–1966, 2021.