pdf short=PDF, long=Probability Density Function, iid short=i.i.d.,
long=Independent Identically Distributed ou short=OU, long=Ornstein-Uhlenbeck,
gktl short=GKTL, long=Giardina-Kurchan-Tailleur-Lecomte, ams short=AMS,
long=Adaptive Multilevel Splitting, tams short=TAMS, long=Trajectory Adaptive
Multilevel Splitting, scgf short=SCGF, long=Scaled Cumulant Generating
Function, lbm short=LBM, long=Lattice Boltzmann Method, lbe short=LBE,
long=Lattice Boltzmann Equation, lgca short=LGCA, long=Lattice Gas Cellular
Automata, lbgk short=LBGK, long=Lattice Bhatnagar-Gross-Krook OU short=O-U,
long=Ornstein–Ulhenbeck dns short=DNS, long=Direct Numerical Simulation, md
short=MD, long=Molecular Dynamics, cfd short=CFD, long=Computational Fluid
Dynamics, pde short=PDE, long=Partial Differential Equation,

#  Numerical study of extreme mechanical force exerted by a turbulent flow on
a bluff body by direct and rare-event sampling techniques

Thibault Lestang12<EMAIL_ADDRESS>Freddy Bouchet1  Emmanuel
Lévêque2  1Univ Lyon, ENS de Lyon, Univ Claude Bernard de Lyon, CNRS,
Laboratoire de Physique, F-69342 Lyon, France 2Univ Lyon, Ecole Centrale de
Lyon, Univ Claude Bernard de Lyon, INSA de Lyon, CNRS, Laboratoire de
Mécanique des Fluides et d’Acoustique, F-69134 Ecully cedex, France

###### Abstract

This study investigates, by means of numerical simulations, extreme mechanical
force exerted by a turbulent flow impinging on a bluff body, and examines the
relevance of two distinct rare-event algorithms to efficiently sample these
events. The drag experienced by a square obstacle placed in a turbulent
channel flow (in two dimensions) is taken as a representative case study.
Direct sampling shows that extreme fluctuations are closely related to the
presence of a strong vortex blocked in the near wake of the obstacle. This
vortex is responsible for a significant pressure drop between the forebody and
the base of the obstacle, thus yielding a very high value of the drag. Two
algorithms are then considered to speed up the sampling of such flow scenarii,
namely the ams and the gktl algorithms. The general idea behind these
algorithms is to replace a long simulation by a set of much shorter ones,
running in parallel, with dynamics that are replicated or pruned, according to
some specific rules designed to sample large-amplitude events more frequently.
These algorithms have been shown to be relevant for a wide range of problems
in statistical physics, computer science, biochemistry. The present study is
the first application to a fluid-structure interaction problem. Practical
evidence is given that the fast sweeping time of turbulent fluid structures
past the obstacle has a strong influence on the efficiency of the rare-event
algorithm. While the ams algorithm does not yield significant run-time savings
as compared to direct sampling, the gktl algorithm appears to be effective to
sample very efficiently extreme fluctuations of the time-averaged drag and
estimate related statistics such as return times. Software used for
simulations and data processing is available at
https://github.com/tlestang/paper_extreme_drag_fluctuations.

##  1 Introduction

Turbulent flows are important in a variety of natural phenomena, industrial
and civil applications. Their characteristic feature is the spontaneous
development of intense and sporadic motions associated with extreme internal
forces (Lesieur, 2011; Donzis & Sreenivasan, 2010; Yeung et al., 2015).
“Extreme” refers here to fluctuations that can deviate from the mean value by
${\cal{O}}(10)$ standard deviations. In engineering, the nature of such
extreme dynamical events and their statistics are of crucial interest to
predict excessive force, which can threaten the structural integrity of
embedded structures (Kanev & van Engelen, 2010).

From the viewpoint of chaotic dynamical systems, turbulence in fluids is
linked to non-linearity and strong departure from statistical equilibrium
(Kraichnan & Chen, 1989). The use of analytical perturbative methods in
identifying resonant interactions (among degrees of freedom) responsible for
extreme fluctuations is unsuccessful. Alternatively, simulation offers a
practical approach to gain physical insight into these events, quantifying
their intensity and estimating their frequency of occurrence. However, this
requires very long simulations since these extreme events are rare. The
computational cost for sampling a fluctuation of very small probability
typically grows as the inverse of this probability (Wouters & Bouchet, 2016).
Rare-event sampling refers to a large body of methods that aim at
_preferentially_ exploring regions of phase space corresponding to events that
would otherwise be accessed with a very low probability through a brute-force
direct sampling. In the present work, a computational study of extreme
mechanical force acting on an immersed bluff body is conducted using both very
long time-series (direct sampling) and rare-event sampling techniques.

In fluid turbulence, rare-event sampling has been approached mainly from the
perspective of simplified dynamics such as the one-dimensional Burgers’
equation with a stochastic forcing (Bec & Khanin, 2007). In this case,
dynamics can be sampled by using a Markov chain Monte-Carlo algorithm (Düben
et al., 2008; Mesterházy & Jansen, 2011; Mesterházy et al., 2013) that
provides a framework for rare-event sampling. An alternative approach is based
on instantons (Gurarie & Migdal, 1996; Grafke et al., 2015) and applies to
stochastically driven systems in the limit of weak noise. Instantons refer to
the most probable trajectories in phase space that achieve a given rare event
(in the limit of weak noise). Suitable numerical schemes can be used to
evaluate instantons as well as the related probabilities of rare events
(Chernykh & Stepanov, 2001; Grafke et al., 2013; Grigorio et al., 2017; Laurie
& Bouchet, 2015; Bouchet et al., 2014). An example is the investigation of the
physics of rogue waves (Dematteis et al., 2018, 2019). However, a drawback of
the aforementioned approaches is their limitation to simple and stochastically
driven dynamics. Here, a more general approach is considered for complex,
possibly deterministic, dynamical systems. It is based on sampling algorithms
relying on _selection rules_ applied to an ensemble of trajectories of the
system. Even though such ideas date back to the early 1950s, they have
received growing interest over the last twenty years with successful
applications in various domains such as chemistry (van Erp & Bolhuis, 2005;
Escobedo et al., 2009; Teo et al., 2016), biophysics (Huber & Kim, 1996;
Zuckerman & Chong, 2017; Bolhuis, 2005), nuclear physics (Louvin et al.,
2017), nonlinear dynamical systems (Tailleur & Kurchan, 2007) and
communication networks simulation (Villen-Altamirano & Villen-Altamirano,
1994). More importantly, such algorithms have been shown to be useful for the
study of rare events in simple deterministic dynamics (Wouters & Bouchet,
2016). Certainly, an original contribution of our study is to test the
application of rare-event sampling algorithms in the context of far-from-
equilibrium dynamics with an irreducible very large number of degrees of
freedom. Two algorithms that are a priori suitable for such dynamics are
considered, namely, the Adaptive Multilevel Splitting (AMS) algorithm (Cérou &
Guyader, 2007) and the Giardina-Kurchan-Tailleur-Lecomte (GKTL) algorithm
(Giardinà et al., 2006). Another problem related to rare events is the
estimation of parameters or statistics from a limited number of samples.
Several authors have proposed original strategies (Mohamad & Sapsis, 2018;
Blonigan et al., 2019) such as sequential sampling. These approaches and rare
event algorithms could also be considered for applications in wave-structure
problems, ship motion and load acting on offshore platforms (Belenky et al.,
2016). The paper is organized in two parts. The first part highlights the
phenomenology of extreme fluctuations of the drag force acting on a square
placed in a two-dimensional turbulent channel flow. This study is based on the
simulation of the flow over a very long duration, made possible by the
relative simplicity of the flow. The motivation for this study is twofold.
Firstly, it provides a detailed description of the statistics and dynamics
related to extreme drag fluctuations. This analysis is informative from the
viewpoint of fluid mechanics and, to the best of our knowledge, has never been
reported before. Secondly, it yields reference results that are required to
validate the outputs of rare-event algorithms and to evaluate the possible
computational gain obtained from them. This assessment is developed in the
second part of the paper.

The flow set-up is introduced and the dynamics related to typical drag
fluctuations is described in section 2. The statistical properties of the drag
are then discussed. In section 3, the phenomenology of extreme drag
fluctuations is investigated based on direct sampling. Both the instantaneous
drag and time-averaged drag are considered. It is found in particular that
sampled extreme fluctuations of the instantaneous drag result from very
similar dynamics. Section 4 examines the applicability of both the ams and
gktl algorithms to the simulation of extreme drag fluctuations in the same
flow configuration. In subsection 4.1, we show that the use of the ams
algorithm is not successful, or at least not straightforwardly. This
difficulty is put in perspective with the phenomenology of extreme drag
fluctuations developed in section 3. Subsection 4.2 presents the computation
of extremes of the time-averaged drag by using the gktl algorithm. This latter
allows for an exceptional reduction of the computational cost needed to
simulate trajectories corresponding to extreme values of the time-averaged
drag. As a specific successful application, the gktl algorithm is used to
compute the return times of extreme fluctuations of the time-averaged drag
acting on the immersed obstacle. Perspectives and conclusion end this work.

##  2 Description of the numerical case study

Figure 1: Our case study is a grid-generated turbulent flow impinging onto a
fixed squared obstacle (of size $R$) located at the centre of a channel in two
dimensions. The flow is artificially damped near the end of the channel. In
the developed flow, turbulent eddies have typically the size of the square,
which results in strong load fluctuations. The vorticity is displayed with an
arbitrary colour map from blue (negative values) to red (positive values).

The drag exerted by a grid-generated turbulent flow onto a fixed squared
obstacle is considered as a representative case study (see Fig. 1). Although
real-world applications would eventually imply three-dimensional dynamics, a
simplified two-dimensional setting has been chosen here to reduce the
computational cost and allow for a systematic study. We believe that this
system embeds the characteristic features that make this study both relevant
and challenging for fluid-structure-interaction problems.

Turbulent eddies generated in the near-wake of the grid are carried
downstream. They interact with each other and grow in size as expected for
two-dimensional turbulent dynamics. The dimension of the grid is such that the
size of the eddies that hit the square is comparable to its size, resulting in
strong fluctuations of the drag acting on the square. The flow dynamics is
integrated numerically by the lbm (Krüger et al., 2017). While traditional
methods in computational fluid dynamics rely on a discretization of the
Navier-Stokes equations, the lbm considers the fluid at an underlying
mesoscopic level. Capturing the dynamics of collections of fluid particles
moving and colliding on a lattice is here preferred to solving nonlinear pde.
Further details about the lbm are given in Appendix A and references therein.
In our context, this numerical method has been chosen principally for its
computational efficiency.

The simulated flow develops in a long plane channel of dimension $513\times
129$ mesh points. The square obstacle has size $R=16$ (in mesh units) and is
located at the centre of the channel. The spacing and bar height of the
entrance grid are both equal to $R/2$ (see Fig. 1). No-slip boundary
conditions are enforced on top and bottom walls of the channel and on the
surface of the obstacle by using an halfway bounce-back procedure (Sukop & Jr,
2006). Upstream of the grid, a constant parabolic velocity profile and a
constant mass density (equal to unity) are imposed as an inlet condition. The
centreline velocity is $0.05$ in lattice units, i.e. normalised by $\Delta x$
and $\Delta t$ referring to the lattice resolution and the time-step
respectively. The initial distributions are imposed at equilibrium (see
Appendix A). In the bulk, the viscosity is adjusted so that grid turbulence is
generated with Reynolds number $\mathrm{Re_{grid}}=1200$. The reference Mach
number is equal to $0.06$ in agreement with the assumption of weak
compressibility of the lbm. Near the end of the channel, the flow is
progressively damped within a sponge zone where the viscosity is artificially
enhanced. Finally, the outlet boundary condition relies on a second-order
extrapolation of the velocity and mass density. The extrapolated distributions
are evaluated through a regularization procedure relying on a finite
difference estimation of the local stress tensor, as introduced in (Latt et
al., 2008). For more details about the implementation of the flow, see the
software repository attached to this paper
(https://github.com/tlestang/paper_extreme_drag_fluctuations).

###  2.1 The drag force

Figure 2: Snapshots of the vorticity related to typical drag fluctuations
(within one standard deviation) over a time interval of length $\tau_{c}\simeq
4T_{0}$; $\tau_{c}$ will later be identified as the correlation time of the
drag signal. The vorticity is given in lattice units.

The incoming turbulent flow exerts a fluctuating force onto the squared
obstacle. The drag is defined as the resulting force in the streamwise
$x$-direction. Formally

$f_{d}(t)=\int_{\mathcal{S}}\boldsymbol{\tau}_{x\beta}(\mathbf{x},t)\leavevmode\nobreak\
\mathrm{d}{\mathcal{S}}_{\beta}(\mathbf{x}),$ (1)

where $\mathcal{S}$ is the surface of the obstacle and $\boldsymbol{\tau}$
denotes the stress tensor (see Appendix A). Here, the viscous stress makes a
negligible contribution to the drag. The latter therefore results mostly from
pressure forces. Since the pressure on the top and bottom sides of the square
applies in the normal direction, they do not contribute to the drag. As a
consequence, the drag can eventually be expressed as the difference

$f_{d}(t)=p_{fb}(t)-p_{base}(t)$ (2)

between the pressure integrated over the upstream side of the obstacle or
forebody, $p_{fb}(t)$, and the downstream side or base $p_{base}(t)$. Pressure
fluctuations are related to the dynamics of the vorticity field. Regions of
strong vorticity correspond to strong local pressure gradients, _e.g._ as
demonstrated analytically with a Rankine vortex.

The typical timescale (turnover time) of drag fluctuations can be estimated
from dimensional analysis as

$T_{0}=\frac{R}{U},$ (3)

where $R$ is the size of the square and $U$ is the averaged velocity in the
channel. Fig. 2 displays the evolution of the vorticity field around the
obstacle for typical drag fluctuations. Vorticity is generated along the
forebody and eventually carried away by the flow. Typical fluctuations of the
drag (within one standard deviation) do not result from some preferred
arrangement of the vorticity around the obstacle.

###  2.2 The drag as a random process

Figure 3:  Temporal evolution of the drag (in lattice units) acting on the
square under the action of the impinging turbulent flow. The time has been
normalised by the turnover time related to the mean-flow velocity and the size
of the obstacle, _i.e._ $T_{0}=R/U$. The return time $r(a)$ is the averaged
waiting time between the occurrence of peak fluctuations of amplitude larger
than $a$. If the threshold $a$ is sufficiently high, one observes that $r(a)$
is much larger than the correlation time of the signal. The selected peak
fluctuations are therefore well separated.

In Fig. 3, the time signal of the drag acting on the square, $f_{d}(t)$,
appears unpredictable in details and exhibits repeated bursts of high
amplitude that deviate significantly from the averaged value. Therefore, it is
natural to model the drag as a (scalar) random process.

Drag fluctuations have been sampled along a simulation of duration
$T_{tot}=4\times 10^{6}\leavevmode\nobreak\ T_{0}$. This long simulation will
be referred to as the control run in the following. It has been made possible
by the relative simplicity of the investigated flow and the computational
efficiency of the lbm. The pdf of drag fluctuations is shown in Fig. 4a. It
deviates from a normal law and shows an exponential tail for large positive
fluctuations, i.e. ${\mathbb{P}}(f_{d})\propto e^{-\ell f_{d}}$. Fig. 4a also
displays the pdf of drag fluctuations acting on a control surface
corresponding to the periphery of the obstacle but in the absence of the
obstacle. In that case, the pdf is quasi-symmetric and does not display
exponential tails. This shows that the asymmetry of the pdf and the
development of a positive exponential tail are closely related to the no-slip
condition on the obstacle boundary.

(a) pdf of normalised drag fluctuations

(b) Autocorrelation of drag fluctuations

Figure 4: (a) pdf of normalised drag fluctuations
($\tilde{f}_{d}=f_{d}-\bar{f_{d}}$) where $\mu=\bar{f_{d}}$ denotes the time-
averaged value and $\sigma=\sqrt{\overline{{\tilde{f}_{d}}^{2}}}$ is the
standard deviation. The drag is evaluated both in the presence (blue) and in
the absence (red) of the obstacle. (b) Autocorrelation function of the drag
defined as
$C(\tau)=\overline{\tilde{f}_{d}(t+\tau)\tilde{f}_{d}(t)}\leavevmode\nobreak\
/\leavevmode\nobreak\ \overline{{\tilde{f}_{d}}^{2}}$. The correlation time
$\tau_{c}\simeq 4T_{0}$ is given by $C(\tau_{c})=0$.

Lastly, the autocorrelation function of the drag $C(\tau)$ is shown in Fig.
4b. It is found that drag fluctuations are correlated over a time interval
$\tau_{c}\simeq 4T_{0}$, illustrating that the drag loses its memory over a
time scale corresponding to the sweeping of a few eddies past the obstacle.
This observation is important for the application of rare-event algorithms as
it will be discussed in section 4. In the following, $\tau_{c}$ will be
referred to as the correlation time of the drag process. The ratio
$T_{0}/\tau_{c}$ may be viewed as a Strouhal number. The value $St=0.25$ is
consistent with common observations for flows past blunt structures at
comparable Reynolds numbers (Rodi, 1998).

##  3 Extreme fluctuations of the drag by means of direct sampling

Figure 5: Amplitude of drag fluctuations as a function of the corresponding
return time. $\tilde{f}_{d}$ denotes the drag with zero mean, i.e.
$\tilde{f}_{d}=f_{d}-\overline{f_{d}}$.

The phenomenology of extreme fluctuations of the drag is first investigated
through brute-force direct sampling applied to the control run. Direct
sampling is here used as opposed to approaches involving rare-events
algorithms discussed in section 4. It will provide a trustworthy baseline for
the validation of rare-events algorithms.

The waiting times $\tau$ are defined as the time between two consecutive
occurrences of peak fluctuations with amplitude $f_{d}\geq a$, as illustrated
in Fig. 3. The mixing time $\tau_{m}$ is the time needed for the dynamics to
lose the memory of its initial condition. As soon as the typical waiting times
are much larger than the mixing time $\tau_{m}$, the occurrences of such
events follow a Poisson process and the distribution of the waiting times is
exponential, _i.e._ $P(\tau)=\lambda(a)\exp(-\lambda(a)\tau)$ where
$r(a)=1/\lambda(a)$ is the averaged waiting time (Lestang et al., 2018);
$r(a)$ is called the return time of the level $a$. For systems without multi-
stability, it is common for the mixing time $\tau_{m}$ to be of the order of
the correlation time $\tau_{c}$.

How rare is a fluctuation $a$ is quantified by the return time $r(a)$. We can
define extreme drag fluctuations as rare events in the sense that the return
time is much larger than the correlation time, _i.e._ $r(a)\gg\tau_{c}$.

If one assumes that $r(a)=t(a)\leavevmode\nobreak\ /\leavevmode\nobreak\
\mathbb{P}(f_{d}\geq a)$ where the time scale $t(a)$ is of order $\tau_{c}$
and varies much more slowly with $a$ than ${\mathbb{P}(f_{d}\geq a)}$, one
might expect that

$r(a)\underset{a\to\infty}{\propto}\exp(\ell a)$ (4)

where $l$ is the rate describing the positive tail of the pdf of the drag
(shown in Fig. 4). Fig. 5 shows the evolution of the return time $r(a)$ with
the amplitude of fluctuation $a$, computed from direct sampling of the drag
signal $f_{d}(t)$ (Lestang et al., 2018). Consistently, it is found that the
return time $r(a)$ is well approximated by an exponential for large levels
$a$. Let us also point out some deviation from the exponential law at the
largest levels, which is probably the consequence of under-sampling.

###  3.1 Extracting extreme drag fluctuations from a very long timeseries

We have extracted the fluctuations of the drag with a return time $r(a)$
greater than $10^{4}\tau_{c}$ from the control time-series
$\\{f_{d}(t)\\}_{0\leq t\leq T_{tot}}$. This set will be considered as
representative of _extreme events_ in the upcoming study. The choice of this
particular threshold has been driven by the need to collect enough events with
large amplitude and to possibly identify generic features. According to Fig.
5, the related amplitude $a$ is found equal to $7.6\leavevmode\nobreak\
\sigma$ with $\sigma=\sqrt{\overline{\tilde{f_{d}}^{2}}}$ being the standard
deviation of the drag process, where $\tilde{f_{d}}$ denotes the drag with
zero-mean, i.e. $\tilde{f}_{d}=f_{d}-\overline{f_{d}}$. Precisely, 104
independent fluctuations with $\tilde{f}_{d}(t)\geq 7.6\sigma$ have been
identified. Each fluctuation is characterized by its maximal value,
$f_{d}^{\star}$, and the time, $t^{\star}$, at which this maximum is reached.
In the following, the phenomenology of extreme drag fluctuations will be
examined on the basis of this set of events.

###  3.2 Extreme fluctuations of the instantaneous drag

####  3.2.1 Contribution of forebody and base pressure fluctuations to the
overall drag fluctuation

Figure 6:  Relative contributions of the forebody and base pressure variations
to extreme amplitudes of the drag fluctuation $\tilde{f}_{d}$. An extreme
event corresponds to an amplitude $\tilde{f}^{\star}_{d}$ and a unique pair
($\tilde{p}^{\star}_{base}$, $\tilde{p}^{\star}_{fb}$).

In section 2, it was pointed out that drag fluctuations within one standard
deviation were not associated with any particular arrangement of the vorticity
around the obstacle. We shall see in the following that the situation is
different in the case of extreme drag fluctuations. Let
$(t^{\star},f_{d}^{\star})$ refer to an extreme-drag event. The (zero-mean)
fluctuation $\tilde{f}_{d}^{\star}=f_{d}^{\star}-\overline{f_{d}}$ can be
decomposed into

$\tilde{f}_{d}^{\star}=\Delta p_{fb}^{\star}-\Delta p_{base}^{\star}$ (5)

where $\Delta p_{fb}^{\star}$ and $\Delta p_{base}^{\star}$ denote the
variations of the forebody and base pressure, respectively. Fig. 6 displays
the relative contributions $\Delta p_{fb}^{\star}/\tilde{f}_{d}^{\star}$ and
$-\Delta p_{base}^{\star}/\tilde{f}_{d}^{\star}$ to the drag fluctuation
$\tilde{f}_{d}^{\star}$. It is found that very large fluctuations of the drag
result mainly ($\sim 80\%$) from a drop in base pressure, whereas the
variation of forebody pressure contributes much less to
$\tilde{f}_{d}^{\star}$. On the contrary, moderate fluctuations arise from
combined variations of the forebody and base pressures without any particular
predominance, which is in agreement with the previous observations (see
Fig.2).

####  3.2.2 Fluid dynamics related to extreme drag fluctuations

Figure 7:  Ensemble average of drag signals centred around extreme
fluctuations occurring at $t=t^{\star}$. The blue line shows the mean profile
whereas the shaded area indicates variations (around the mean profile) within
one standard deviation. Extreme-drag events exhibit a typical lifetime of one
correlation time $\tau_{c}$. The profile is slightly skewed indicating that
the step up is slower than the return to typical values.

The focus is now on the flow scenarii that yield extreme values of the drag.
Fig. 7 displays the mean profile (in time) of the drag signal around extreme
events. A peaked profile is observed with a width roughly corresponding to one
correlation time $\tau_{c}$. This shows that the duration of extreme events
corresponds typically to the sweeping time of the flow past the obstacle.
Interestingly, the profile is also slightly skewed, indicating that the step
up of the drag is slower than the return to typical values past the peak
value. This asymmetry (under time reversal) is closely linked to the symmetry
breaking in what happens before and after the obstacle. To better understand
the flow scenarii leading to these events, the vorticity fields around the
obstacle are now examined.

Figure 8:  Vorticity field (in lattice units) around the obstacle at
$t=t^{\star}$ for the highest drag amplitudes recorded in the control run.
Figure 9:  Pressure field (in lattice units) and velocity streamlines at
$t=t^{\star}$. In the near wake of the square, a (blocking) vortex blocks an
intense vortex against the base of the obstacle. Figure 10:  Snapshots of the
vorticity field (in lattice units) around $t=t^{\star}$.

Fig. 8 displays the vorticity field (in lattice units) around the obstacle for
the highest amplitudes of the drag during the control run. In each case, an
intense vortical structure is visible near the base of the obstacle. The
vorticity level of this structure is typically twice the amplitude of typical
vorticity fluctuations observed in Fig. 2. The formation of this vortex
originates from an intense negative (or positive) vorticity layer at the top
(or bottom) boundary of the obstacle. This high vorticity is responsible for a
significant pressure drop at the base of the obstacle and therefore a strong
drag. In contrast, nothing special happens near the forebody of the obstacle
during extreme-drag events.

The high pressure drop near the base of the obstacle appears to be closely
related to the presence of a strong vortex blocked against the base. As
illustrated in Fig. 9, this blockage is enforced by the presence of opposite
vorticity in the near wake, which holds the vortex against the base of the
obstacle and prevents it from being swept away for a while. This scenario is
better evidenced by Fig. 10, where the time history of the vorticity field
around $t=t^{\star}$ for the same event is shown. Before the occurrence of the
extreme event, positive vorticity originating from the bottom boundary layer
develops in the near wake of the square. This positive vorticity prevents the
shedding of negative vorticity and enforces the development of an intense
vortex against the base of the square. As the blocking vortex is in turn
advected downstream, the vortex against the base is released. Consistently,
one can argue that the typical duration of this scenario is related to the
sweeping time of the flow past the obstacle, and is therefore of the order of
$\tau_{c}$. This is in full agreement with the typical duration obtained from
statistical consideration on the mean profile of large-drag fluctuations in
Fig. 7. This scenario is generic and has been observed for most extreme events
sampled in the control run.

Figure 11:  Evolution of the (integrated) shear along the top or bottom sides
of the obstacle as a function of the drag for $t^{\star}-2\tau_{c}\leq t\leq
t^{\star}+2\tau_{c}$. Each trajectory corresponds to a single event. The blue
line is the mean path averaged over the set of extreme events sampled in the
control run.

Since the occurrence of large drag amplitudes arises from the production of
vorticity along the top or bottom side of the square, it is proposed to
characterize the dynamics of extreme events by their trajectory in the
parameter space $(f_{d}(t),\bar{\gamma}(t))$ where $\bar{\gamma}(t)$ is the
averaged shear along the top or bottom boundary of the square:

$\overline{\gamma}=\frac{1}{R}\int_{\mathcal{S}_{\parallel}}\frac{\partial
u(\mathbf{x})}{\partial y}\mathrm{d}\mathbf{x},$ (6)

where $R$ denotes the size of the square, $u$ is the streamwise component of
the velocity field and $\mathcal{S}_{\parallel}$ is the surface of either the
top or the bottom boundary. Fig. 11 shows $\overline{\gamma}(t)$ as a function
of the instantaneous drag $f_{d}$(t) for $t^{\star}-2\tau_{c}\leq t\leq
t^{\star}+2\tau_{c}$ for the 104 sampled extreme events. Before and after the
extremal fluctuation, i.e. for $t^{\star}-2\tau_{c}\leq t\leq
t^{\star}-\tau_{c}$ and $t^{\star}+\tau_{c}\leq t\leq t^{\star}+2\tau_{c}$,
paths wander in a region related to typical values of both $\overline{\gamma}$
and $f_{d}$. On the contrary, the drag abruptly varies for
$t^{\star}-\tau_{c}\leq t\leq t^{\star}+\tau_{c}$ near the extremal amplitude.
These excursions always go clockwise, that is, $\overline{\gamma}$ attains its
maximum value before $f_{d}$ does. This is consistent with an increase of
$\overline{\gamma}$ acting as a precursor for extreme drag fluctuations. In
this representation, we also observe that the path related to the increase of
the drag is longer than the path related to the return to typical values,
which may be related to the asymmetry of the mean profile displayed in Fig. 7.

###  3.3 Extreme fluctuations of the time-averaged drag

Figure 12: Instantaneous drag signals $f_{d}(t)$ corresponding to the highest
fluctuations of the averaged drag $F_{T}$ with a time window $T=10\tau_{c}$;
$\sigma$ and $\sigma_{T}$ denote the standard deviations of the instantaneous
and the averaged drag, respectively. Figure 13: pdf of the time-averaged drag
$F_{T}$ with $T=10\tau_{c}$.

We discussed previously the phenomenology of extreme fluctuations of the
_instantaneous_ drag, and identified the sweeping time of the flow past the
obstacle as the characteristic lifetime of these events. In applications, this
duration may be much smaller than the response time of the material structure
subject to these fluctuations, justifying a practical interest in the averaged
(in time) drag force. Therefore, another relevant observable is the time-
averaged drag defined as

$F_{T}(t)=\frac{1}{T}\int_{t}^{t+T}f_{d}(t)\mathrm{d}t,$ (7)

where $f_{d}(t)$ denotes the instantaneous drag and $T$ is the investigated
timescale (response time). In the following, we shall consider $T=10\tau_{c}$,
where $\tau_{c}$ is the correlation time of the instantaneous drag. In that
case, the PDF of $F_{T}$ is found nearly Gaussian as a consequence of the
Central Limit theorem (see Fig. 13).

During a time interval $[t;t+T]$, a fluctuation of $F_{T}(t)$ may be roughly
viewed as the overall contribution of $T/\tau_{c}$ independent fluctuations of
the instantaneous drag $f_{d}$. It is thus legitimate to ask whether a large
value of the averaged drag results from a single outstanding fluctuation of
the instantaneous drag (case (1)), or from an unusual succession of moderate
positive fluctuations (case (2)). In the same way as in section 3.1, one can
identify extreme fluctuations of $F_{T}$ exceeding some fixed threshold $a$,
and sample a set of extreme events. Setting $a=5.2\sigma_{T}$ with
$\sigma_{T}$ being the standard deviation of $F_{T}$, $84$ independent extreme
events were sampled. Again, this choice results from a compromise between the
need to consider large deviations from the mean value and the requirement to
sample a sufficient number of events for meaningful statistics. As a rule of
thumb, the threshold has been set so that about one hundred events are
sampled.

Fig. 12 displays the time-series $\\{f_{d}(t)\\}_{t^{\star}\leq t\leq
t^{\star}+T}$ for several extreme fluctuations of $F_{T}$ occurring at
$t=t^{\star}$. Qualitatively, it is found that extreme fluctuations of the
time-averaged drag can neither be reduced to case (1) nor case (2). Indeed,
both cases are featured in Fig. 12; very large value of the averaged-drag
appear to result from either a very large fluctuation, or a significant
succession of moderate (positive) fluctuations of the instantaneous drag.

##  4 Rare-event sampling algorithms

In the limit of very rare events and complex dynamics such as turbulent flows,
the computational cost of direct sampling becomes prohibitive. Indeed, the
return time of fluctuations $f_{d}\geq a$ scales as $r(a)\propto e^{\ell
a}\leavevmode\nobreak\ (\ell>0)$ according to Eq. (4). As a consequence, the
computational cost required to sample events of amplitude $f_{d}\geq a$
through brute-force sampling diverges typically as $e^{\ell a}$. The aim of
rare-event sampling algorithms is therefore to sample extreme fluctuations for
a computational cost much lower than their return time. Rare-event (sampling)
algorithms compute an ensemble of $N$ trajectories
$\\{\mathbf{x}_{n}(t)\\}_{0\leq t\leq T_{a}}$ with $n=1\cdots N$, where
$\\{\mathbf{x}_{n}(t)\\}_{0\leq t\leq T_{a}}$ refers to a trajectory of
duration $T_{a}$ in the phase space of the system. At each step of the
algorithm, some trajectories are discarded and others are replicated in order
to preferably sample trajectories with extreme fluctuations. The algorithm
tracks the ratio of the probability of the new ensemble with the probability
of the previous one, allowing the estimation of the statistical bias and
therefore the inference of the statistics of extreme events.

Different algorithms use different _selection rules_. The success of the
algorithm essentially depends on the quality of the selection rule to detect
the precursors of the extreme event (Rolland & Simonnet, 2015). However, with
the exception of rare analytical cases, the optimal selection rules are not a
priori known (Lestang et al., 2018). In the following we consider the ams
(ams) and gktl (gktl) algorithms. They both proved to be efficient for various
dynamics but adopt opposite strategies: the ams algorithm uses as a selection
rule the value of a predetermined score function, whereas the gktl uses a
selection rule based on the increment of a score function. A complete
description of these two algorithms and their operating principles are
provided in appendices B and C.

The ams algorithm (Cérou & Guyader, 2007) builds on previous ideas about
splitting approaches (Kahn & Harris, 1951; Glasserman et al., 1999). In recent
years, it has allowed for the computation of rare events in problems involving
a large number of degrees of freedom such as molecular dynamics simulations
(Aristoff et al., 2015; Teo et al., 2016) or stochastic partial differential
equations for the computation of rare trajectories in the Allen-Cahn equations
(Rolland et al., 2016). More recently it has been applied to rare events in
stochastic models of wall-turbulence (Rolland, 2018) and atmospheric dynamics
(Bouchet et al., 2019). A review of the ams algorithm, its history and
applications is also available in (Cérou et al., 2019).

During the last decade, the main theoretical framework for the study of rare
events in statistical physics has been the theory of large deviations
(Touchette, 2009). Alongside, numerical methods have been developed to sample
rare events (Del Moral, 2004). Among them the gktl algorithm (Giardinà et al.,
2006) is particularly suitable for estimating the probability of observables
which are temporal integrals over a very long period (Giardinà et al., 2011;
Laffargue et al., 2013). Recently the gktl algorithm has proven to be
extremely efficient to simulate extreme heat waves in a comprehensive climate
model (Ragone et al., 2018), with a gain of a factor 100 to 1000 for the
computation time. This achievement already represents a significant leap in
the applicability of rare-event sampling to complex dynamical systems. The aim
of this part of our work is to test rare-event algorithms for fluid-structure
interaction in a turbulent flow, which is an unexplored area of application.

###  4.1 Extreme instantaneous drag forces with the ams algorithm

The ams algorithm is here used with trajectories of fixed duration $T_{a}$
(Lestang et al., 2018). At each iteration, trajectories with the lowest maxima
of the score function $\xi(\mathbf{x}(t),t)$ during the time interval $0\leq
t\leq T_{a}$ are discarded. These trajectories are re-sampled by branching
them from the remaining trajectories. The operating principle of the ams
algorithm is detailed and sketched in Appendix B. The objective is here to
sample flow evolutions which exhibit extreme fluctuations of the drag
$f_{d}(t)$ acting on the square obstacle. The observable itself is used as the
score function, _i.e._ $\xi(\mathbf{x}(t),t)=f_{d}(t)$. Since the
Navier–Stokes dynamics is deterministic, small random perturbations are
artificially introduced at branching points; this procedure is detailed in
Appendix D. The chaotic dynamics then ensures that re-sampled trajectories
separate from their parents over a time interval $\tau_{L}$, usually referred
to as the Lyapunov timescale. Based on linear response theory, we expect this
small perturbation to have a negligible impact on the statistics of the
sampled rare events. This has been tested by performing a long simulation of
the dynamics and by regularly perturbing the flow. We checked that the
obtained statistics of the drag were consistent, within an error of the same
amplitude as the perturbation amplitude, with the statistics computed from the
(unperturbed) control simulation (see Appendix D).

In addition to the score function, two important parameters are the number of
trajectories $N$ and their duration $T_{a}$. The size of the ensemble $N$
governs the statistical error affecting quantities averaged over the sampled
set of trajectories. Therefore, $N$ should be taken as large as possible to
reduce these errors with, nevertheless, a practical limit given by the
available computational resources. In this work, we have performed two
numerical simulations with $N=32$ and $N=256$. The duration $T_{a}$ should be
much larger than $\tau_{c}$ (Lestang et al., 2018) but, again, kept small
enough to limit the computational cost. In practice, the duration of the
trajectories was eventually set to $T_{a}=5\tau_{c}$ in both simulations; we
checked in particular that larger values of $T_{a}$ did not improve the
results.

Figure 14:  Ensemble of $N=32$ trajectories after $181$ iterations of the ams
algorithm. In this experiment, the algorithm is used with the instantaneous
drag $f_{d}(t)$ as score function. Each trajectory has a duration
$T_{a}=5\tau_{c}$ where $\tau_{c}$ is the correlation time of the
instantaneous drag. Figure 15:  Maxima of the instantaneous drag throughout
$256$ re-sampled trajectories (vertical axis) as a function of the
corresponding computational cost $C_{AMS}$ in correlation-time unit. The ams
algorithm fails to efficiently sample rare trajectories associated to drag
fluctuations higher than the largest fluctuation already captured in the
initial ensemble.

Fig. 14 displays the ensemble of drag values after many iterations of the ams
algorithm with $N=32$. Interestingly, one obtains that all trajectories are
eventually re-sampled from the same trajectory displaying the highest maximum
in the initial ensemble, and overlap over most of their duration. The
algorithm thus fails to enhance the sampling of extreme events. This result is
also confirmed by increasing the number of initial trajectories to $N=256$.
The maximum drag achieved by the re-sampled trajectories is displayed as a
function of the computational cost in Fig. 15; the distribution of the maximal
drag for the initial trajectories is also shown. We observe that the
trajectories with the lowest maxima of the score function are discarded after
a few iterations, and new trajectories with higher maxima are re-sampled.
However, the re-sampled trajectories never exceed the amplitude of the highest
maximum already attained in the initial set of trajectories. A
phenomenological explanation is developed in the next paragraph.

It takes a time $\tau_{L}$ (Lyapunov timescale) before a re-sampled trajectory
separates from its parent. In our situation, this “memory effect” is related
to the fact that the score function is of dimension much smaller (one) than
the dimension of the phase space or, in other words, that the score function
results from the contribution of a very large number of degrees of freedom. As
shown in section 3.2.2, extreme drag fluctuations have a lifetime $\tau_{c}$
related to the timescale over which a vortex remains trapped against the base
of the obstacle. After $\tau_{c}$, the vortex is swept away by the flow and
further large fluctuations of the drag can only result from the trapping of
new vortices. Fig. 14 shows that $\tau_{c}$ is shorter than the Lyapunov’s
timescale $\tau_{L}$. Therefore, the re-sampling of a trajectory branched
close to $t=t^{\star}$ (when the maximum drag occurs) cannot lead to larger
values at $t^{\star}\leq t\leq t^{\star}+\tau_{L}$. For
$t-t^{\star}>\tau_{L}$, the drag process has lost the memory of the drag
fluctuations on which the re-sampling was based and, thus, the probability of
observing a new extreme fluctuation is also very low.

The difference between the typical duration of drag fluctuations $\tau_{c}$
and the Lyapunov timescale $\tau_{L}$ may be heuristically associated with the
so-called _turbulence rate_ (Frisch, 1995). As discussed previously, the
duration of extreme fluctuations of the drag is closely related to the
sweeping time of the flow past the obstacle, and consequently to the mean-flow
velocity $U$. On the contrary, the Lyapunov timescale is rather associated
with the intrinsic evolution of turbulent fluctuations in the reference frame
of the mean flow, i.e. with the root mean square velocity $u_{rms}$. The ratio
$u_{rms}/U$ (turbulence rate) is much lower than one in our case of grid-
generated turbulence, which thus justifies that $\tau_{L}>\tau_{c}$.

In summary, a straightforward application of the ams algorithm with the score
function being the drag itself does not allow us to efficiently sample extreme
fluctuations. This behaviour is independent of the choice of $N$ and $T_{a}$.
Increasing the size of the initial ensemble, or the duration of each
trajectory, can only increase the amplitude of the global maximum reached
initially but does not solve the issue of overlapping trajectories.

###  4.2 Extreme time-averaged drag forces with the gktl algorithm

The sampling of extreme fluctuations of the time-averaged drag $F_{T}$ is now
examined. The ams algorithm could be used in the same way as before by taking
the time-averaged observable itself as the score function. However, this would
lead to similar unsatisfactory results. For a time-averaged observable, an
alternative approach is provided by the gktl algorithm (Giardinà et al., 2006;
Tailleur & Kurchan, 2007; Giardinà et al., 2011). Similar to the ams
algorithm, the gktl algorithm relies on the simulation of an ensemble of
trajectories. At regular time intervals, some elements of the ensemble are
killed and others are cloned according to a weight that depends on the history
of the element itself. The weights are chosen so that, after several
iterations of the algorithm, the trajectories in the ensemble are distributed
according to a biased probability distribution that favors trajectories
related to large values of the time average of the observable. The gktl
algorithm belongs to a family of algorithms known as “go-with-the-winners”
(Aldous & Vazirani, 1994; Grassberger, 2002). Similar ideas have already been
applied in a wide range of fields such as polymer physics (Grassberger et al.,
1998), out of equilibrium statistical physics (Nemoto et al., 2017), computer
science (Aldous & Vazirani, 1994), dynamical systems (Tailleur & Kurchan,
2007), quantum mechanics (Kosztin et al., 1996). The application of a go-with-
the-winners approach to the computation of large deviations in non-equilibrium
systems has first been proposed in (Giardinà et al., 2006). Over the last ten
years, it has been successfully applied to investigate rare events in both
stochastic (Giardinà et al., 2006; Lecomte & Tailleur, 2007; Garrahan et al.,
2007) and deterministic systems (Giardinà et al., 2006; Tailleur & Kurchan,
2007). The operating principle of the gktl algorithm is developed in Appendix
C.

Figure 16:  Rare-event sampling of the (zero-mean) time-averaged drag
$\tilde{F}_{T}=F_{T}-\overline{F_{T}}$ with $T=10\tau_{c}$; $\tau_{c}$ is the
correlation time of the instantaneous drag. The shaded pdfs are estimated from
the biased ensemble resulting from the gktl algorithm applied to $N=1024$
trajectories of duration $T_{a}=T$ with the cloning period $\tau=\tau_{c}/2$.
The dashed line refers to the unbiased pdf of $\tilde{F}_{T}$, i.e. obtained
from direct sampling or with $k=0$ (no bias) in the gktl algorithm.

The application of the gktl algorithm is considered for the flow dynamics
introduced in section 2. The purpose is to speed-up the sampling of
trajectories with extreme fluctuations of the time-averaged drag, $F_{T}$. Our
observable of interest is therefore the drag $f_{d}(t)$ and the duration
$T_{a}$ of each trajectory corresponds to the period of averaging, _i.e._
$T_{a}=T$. In a nutshell, trajectories are first evolved in time independently
up to $t=\tau$, with $\tau<T_{a}$ referring to a cloning period. Then, the
selection (and cloning) rules apply according to the average of the observable
of interest, here $f_{d}(t)$, over the interval $[0,\tau]$. This procedure is
repeated $n-1$ times over the intervals $[\tau,2\tau],\leavevmode\nobreak\
...,\leavevmode\nobreak\ [(n-1)\tau,n\tau=T_{a}]$. As a result, the resampled
trajectories are distributed according to a probability distribution that is
tilted towards large values of the averaged observable. Further details are
provided in Appendix C.

The algorithm depends on three parameters: the cloning strength $k$, the
number of trajectories $N$ and the resampling period $\tau$. The higher $k$,
the larger is the bias involved in the statistical resampling. Similar to the
ams algorithm, $N$ governs the error affecting the averaged quantities
evaluated from the biased ensemble of trajectories, and should be chosen as
high as possible. The cloning period $\tau$ determines how often the
resampling is performed. A small cloning period can result in a loss of
information if the clones do not separate from their parents between two
cloning steps. On the contrary, choosing $\tau$ much higher than $\tau_{c}$
results in insufficient cloning steps. As a result, a rule of thumb is
$\tau\approx\tau_{c}$.

In the following experiments, $N=1024$ and $T_{a}=10\tau_{c}$, which yields a
computational cost $C_{gktl}=N\times T_{a}\approx 10^{4}\tau_{c}$. We fixed
$\tau=\tau_{c}/2$ as a satisfactory compromise in order to ensure an efficient
sampling (Lestang, 2018). Three numerical experiments corresponding to three
different values of the strength parameter $k$ have been carried out. Fig. 16
shows the histograms of the (zero-mean) time-averaged drag for the three
ensembles of trajectories, in addition to the unbiased histogram based on a
gaussian approximation of the pdf of the time-averaged drag. As expected, the
algorithm samples preferentially trajectories with a higher value of the
averaged drag. Furthermore, the higher $k$, the stronger the bias.
Nevertheless, it should be mentioned that for a given number $N$ of
trajectories, there is necessarily an upper limit $k_{max}$ of the strength
parameter over which the finite number of trajectories becomes detrimental to
the efficiency and accuracy of the selection procedure. For $k\gtrsim
k_{max}$, the re-sampling relies only on a small number of “independent
trajectories” and most of the trajectories in the biased ensemble overlap.
This effect is highlighted in Fig. 16 where the histogram corresponding to
$k=0.03$ becomes artificially peaked. In the present simulations with $N=1024$
trajectories, one can empirically estimate that $k$ should be kept smaller
than $k_{max}\approx 0.03$ to ensure the independence of the trajectories in
the biased ensemble. Fig. 17 displays the drag signal for the extreme
trajectories sampled by the gktl algorithm.

Figure 17: Drag timeseries corresponding to the highest fluctuations of the
averaged drag in the ensemble of trajectories sampled by the gktl algorithm.
The algorithm was applied with $N=1024$, $\tau=\tau_{c}/2$ and $k=0.03$.

In addition, Fig. 18 displays the vorticity field related to the maximum of
the drag in sampled trajectories. One observes that the extreme events are
consistent with the picture of strong vorticity being trapped in the vicinity
of the base of the obstacle, as pointed out in section 3.2. Qualitatively, the
related drag signals are also found very similar to those obtained from direct
sampling. However, it should be stressed that the computational cost for
sampling the events shown in Fig. 17 and Fig. 18 is roughly one hundred time
lower than the computational cost required to capture by direct sampling the
events displayed in Fig. 8 and Fig. 12.

Figure 18: Vorticity field associated with the maximum of the instantaneous
drag $f_{d}$ in trajectories sampled by the gktl algorithm.

####  4.2.1 Computation of return times

Importantly, the gktl algorithm provides an ensemble of trajectories over
which statistics of rare events can be evaluated. In this section, we show in
particular that this ensemble allows for the estimate of return times for
fluctuation amplitudes that would be unreachable by direct sampling (with the
same computational cost). The method for estimating return times is described
in Appendix C.2.

Fig. 19 displays the return times for extreme fluctuations of the time-
averaged drag obtained by using the gktl algorithm with two different values
of the strength parameter $k$. Both estimates have been obtained with the same
computational cost $N\times T_{a}$. An estimate given by direct sampling with
a time-series of the same effective duration $T_{tot}=N\times T_{a}$ is also
shown. Whilst a direct approach cannot access (by definition) events with a
return time greater than $T_{tot}$, the gktl algorithm allows us to estimate
the statistics of time-averaged drag fluctuations having a return time several
orders of magnitude above $T_{tot}$. Alternatively, for a fixed target return
time, the use of the gktl algorithm can reduce the computational cost of the
estimation by several orders of magnitude. This is obviously a major advantage
of this rare-event sampling algorithm.

Figure 19:  Return times for the time-averaged drag acting on the square
obstacle. $\tilde{F}_{T}$ denotes the time-averaged drag with zero mean. The
blue and red lines are obtained from the biased ensemble of trajectories
generated by the gktl algorithm with $N=1024$, $T_{a}=30\tau_{c}$ and a
cloning period $\tau=\tau_{c}/2$. The green line is the return times obtained
from a single timeseries of duration equal to the computational cost of both
gktl experiments. Uncertainty ranges for the gktl estimates are computed as
the standard deviation over a set of 10 independent experiments. Uncertainty
ranges for the direct estimation are computed as the standard deviation over a
ensemble of direct estimates resulting from 60 independent timeseries.

##  5 Conclusion

The dynamics and statistics of extreme fluctuations of the drag acting on a
squared obstacle mounted in a turbulent channel flow (in two dimensions) have
been investigated numerically through both direct and rare-event sampling
methods. By means of direct sampling based on a very long simulation, it was
first observed that such extreme events are generically related to the
trapping of a strong vortex against the base of the obstacle by a streamline
cap. This arrangement does not persist over time, however, since the main flow
eventually sweeps away the surrounding fluid structures responsible for this
trapping. Therefore, the lifetime of extreme drag fluctuations is found of the
order of the sweeping time of the flow past the obstacle, and the
corresponding drag signal is very peaked around extreme fluctuations. In
addition, it was found that extreme fluctuations of the time-averaged drag do
not preferentially result from a small number of very large fluctuations or an
exceptional succession of moderate fluctuations that pile up to yield a large
value of the average; both configurations are observed. A second part of this
study has been dedicated to the application of two representative rare-event
algorithms, namely the ams and the gktl algorithms, to our fluid-mechanical
problem. These algorithms rely on selection rules that determine how an
initial ensemble of trajectories is evolved to possibly enhance the sampling
of extreme fluctuations. On the one hand, the ams algorithm fails (for this
specific application) to generate trajectories exhibiting extreme events at a
better rate than a direct sampling. This result can be related to the
phenomenology of extreme drag fluctuations, whose lifetime is shorter than the
timescale over which the replicated trajectories manage to separate from their
duplicate. Therefore, the algorithm is unable to benefit from the precursors
of extreme fluctuations to enhance the realization of new rare-event
trajectories. On the other hand, the gktl algorithm leads to a computational
gain of several orders of magnitude. The latter is based on the cumulative
evolution of the system rather than its instantaneous behaviour, as opposed to
the ams algorithm. The downside is that the algorithm provides only statistics
of extreme values of the time-averaged fluctuations. This can nevertheless be
of great interest for computing return times of extreme (time-averaged)
fluctuations, for instance.

Selection rules at the heart of rare-event algorithms rely heavily on the
choice of a score function that drives the selection and replication of
trajectories. In this study, the observable itself (the drag) has been chosen
as the score function. Optimizing the choice of this score function is
desirable (especially for the ams algorithm) but difficult in practice, since
it should account for the phenomenology of the extreme events themselves, for
example how these events build up. Our study shows this quite clearly. A
possible direction would then be to take advantage of recent advances in
_learning methods_ to dynamically optimise the score function.

##  6 Acknowledgements

The authors thank Francesco Ragone, Corentin Herbert, Charles-Edouard Bréhier
and Eric Simonnet for useful discussions and suggestions on various aspects of
this work. T.L and F.B acknowledge support from the European Research Council
under the European Union’s seventh Framework Program (FP7/2007-2013 Grant
Agreement No. 616811). Simulations have been performed on the local HPC
facilities at École Normale Supérieure de Lyon (PSMN) and École Centrale de
Lyon (PMCS2I). The facilities at PSMN are supported by the Auvergne-Rhône-
Alpes region (GRANT CPRT07-13 CIRA) and the national Equip@Meso grant
(ANR-10-EQPX-29-01).

##  Appendix A The lbm

In the LB method, the fluid is viewed as a population of particles that
collide, redistribute and propagate along the different links of a discrete
lattice (see (Krüger et al., 2017) for a comprehensive introduction). In our
two-dimensional situation, the so-called D2Q9 lattice with only nine possible
velocities $\\{\mathbf{c_{i}}\\}_{i=0...8}$ at each node has been adopted (see
Fig. 20). Locally, the macroscopic flow variables (per unit volume) are
recovered by summing over the densities of particles $\\{f_{i}\\}_{i=0...8}$
moving with the different velocities, i.e.

$\rho(\mathbf{x},t)=\sum_{i}f_{i}(\mathbf{x},t)\quad\mathrm{and}\quad\rho(\mathbf{x},t)\mathbf{u}(\mathbf{x},t)=\sum_{i}f_{i}(\mathbf{x},t)\mathbf{c_{i}}$

for the mass density and the fluid momentum respectively. The assumption of
weak compressibility (for an ideal gas) is made so that the pressure is
directly proportional to the mass density: $p=c_{s}^{2}\rho$ where $c_{s}$ is
interpreted as a speed of sound.

The complexity of the flow emerges from the repeated application of simple
rules of streaming and collision. The lbm advances the local densities of
particles $f_{i}(\mathbf{x},t)$ moving with velocities $\mathbf{c}_{i}$ in a
two-step procedure. Namely, an _exact_ streaming step

$f_{i}(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta
t)=f_{i}^{\mathrm{out}}(\mathbf{x},t)$

during which particles move with their own velocity to a neighbouring node,
and an instantaneous collision step

$f_{i}^{\mathrm{out}}(\mathbf{x},t)=-\frac{1}{\tau_{\nu}}\left(f_{i}(\mathbf{x},t)-f_{i}^{\mathrm{eq}}(\mathbf{x},t)\right)$

which achieves a relaxation of local densities towards an absolute equilibrium
(at the macroscopic level). The time-scale $\tau_{\nu}$ (in lattice unit) is
related to the kinematic viscosity of the fluid by

$\nu=\left({\tau_{\nu}}-\frac{1}{2}\right)c_{s}^{2}\leavevmode\nobreak\ \Delta
t$

This simplification of the collision kernel is known as the BGK approximation
in the kinetic theory of gas (Bhatnagar et al., 1954). The equilibrium
function is given by

$f_{i}^{\mathrm{eq}}(\mathbf{x},t)=w_{i}\rho(\mathbf{x},t)\left(1+\frac{\mathrm{u}(\mathbf{x},t)\cdot\mathbf{c_{i}}}{c_{s}^{2}}+\frac{u_{\alpha}(\mathbf{x},t)u_{\beta}(\mathbf{x},t)({c_{i}}_{\alpha}{c_{i}}_{\beta}-c_{s}^{2}\delta_{\alpha\beta})}{2c_{s}^{4}}\right)$

with the weight factors $w_{0}=4/9,\leavevmode\nobreak\ w_{1...4}=1/9$ and
$w_{5...8}=1/36$ for the D2Q9 lattice. This discrete Lattice Boltzmann scheme
is second-order accurate in $\Delta x$ and compliant to the weakly-
compressible Navier-Stokes equations with a third-order error in
$\mathrm{Ma}=|\mathbf{u}|/c_{s}$ as the lattice spacing vanishes, i.e. $\Delta
x\to 0$ (Succi, 2018).

Figure 20: Sketch of the D2Q9 lattice. Particles move exactly from a lattice
node towards one of its nine neighbours (including the node itself) during one
time step. By definition, the lattice spacing is related to the time step by
$\Delta x/\Delta t=\sqrt{3}c_{s}$ where $c_{s}$ is interpreted as a speed of
sound.

As mentioned before, the pressure is directly accessible from the mass
density: $p=\rho c_{s}^{2}$. The viscous stress is also obtained easily from
the densities of particles by

$\tau^{\mathrm{visc.}}_{\alpha\beta}=-\frac{\nu}{\tau_{\nu}\leavevmode\nobreak\
c_{s}^{2}\Delta
t}\sum_{i}{c_{i}}_{\alpha}{c_{i}}_{\beta}(f_{i}-f_{i}^{\mathrm{eq}})$

so that the total stress expresses as

$\tau_{\alpha\beta}=-c_{s}^{2}\sum_{i}f_{i}\leavevmode\nobreak\
\delta_{\alpha\beta}-\frac{\nu}{\tau_{\nu}\leavevmode\nobreak\ c_{s}^{2}\Delta
t}\sum_{i}{c_{i}}_{\alpha}{c_{i}}_{\beta}(f_{i}-f_{i}^{\mathrm{eq}})$ (8)

Finally, let us mention that in the present context of turbulent flows, the
single-relaxation-time BGK collision has been replaced by a multi-relaxation-
time procedure based on central moments with an improved stability (Rosis,
2016).

##  Appendix B The ams algorithm

###  B.1 Operating principle

Figure 21:  Operating principle of the ams algorithm with three trajectories
of duration $T_{a}$. $\xi(x(t))$ is a prescribed score function defined along
each trajectory. Selection and branching rules are applied in order to reach a
maximum score value within the time interval $[0,T_{a}]$.

The operating principe of the AMS algorithm is sketched in Fig. 21. The
trajectories 1, 2 and 3 represent the evolution of the _score function_ for
the current ensemble of trajectories. On the basis of their respective
maximum: $Q_{1}$, $Q_{2}$ and $Q_{3}$, the trajectory with the lowest maximum
is discarded in the ensemble (dashed line). Among the two remaining
trajectories, trajectory 3 is chosen randomly and copied until it reaches the
value $Q_{1}$. It is then simulated from this branching point to the final
time $T_{a}$. In the present case of deterministic dynamics, a small
perturbation is introduced at the branching to separate the trajectories. This
re-sampling procedure can be iterated $J$ times or until all trajectories do
exceed a fixed threshold $Q$.

###  B.2 Application to a simple case: the ou process

Figure 22: Efficiency of the ams algorithm with respect to direct sampling in
the case of an Ornstein-Uhlenbeck process (Lestang et al., 2018). The red line
represents the evolution of the maximum obtained from re-sampled trajectories
as a function of the computational cost $C_{AMS}$. The blue line is the
analytical solution for the return time of amplitude $a$.

A one-dimensional ou process is considered:

$\dot{x}=-x+\eta(t),$ (9)

where $\eta$ is a Gaussian noise with
$\langle\eta(t)\eta(t-t^{\prime})\rangle=\delta(t-t^{\prime})$.

The ams is applied to a set of $N=32$ trajectories $\\{x_{n}(t)\\}_{0\leq
t\leq T_{a}}$ with $T_{a}=5\tau_{c}$. Let us note that the correlation time is
$\tau_{c}=1$ for the process defined by Eq. (9). Our objective is to sample
fluctuations $x\geq a$ with $a$ being very large compared to the typical
values of $x$. The score function is simply $x(t)$ and only one trajectory is
re-sampled at each iteration. The computational cost of the algorithm after
$J$ iterations is therefore related to the simulation of the $N$ initial
trajectories and the number of re-sampled trajectories. Fig. 22 compares the
computational cost of the ams algorithm with that of a direct sampling. In the
latter, the typical computational cost is simply the return time $r(a)$. One
can see that the successive re-samplings of the ams algorithm lead rapidly to
trajectories exhibiting extreme fluctuations. For large $a$, the computational
cost is many orders of magnitude lower than that obtained by direct sampling.

The ou process showcases the efficiency of the ams algorithm. However, the
state space is here one-dimensional and the choice of the score function is
straightforward: it is $x$ itself. In addition, the noise term in Eq. (9) has
no correlation in time, which implies that newly generated trajectories
quickly separate from their parents. Such favourable features do not persist
in the case of fluid dynamics.

##  Appendix C The gktl algorithm

###  C.1 The operating principle

The gktl algorithm is based on the simulation of an ensemble of $N$
trajectories $\left\\{\mathbf{x}_{n}(t)\right\\}_{0\leq t\leq T_{a}}$ with
$n=1\cdots N$ starting from independent random initial conditions. Let us
consider a real-valued observable of interest $A(\mathbf{x}(t))$, _e.g._ the
drag $f_{d}(t)$, and introduce a cloning period $\tau$. At time instants
$t_{i}=i\tau$ with $i=1,\leavevmode\nobreak\ 2,\leavevmode\nobreak\
...,\leavevmode\nobreak\ T_{a}/\tau$ ($T_{a}$ is a multiple of $\tau$) a
weight $W_{n}^{i}$ is assigned to each trajectory. This weight is defined
($t_{0}=0$) by

$W_{n}^{i}=\frac{e^{k\intop_{t_{i-1}}^{t_{i}}A(\mathbf{x}_{n}(t))dt}}{R_{i}}\quad\mbox{with
the normalisation factor}\quad
R_{i}=\frac{1}{N}\sum_{n=1}^{N}e^{k\int_{t_{i-1}}^{t_{i}}A(\mathbf{x}_{n}(t))dt}$
(10)

so that $\sum_{n=1}^{N}W_{n}^{i}=N$. The weights $\\{W_{n}^{i}\\}_{n=1\cdots
N}$ determine how many copies of each trajectory are made at time $t=t_{i}$.
The parameter $k$ characterizes the amplitude of the statistical bias involved
in the algorithm (see Fig. 16). For more information about the practical
implementation of the algorithm, the interested reader can refer to (Brewer et
al., 2018; Lestang, 2018). The application of this re-sampling at each step
$t_{i}$ eventually leads to a biased sampling in the trajectory space; the
trajectories corresponding to extreme values of
$\int_{0}^{T_{a}}A(\mathbf{x}_{n}(t))dt$ have a higher probability. The
sampled biased distribution writes

$\displaystyle\mathbb{P}_{k}\left(\left\\{\mathbf{X}(t)\right\\}_{0\leq t\leq
T_{a}}=\left\\{\mathbf{x}(t)\right\\}_{0\leq t\leq T_{a}}\right)$
$\displaystyle\underset{N\rightarrow\infty}{\sim}\frac{e^{k\int_{0}^{T_{a}}A(\mathbf{x}(t))dt}}{Z(k,T_{a})}\mathbb{\mathbb{P}}_{0}\left(\left\\{\mathbf{X}(t)\right\\}_{0\leq
t\leq T_{a}}=\left\\{\mathbf{x}(t)\right\\}_{0\leq t\leq T_{a}}\right),$ (11)

where $\mathbb{P}_{0}\left(\left\\{\mathbf{X}(t)\right\\}_{0\leq t\leq
T_{a}}=\left\\{\mathbf{x}(t)\right\\}_{0\leq t\leq T_{a}}\right)$ refers
formally to the probability of observing the trajectory
$\left\\{\mathbf{x}(t)\right\\}_{0\leq t\leq T_{a}}$. The normalisation factor
is given by $Z(k,T_{a})=\prod_{i=1}^{T_{a}/\tau}R_{i}$. One can mention that

$Z(k,T_{a})\underset{N\to\infty}{\sim}\mathbb{E}_{0}\left[e^{k\int_{0}^{T_{a}}A(\mathbf{X}(t))dt}\right],$
(12)

with $\mathbb{E}_{0}$ being the expectation value with respect to the
distribution $\mathbb{P}_{0}$. This result relies on the mean-field
approximation

$R_{i}=\frac{1}{N}\sum_{n=1}^{N}e^{k\int_{t_{i-1}}^{t_{i}}A(\mathbf{X}_{n}(t))dt}\underset{N\rightarrow\infty}{\sim}Z(k,t_{i})=\mathbb{E}_{i}\left[e^{k\int_{t_{i-1}}^{t_{i}}A(\mathbf{X}(t))dt}\right],$
(13)

where $\mathbb{E}_{i}[.]$ denotes the expectation value with respect to the
biased distribution $\mathbb{P}_{k}^{(i)}$ obtained after $i$ cloning steps.
The typical relative error related to this approximation can be shown to be of
order $1/\sqrt{N}$ for a family of rare-event algorithms including the gktl
algorithm (Del Moral, 2004, 2013). Rejected trajectories are discarded from
the statistics. Eventually, an effective ensemble of $N$ trajectories of
duration $T_{a}$ is obtained, distributed according to $\mathbb{P}_{k}$.

A key feature of the gktl algorithm is that all resampled trajectories are
solutions of the original dynamics. Nevertheless, it should be noted that a
small random perturbation is introduced in the cloning procedure to force
clones from the same trajectory to separate, _i.e._ artificial randomness is
introduced so that the cloning procedure is effective for deterministic
dynamics. As for the ams algorithm, it has been checked that this perturbation
did not affect the statistics of the sampled trajectories. Eventually, the
sampled trajectories obtained with the gktl algorithm can be used to compute
the statistical properties of any observable with respect to the distribution
$\mathbb{P}_{0}$ from the distribution $\mathbb{P}_{k}$ by using Eq. (11).

###  C.2 Computation of return times

Each trajectory in the biased ensemble results in a timeseries of the time-
averaged drag

$F_{T}^{(j)}(t)=\int_{t-T}^{t}f_{d}^{(j)}(\tau)d\tau,\quad t\in[T,T_{a}]$ (14)

and the return time of a fluctuation $F_{T}\geq a$ is given by (Lestang et
al., 2018)

$r(a)=-\frac{T_{a}-T}{\ln(1-\mathbb{P}(F_{T}\geq a))}.$ (15)

The probability $\mathbb{P}(F_{T}\geq a)$ can be estimated from the biased
ensemble by inverting Eq. (11)

$\mathbb{P}(f_{d}\geq
a)\approx\frac{1}{N}\sum_{j=1}^{N}e^{T_{a}\lambda(k)}e^{kT_{a}F_{T}^{(j)}}s_{j}(a)$
(16)

with $s_{j}(a)=1$ if $\max_{T\leq t\leq T_{a}}[F_{T}^{(j)}]\geq a$ and
$s_{j}(a)=0$ otherwise, _i.e._ by summing the weights of the timeseries which
maximum is larger than $a$.

##  Appendix D Perturbation at branching time

The application of the ams and gktl algorithms to a _deterministic_ dynamical
system requires a perturbation at initial time of the resampled trajectories.
In the absence of such perturbation, the resampling would not yield new
trajectories but exact copies of the original trajectories. In the framework
of Lattice Boltzmann simulations, the state of the system is described at a
mesoscopic level by the particle densities
$\\{f_{i}(\mathbf{x},t)\\}_{i=0\cdots 8}$ (see appendix A). Therefore, the
perturbation of the solution at time $t_{0}$ applies directly to the $f_{i}$’s
with

$f_{i}(\mathbf{x},t_{0})\longrightarrow
f_{i}(\mathbf{x},t_{0})+\epsilon\sum_{n=1}^{N_{s}}\alpha_{n}f_{i}^{(n)}(\mathbf{x}),$
(17)

where the $\alpha_{n}$’s are random numbers uniformly picked in the interval
$[0,1]$, $\epsilon$ is the relative amplitude of the perturbation and
$\\{f_{i}^{(n)}\\}_{n=1\cdots N_{s}}$ is a set of arbitrary snapshots of the
flow. In other words, the perturbation is a random linear combination of
$N_{s}$ snapshots. Eventually, the perturbed densities are rescaled so that
the mass of the system is preserved. In practice, we chose $\epsilon=0.002$
and $N_{s}=10$ to ensure that the perturbations remain sufficiently small,
random and independent. In order to check that this perturbation does not
impact the statistics of drag fluctuations, a long simulation of duration
$T_{tot}=10^{5}\tau_{c}$ has been performed with a periodic perturbation (with
period $\tau_{c}/2$) mimicking the perturbation of the clones in the algorithm
(see section 4). Fig. 23 shows that the statistics of the perturbed and
unperturbed simulations are equivalent.

Figure 23: pdfs of (normalised) drag fluctuations obtained from the perturbed
(at branching points) and unperturbed numerical simulations; $\epsilon$ is the
relative amplitude of the perturbation.

## References

  * Aldous & Vazirani (1994) Aldous, D. & Vazirani, U. 1994 "go with the winners" algorithms. In Proceedings 35th Annual Symposium on Foundations of Computer Science, pp. 492–501.
  * Aristoff et al. (2015) Aristoff, David, Lelièvre, Tony, Mayne, Christopher G. & Teo, Ivan 2015 Adaptive Multilevel Splitting in Molecular Dynamics Simulations. ESAIM: Proceedings and Surveys 48, 215–225.
  * Bec & Khanin (2007) Bec, Jeremie & Khanin, Konstantin 2007 Burgers Turbulence. Physics Reports 447 (1-2), 1–66, arXiv: 0704.1611.
  * Belenky et al. (2016) Belenky, Vadim, Weems, Kenneth & Lin, Woei-Min 2016 Split-time method for estimation of probability of capsizing caused by pure loss of stability. Ocean Engineering 122.
  * Bhatnagar et al. (1954) Bhatnagar, P. L., Gross, E. P. & Krook, M. 1954 A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems. Physical Review 94 (3), 511.
  * Blonigan et al. (2019) Blonigan, Patrick J., Farazmand, Mohammad & Sapsis, Themistoklis P. 2019 Are extreme dissipation events predictable in turbulent fluid flows? Phys. Rev. Fluids 4, 044606.
  * Bolhuis (2005) Bolhuis, Peter G 2005 Kinetic pathways of $\beta$-hairpin (un) folding in explicit solvent. Biophysical journal 88 (1), 50–61.
  * Bouchet et al. (2014) Bouchet, Freddy, Laurie, Jason & Zaboronski, Oleg 2014 Langevin dynamics, large deviations and instantons for the quasi-geostrophic model and two-dimensional euler equations. Journal of Statistical Physics 156 (6), 1066–1092.
  * Bouchet et al. (2019) Bouchet, Freddy, Rolland, Joran & Simonnet, Eric 2019 Rare event algorithm links transitions in turbulent flows with activated nucleations. Physical review letters 122 (7), 074502\.
  * Brewer et al. (2018) Brewer, Tobias, Clark, Stephen R, Bradford, Russell & Jack, Robert L 2018 Efficient characterisation of large deviations using population dynamics. Journal of Statistical Mechanics: Theory and Experiment 2018 (5), 053204.
  * Cérou & Guyader (2007) Cérou, Frédéric & Guyader, Arnaud 2007 Adaptive Multilevel Splitting for Rare Event Analysis. Stochastic Analysis and Applications 25 (2), 417–443.
  * Cérou et al. (2019) Cérou, Frédéric, Guyader, Arnaud & Rousset, Mathias 2019 Adaptive multilevel splitting: Historical perspective and recent results. Chaos: An Interdisciplinary Journal of Nonlinear Science 29 (4), 043108.
  * Chernykh & Stepanov (2001) Chernykh, A. I. & Stepanov, M. G. 2001 Large negative velocity gradients in Burgers turbulence. Physical Review E 64 (2), 026306.
  * Del Moral (2004) Del Moral, Pierre 2004 Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications. Springer-Verlag, New York.
  * Del Moral (2013) Del Moral, Pierre 2013 Mean field simulation for Monte Carlo integration, Monographs on Statistics and Applied Probability, vol. 126. CRC Press, Boca Raton, FL.
  * Dematteis et al. (2019) Dematteis, Giovanni, Grafke, Tobias, Onorato, Miguel & Vanden-Eijnden, Eric 2019 Experimental evidence of hydrodynamic instantons: The universal route to rogue waves. arXiv preprint arXiv:1907.01320 .
  * Dematteis et al. (2018) Dematteis, Giovanni, Grafke, Tobias & Vanden-Eijnden, Eric 2018 Rogue waves and large deviations in deep sea. Proceedings of the National Academy of Sciences 115 (5), 855–860.
  * Donzis & Sreenivasan (2010) Donzis, D. A. & Sreenivasan, K. R. 2010 Short-term forecasts and scaling of intense events in turbulence. Journal of Fluid Mechanics 647, 13–26.
  * Düben et al. (2008) Düben, P., Homeier, D., Jansen, K., Mesterhazy, D., Münster, G. & Urbach, C. 2008 Monte Carlo simulations of the randomly forced Burgers equation. EPL (Europhysics Letters) 84 (4), 40002.
  * van Erp & Bolhuis (2005) van Erp, Titus S. & Bolhuis, Peter G. 2005 Elaborating transition interface sampling methods. Journal of Computational Physics 205 (1), 157–181.
  * Escobedo et al. (2009) Escobedo, Fernando A., Borrero, Ernesto E. & Araque, Juan C. 2009 Transition path sampling and forward flux sampling. Applications to biological systems. Journal of Physics: Condensed Matter 21 (33), 333101.
  * Frisch (1995) Frisch, Uriel 1995 Turbulence: The Legacy of A. N. Kolmogorov. Cambridge University Press.
  * Garrahan et al. (2007) Garrahan, Juan P, Jack, Robert L, Lecomte, Vivien, Pitard, Estelle, van Duijvendijk, Kristina & van Wijland, Frédéric 2007 Dynamical first-order phase transition in kinetically constrained models of glasses. Physical review letters 98 (19), 195702.
  * Giardinà et al. (2011) Giardinà, Cristian, Kurchan, Jorge, Lecomte, Vivien & Tailleur, Julien 2011 Simulating Rare Events in Dynamical Processes. Journal of Statistical Physics 145 (4), 787–811.
  * Giardinà et al. (2006) Giardinà, Cristian, Kurchan, Jorge & Peliti, Luca 2006 Direct Evaluation of Large-Deviation Functions. Physical Review Letters 96 (12), 120603.
  * Glasserman et al. (1999) Glasserman, Paul, Heidelberger, Philip, Shahabuddin, Perwez & Zajic, Tim 1999 Multilevel Splitting for Estimating Rare Event Probabilities. Operations Research 47 (4), 585–600.
  * Grafke et al. (2013) Grafke, Tobias, Grauer, Rainer & Schäfer, Tobias 2013 Instanton filtering for the stochastic Burgers equation. Journal of Physics A: Mathematical and Theoretical 46 (6), 062002\.
  * Grafke et al. (2015) Grafke, Tobias, Grauer, Rainer & Schäfer, Tobias 2015 The instanton method and its numerical implementation in fluid mechanics. Journal of Physics A: Mathematical and Theoretical 48 (33), 333001.
  * Grassberger (2002) Grassberger, Peter 2002 Go with the winners: A general monte carlo strategy. Computer Physics Communications 147 (1-2), 64–70.
  * Grassberger et al. (1998) Grassberger, Peter & others 1998 Perm: a monte carlo strategy for simulating polymers and other things. arXiv preprint cond-mat/9806321 .
  * Grigorio et al. (2017) Grigorio, L. S., Bouchet, F., Pereira, R. M. & Chevillard, L. 2017 Instantons in a Lagrangian model of turbulence. Journal of Physics A: Mathematical and Theoretical 50 (5), 055501.
  * Gurarie & Migdal (1996) Gurarie, Victor & Migdal, Alexander 1996 Instantons in the Burgers equation. Physical Review E 54 (5), 4908–4914.
  * Huber & Kim (1996) Huber, G. A. & Kim, S. 1996 Weighted-ensemble Brownian dynamics simulations for protein association reactions. Biophysical Journal 70 (1), 97–110.
  * Kahn & Harris (1951) Kahn, H. & Harris, T. E. 1951 Estimation of Particle Transmission by Random Sampling. National Bureau of Standards 12, 27–30.
  * Kanev & van Engelen (2010) Kanev, Stoyan & van Engelen, Tim 2010 Wind turbine extreme gust control. Wind Energy 13 (1), 18–35, arXiv: https://onlinelibrary.wiley.com/doi/pdf/10.1002/we.338.
  * Kosztin et al. (1996) Kosztin, Ioan, Faber, Byron & Schulten, Klaus 1996 Introduction to the diffusion monte carlo method. American Journal of Physics 64 (5), 633–644.
  * Kraichnan & Chen (1989) Kraichnan, Robert H. & Chen, Shiyi 1989 Is there a statistical mechanics of turbulence? Physica D: Nonlinear Phenomena 37 (1), 160 – 172.
  * Krüger et al. (2017) Krüger, Timm, Kusumaatmaja, Halim, Kuzmin, Alexandr, Shardt, Orest, Silva, Goncalo & Viggen, Erlend Magnus 2017 The Lattice Boltzmann Method: Principles and Practice. Graduate Texts in Physics . Springer International Publishing.
  * Laffargue et al. (2013) Laffargue, Tanguy, Lam, Khanh-Dang Nguyen Thu, Kurchan, Jorge & Tailleur, Julien 2013 Large deviations of lyapunov exponents. Journal of Physics A: Mathematical and Theoretical 46 (25), 254002.
  * Latt et al. (2008) Latt, Jonas, Chopard, Bastien, Malaspinas, Orestis, Deville, Michel & Michler, Andreas 2008 Straight velocity boundaries in the lattice boltzmann method. Phys. Rev. E 77, 056703.
  * Laurie & Bouchet (2015) Laurie, Jason & Bouchet, Freddy 2015 Computation of rare transitions in the barotropic quasi-geostrophic equations. New Journal of Physics 17 (1), 015009.
  * Lecomte & Tailleur (2007) Lecomte, Vivien & Tailleur, Julien 2007 A numerical approach to large deviations in continuous time. Journal of Statistical Mechanics: Theory and Experiment 2007 (03), P03004.
  * Lesieur (2011) Lesieur, M. 2011 Turbulence in Fluids : Stochastic and Numerical Modelling. Springer.
  * Lestang (2018) Lestang, Thibault 2018 Numerical simulation and rare events algorithms for the study of extreme fluctuations of the drag force acting on an obstacle immersed in a turbulent flow. Theses, Université de Lyon.
  * Lestang et al. (2018) Lestang, Thibault, Ragone, Francesco, Bréhier, Charles-Edouard, Herbert, Corentin & Bouchet, Freddy 2018 Computing return times or return periods with rare event algorithms. Journal of Statistical Mechanics: Theory and Experiment 2018 (4), 043213.
  * Louvin et al. (2017) Louvin, Henri, Dumonteil, Eric & Lelièvre, Tony 2017 Three-dimensional neutron streaming calculations using adaptive multilevel splitting.
  * Mesterházy et al. (2013) Mesterházy, David, Biferale, Luca, Jansen, Karl & Tripiccione, Raffaele 2013 Lattice monte carlo methods for systems far from equilibrium. arXiv preprint arXiv:1311.4386 .
  * Mesterházy & Jansen (2011) Mesterházy, David & Jansen, Karl 2011 Anomalous scaling in the random-force-driven burgers’ equation: a monte carlo study. New Journal of Physics 13 (10), 103028.
  * Mohamad & Sapsis (2018) Mohamad, Mustafa A. & Sapsis, Themistoklis P. 2018 Sequential sampling strategy for extreme event statistics in nonlinear dynamical systems. Proceedings of the National Academy of Sciences 115 (44), 11138–11143, arXiv: https://www.pnas.org/content/115/44/11138.full.pdf.
  * Nemoto et al. (2017) Nemoto, Takahiro, Jack, Robert L. & Lecomte, Vivien 2017 Finite-size scaling of a first-order dynamical phase transition: Adaptive population dynamics and an effective model. Phys. Rev. Lett. 118, 115702.
  * Ragone et al. (2018) Ragone, Francesco, Wouters, Jeroen & Bouchet, Freddy 2018 Computation of extreme heat waves in climate models using a large deviation algorithm. Proceedings of the National Academy of Sciences 115 (1), 24–29.
  * Rodi (1998) Rodi, W. 1998 Large-eddy simulations of the flow past bluff bodies: State-of-the art. JSME 41, 361.
  * Rolland (2018) Rolland, Joran 2018 Extremely rare collapse and build-up of turbulence in stochastic models of transitional wall flows. Physical Review E 97 (2), 023109.
  * Rolland et al. (2016) Rolland, Joran, Bouchet, Freddy & Simonnet, Eric 2016 Computing transition rates for the 1-d stochastic ginzburg–landau–allen–cahn equation for finite-amplitude noise with a rare event algorithm. Journal of Statistical Physics 162 (2), 277–311.
  * Rolland & Simonnet (2015) Rolland, Joran & Simonnet, Eric 2015 Statistical behavior of adaptive multilevel splitting algorithms in simple models. Journal of Computational Physics 283, 541–558, arXiv: 1412.3362.
  * Rosis (2016) Rosis, Alessandro De 2016 Non-orthogonal central moments relaxing to a discrete equilibrium: A d2q9 lattice boltzmann model. EPL (Europhysics Letters) 116 (4), 44003.
  * Succi (2018) Succi, S. 2018 The Lattice Boltzmann Equation: For Complex States of Flowing Matter. Oxford University Press.
  * Sukop & Jr (2006) Sukop, Michael C. & Jr, Daniel T. Thorne 2006 Lattice Boltzmann Modeling: An Introduction for Geoscientists And Engineers. Springer-Verlag.
  * Tailleur & Kurchan (2007) Tailleur, Julien & Kurchan, Jorge 2007 Probing rare physical trajectories with Lyapunov weighted dynamics. Nature Physics 3 (3), 203–207.
  * Teo et al. (2016) Teo, Ivan, Mayne, Christopher G., Schulten, Klaus & Lelièvre, Tony 2016 Adaptive Multilevel Splitting Method for Molecular Dynamics Calculation of Benzamidine-Trypsin Dissociation Time. Journal of Chemical Theory and Computation 12 (6), 2983 – 2989.
  * Touchette (2009) Touchette, Hugo 2009 The large deviation approach to statistical mechanics. Physics Reports 478 (1-3), 1–69, arXiv: 0804.0327.
  * Villen-Altamirano & Villen-Altamirano (1994) Villen-Altamirano, M. & Villen-Altamirano, J. 1994 RESTART: a straightforward method for fast simulation of rare events. In Proceedings of Winter Simulation Conference, pp. 282–289.
  * Wouters & Bouchet (2016) Wouters, Jeroen & Bouchet, Freddy 2016 Rare event computation in deterministic chaotic systems using genealogical particle analysis. Journal of Physics A: Mathematical and Theoretical 49 (37), 374002.
  * Yeung et al. (2015) Yeung, P. K., Zhai, X. M. & Sreenivasan, Katepalli R. 2015 Extreme events in computational turbulence. Proceedings of the National Academy of Sciences 112 (41), 12633–12638, arXiv: https://www.pnas.org/content/112/41/12633.full.pdf.
  * Zuckerman & Chong (2017) Zuckerman, Daniel M & Chong, Lillian T 2017 Weighted ensemble simulation: review of methodology, applications, and software. Annual review of biophysics 46, 43–57.