# An Energy-Based Lengthscale for Reduced Order Models of Turbulent Flows

Changhong Mou  Elia Merzari  Omer San  Traian Iliescu<EMAIL_ADDRESS>[
Department of Mathematics, University of Wisconsin-Madison, Madison, WI 53706,
USA  Department of Nuclear Engineering, The Pennsylvania State University,
University Park, PA 16802, USA  School of Mechanical and Aerospace
Engineering, Oklahoma State University, Stillwater, OK 74078, USA  Department
of Mathematics, Virginia Tech, Blacksburg, VA 24061, USA

###### Abstract

In this paper, we propose a novel reduced order model (ROM) lengthscale that
is constructed by using energy distribution arguments. The new energy-based
ROM lengthscale is fundamentally different from the current ROM lengthscales,
which are built by using dimensional arguments. To assess the novel, energy-
based ROM lengthscale, we compare it with a standard, dimensionality-based ROM
lengthscale in two fundamentally different types of models: (i) the mixing-
length ROM (ML-ROM), which is a ROM closure model; and (ii) the evolve-filter-
relax ROM (EFR-ROM), which is a regularized ROM. We test the four combinations
(i.e., ML-ROM and EFR-ROM equipped with the energy-based and dimensionality-
based lengthscales) in the numerical simulation of the turbulent channel flow
at $Re_{\tau}=395$. The numerical investigation yields the following
conclusions: (i) The new energy-based ROM lengthscale is significantly (almost
two orders of magnitude) larger than the standard dimensionality-based ROM
lengthscale. As a result, the energy-based lengthscale yields more stable ML-
ROMs and EFR-ROMs than the dimensionality-based lengthscale. (ii) The energy-
based lengthscale displays the correct asymptotic behavior with respect to the
ROM dimension, whereas the dimensionality-based lengthscale does not. (iii)
The energy-based lengthscale yields ML-ROMs and (when significant filtering is
effected) EFR-ROMs whose parameters are less sensitive (i.e., more robust)
than the parameters of the ML-ROMs and EFR-ROMs based on the dimensionality-
based lengthscale. The novel energy-based lengthscale could enable the
development of better scale-aware ROM strategies for flow-specific
applications and is expected to have long term applications in nuclear reactor
thermal-hydraulics.

###### keywords:

Reduced order model , lengthscale , mixing-length , evolve-filter-relax ,
turbulent channel flow

††journal: arXiv.org

url]sites.google.com/view/iliescu/

##  1 Introduction

Reduced order models (ROMs) are models whose dimensions are orders of
magnitude lower than the dimensions of full order models (FOMs), i.e.,
computational models constructed by using classical numerical discretizations
(e.g., finite element or finite volume methods). Because of their relatively
low-dimensionality, ROMs can be used as efficient alternatives to FOMs in
computationally intensive applications, e.g., flow control, shape
optimization, and uncertainty quantification.

In recent years, ROMs have received a great deal of interest in nuclear
engineering applications to generate models that can account for the
fundamental physics of key phenomena while maintaining a low computational
cost [1, 2]. Examples are efforts to develop modal and machine-learning based
ROMs for thermal stratification [3], which is recognized as critical for the
licensing of liquid metal reactors. ROMs have also found applications in the
modeling of parameterized coupled thermal-hydraulics and reactor physics
problems [4, 5]. Finally, another notable effort is the use of ROMs to develop
flow acceleration for the advection-diffusion equation [6].

The Galerkin ROM (G-ROM) framework has been often used in the numerical
simulation of fluid flows [7, 8]. The G-ROM is constructed as follows: First,
in an offline phase, the FOM is used to produce snapshots, which are then
utilized to construct a low-dimensional (i.e., $r\ll N$) ROM basis
$\\{\varphi_{1},\ldots,\varphi_{r}\\}$, where $r$ is the ROM dimension and $N$
is the FOM dimension. There are several strategies for constructing the ROM
basis, e.g., the proper orthogonal decomposition (POD) [7, 9] and the reduced
basis method [10, 11]. In this paper, we exclusively use the POD to construct
the ROM basis. Next, the ROM basis is used together with a Galerkin projection
to build the G-ROM, which can be written as follows:

$\displaystyle\overset{\bullet}{\boldsymbol{a}}={\bf F}(\boldsymbol{a}),$ (1)

where $\boldsymbol{a}$ is the vector of coefficients in the ROM approximation
$\sum_{i=1}^{r}a_{i}(t)\boldsymbol{\varphi}_{i}(\boldsymbol{x})$ of the
variable of interest, $\overset{\bullet}{\boldsymbol{a}}$ denotes the vector
of time derivatives of $\boldsymbol{a}$, and the vector $\bf F$ comprises the
ROM operators that are preassembled in the offline phase. In the online phase,
the G-ROM (1) is employed for parameter values and/or time intervals that are
different from those used in the training stage.

We emphasize that, when FOMs (i.e., classical numerical methods, e.g., finite
element, finite volume, or finite difference methods) are used for the spatial
discretization of the fluid flow equations (e.g., the Navier-Stokes
equations), the FOM lengthscale is generally defined as the meshsize, $h$, of
the spatial mesh. In contrast, the following natural question is, to our
knowledge, still open:

What is the ROM lengthscale, $\delta$?

To formulate the above question mathematically, we first assume that the
following FOM and ROM variables are given, as is generally the case when ROMs
are applied in practical settings:

  * 1.

ROM variables, e.g., the ROM dimension, $r$, the total number of available ROM
basis functions, $R$, the ROM basis functions, $\boldsymbol{\varphi}_{i}$, and
the corresponding eigenvalues, $\lambda_{i}$, in the eigenvalue problem solved
to construct the POD basis.

  * 2.

FOM variables, e.g., the FOM mesh size, $h$, the FOM solution,
$\boldsymbol{u}^{FOM}$, and the computational domain characteristic
lengthscale, $L$.

Given these input FOM and ROM variables, we then try to answer the above
question, i.e., to express the ROM lengthscale, $\delta$, as a function of the
given FOM and ROM variables.

To motivate the need for a ROM lengthscale, we point out that the lengthscale
is a fundamental notion in engineering, geophysical, and biomedical
applications, where it is used to characterize the resolved spatial scales
(i.e., to determine the size of the spatial scales approximated by the
computational model) in flows around cars or airplanes, ocean or atmospheric
flows, or blood flow in an artery, respectively. To further motivate the need
for a ROM lengthscale, we emphasize that there are ROMs that use a lengthscale
in their very definition. For example, in under-resolved simulations, i.e.,
when the number of ROM basis functions is not enough to accurately represent
the turbulent flow dynamics, ROM closures and stabilizations are often used
(see [8] for a review). There are several examples of ROM closure models that
are defined by using a ROM lengthscale: the mixing-length ROM [12, 7, 13], the
Smagorinsky ROM [13], the dynamic SGS ROM [13], and the eddy viscosity
variational multiscale ROM [13]. There are also several examples of ROM
stabilizations that are defined by using a ROM lengthscale: the evolve-filter-
relax ROM [14, 15, 16, 17, 18] and the Leray ROM [14, 19, 20, 21, 22]. A ROM
lengthscale is needed to construct any of these two fundamentally different
types of ROMs (i.e., closures and stabilizations). A ROM lengthscale could
also be useful in the preprocessing strategy advocated in [23, 24] as a means
to filter out the noise in the input data.

In this paper, we propose a novel ROM lengthscale, which is constructed by
using energy distribution arguments. To assess the new energy-based ROM
lengthscale, we compare it with the classical ROM lengthscale used in [12, 7,
13], which is based on fundamentally different, dimensional arguments. To
compare the two ROM lengthscales, we utilize them to build two different types
of ROMs for under-resolved simulations: (i) the mixing length ROM [12, 7, 13],
which is a ROM closure model that augments the standard G-ROM with a
correction term; and (ii) the evolve-filter-relax ROM [14], which is a
regularized ROM (Reg-ROM) that leverages ROM spatial filtering to increase the
ROM stability and accuracy. We test both the ROM closure and the Reg-ROM
equipped with both ROM lengthscales in the numerical simulation of the
turbulent channel flow at $Re_{\tau}=395$. We note that a preliminary
numerical investigation of the new ROM lengthscale was performed in [25].

The rest of the paper is organized as follows: In Section 2, we outline the
standard G-ROM, ML-ROM, and EFR-ROM. In Section 3, we define the new, energy-
based lengthscale and the standard dimensionality-based lengthscale. In
Section 4, we present results for our investigation of the ML-ROM and EFR-ROM
equipped with the two lengthscales in numerical simulation of the turbulent
channel flow at $Re_{\tau}=395$. Finally, in Section 5, we draw conclusions
and outline directions of future research.

##  2 Reduced Order Models

In this section, we outline the construction of the standard Galerkin ROM
(G-ROM), mixing-length ROM (ML-ROM), and evolve-filter-relax ROM (EFR-ROM). As
a mathematical model, we consider the incompressible Navier-Stokes equation
(NSE):

$\displaystyle\frac{\partial\boldsymbol{u}}{\partial
t}-Re^{-1}\Delta\boldsymbol{u}+\bigl{(}\boldsymbol{u}\cdot\nabla\bigr{)}\boldsymbol{u}+\nabla
p\,$ $\displaystyle=\,\boldsymbol{f},$ (2)
$\displaystyle\nabla\cdot\boldsymbol{u}\,$ $\displaystyle=\,0,$ (3)

where $\boldsymbol{u}=[u_{1},u_{2},u_{3}]^{\top}$ is the velocity vector
field, $p$ the pressure field, $Re$ the Reynolds number, and $\boldsymbol{f}$
the forcing term. The NSE are equipped with appropriate boundary and initial
conditions.

###  2.1 Galerkin ROM (G-ROM)

To build the G-ROM, we consider the centering trajectory of the flow,

$\displaystyle\boldsymbol{U}(\boldsymbol{x})=\frac{1}{T}\int_{t}^{t+T}\boldsymbol{u}(\boldsymbol{x},t)dt,$
(4)

and assume that the ROM velocity approximation can be written as follows:

$\displaystyle\boldsymbol{u}_{r}(\boldsymbol{x},t)\,=\,\boldsymbol{U}(\boldsymbol{x})+\sum_{j=1}^{r}a_{j}(t)\boldsymbol{\varphi}_{j}(\boldsymbol{x})\,,$
(5)

where $\\{\boldsymbol{\varphi}_{j}\\}_{j=1}^{r}$ are the ROM basis functions
and $\boldsymbol{a}=[a_{1},\cdots,a_{r}]^{\top}$ are the sought ROM
coefficients. In our numerical experiments, we use the proper orthogonal
decomposition (POD) [7] to construct the ROM basis, but other ROM bases could
be used [26, 10, 11]. The next step in the G-ROM construction is to replace
$\boldsymbol{u}$ with $\boldsymbol{u}_{r}$ in (2) and project the resulting
equations onto the space spanned by the ROM basis,
$\\{\boldsymbol{\varphi}_{j}\\}_{j=1}^{r}$. This yields the G-ROM:

$\displaystyle\left(\frac{\partial\boldsymbol{u}_{r}}{\partial
t},\boldsymbol{\varphi}_{i}\right)+\bigl{(}(\boldsymbol{u}_{r}\cdot\nabla)\boldsymbol{u}_{r},\boldsymbol{\varphi}_{i}\bigr{)}+Re^{-1}\bigl{(}\nabla\boldsymbol{u}_{r},\nabla\boldsymbol{\varphi}_{i}\bigr{)}=\bigl{(}\boldsymbol{f},\boldsymbol{\varphi}_{i}\bigr{)}\,,\
\ i=1,\cdots,r\,,$ (6)

where $(\cdot,\cdot)$ denotes the $L^{2}$ inner product. The G-ROM can be
written as the following dynamical system for the vector of time coefficients,
$\boldsymbol{a}(t)$:

$\displaystyle\overset{\bullet}{\boldsymbol{a}}=\boldsymbol{b}+A\boldsymbol{a}+\boldsymbol{a}^{\top}B\boldsymbol{a}\,,$
(7)

where

$\displaystyle\boldsymbol{b}_{i}\,=\,\bigl{(}\boldsymbol{\varphi}_{i},\boldsymbol{f}\bigr{)}-\bigl{(}\boldsymbol{\varphi}_{i},\boldsymbol{U}\cdot\nabla\boldsymbol{U}\bigr{)}-Re^{-1}\bigl{(}\nabla\boldsymbol{\varphi}_{i},\nabla\boldsymbol{U}\bigr{)}\,,$
(8)
$\displaystyle\boldsymbol{A}_{im}\,=\,-\bigl{(}\boldsymbol{\varphi}_{i},\boldsymbol{U}\cdot\nabla\boldsymbol{\varphi}\bigr{)}-\bigl{(}\boldsymbol{\varphi}_{i},\boldsymbol{\varphi}_{m}\cdot\nabla\boldsymbol{U}\bigr{)}-Re^{-1}\bigl{(}\nabla\boldsymbol{\varphi}_{i},\nabla\boldsymbol{\varphi}_{m}\bigr{)}\,,$
(9)
$\displaystyle\boldsymbol{B}_{imn}\,=\,-\bigl{(}\boldsymbol{\varphi}_{i},\boldsymbol{\varphi}_{m}\cdot\nabla\boldsymbol{\varphi}_{n}\bigr{)}\,.$
(10)

###  2.2 Mixing-Length ROM (ML-ROM)

The G-ROM (7) is computationally efficient and relatively accurate in the
numerical simulation of laminar flows. However, the G-ROM generally yields
inaccurate results in the numerical simulation of turbulent flows. The main
reason for the G-ROM’s inaccuracy is that it is used in the under-resolved
regime, i.e., when the number of ROM basis functions, $r$, is not large enough
to accurately represent the complex dynamics of the turbulent flow. Thus, in
those cases, the G-ROM can be equipped with a ROM closure model, which models
the effect of the discarded ROM modes
$\\{\boldsymbol{\varphi}_{r+1},\ldots\\}$ on the G-ROM dynamics. In general,
the G-ROM with a closure model can be written as

$\displaystyle\overset{\bullet}{\boldsymbol{a}}=\boldsymbol{b}+A\boldsymbol{a}+\boldsymbol{a}^{\top}B\boldsymbol{a}+\boldsymbol{\tau}\,,$
(11)

where $\boldsymbol{\tau}$ is the ROM closure model. The current ROM closure
models are carefully surveyed in [8]. Some of these ROM closure models are
inspired from classical large eddy simulation (LES) closure modeling [27, 28].
To construct these LES-ROM closure models, one needs to define a ROM
lengthscale, which represents the size of the spatial scales modeled in the
LES-ROM. There are only a few ROM lengthscales in current use. In Section 3,
we define a novel ROM lengthscale. To assess this new ROM lengthscale, we
consider one of the simplest ROM closure models, the ML-ROM [7, 13], in which
the ROM closure term $\boldsymbol{\tau}$ in (11) is written as

$\displaystyle\boldsymbol{\tau}=-\bigl{(}\alpha\,U_{ML}\,\delta\bigr{)}\,S_{r}\,\boldsymbol{a}\,,$
(12)

where $\delta$ is one of the two ROM lengthscales defined in Section 3,
$U_{ML}$ is a characteristic velocity scale, $\alpha$ is a constant, and
$S_{r}$ is the ROM stiffness matrix with entries
$(S_{r})_{ij}=\bigl{(}\nabla\boldsymbol{\varphi}_{i},\nabla\boldsymbol{\varphi}_{j}\bigr{)},\
i,j=1,\ldots,r$. The ML-ROM model (12) is a functional closure model, which
aims at increasing the ROM viscosity in order to dissipate energy and mimic
the effect of the discarded modes [29]. The ML-ROM (12) was first used in [12,
7] and was further investigated in [13].

###  2.3 Evolve-Filter-Relax ROM (EFR-ROM)

Regularized ROMs (Reg-ROMs) [20, 14] represent an alternative to ROM closures
(e.g., the ML-ROM outlined in Section 2.2) in under-resolved simulations of
turbulent flows. Instead of adding a closure term, $\boldsymbol{\tau}$, as in
ROM closure modeling (see (11)), Reg-ROMs are constructed by using ROM spatial
filtering of various terms in the NSE to increase the ROM numerical stability.
Although regularized models have been used for decades in classical CFD [30,
31, 32], Reg-ROMs have only been recently developed [14].

The evolve-filter-relax ROM (EFR-ROM) is one of the most popular Reg-ROMs.
EFR-ROM is a modular ROM stabilization strategy that consists of three steps:
In the first step, which is called the evolve step, the standard G-ROM is used
to advance the current EFR-ROM time iteration, $\boldsymbol{a}^{n}$, to an
intermediate approximation, $\boldsymbol{w}^{n+1}$. In the second step, which
is called the filter step, the intermediate approximation,
$\boldsymbol{w}^{n+1}$, is filtered with the ROM differential filter [14],
which yields the filtered intermediate approximation,
$\overline{\boldsymbol{w}}^{n+1}$. In the third step, which is called the
relax step, the EFR-ROM approximation at the next time step,
$\boldsymbol{a}^{n+1}$, is calculated as the convex combination of the
intermediate approximation, $\boldsymbol{w}^{n+1}$, and the filtered
intermediate approximation, $\overline{\boldsymbol{w}}^{n+1}$. The EFR-ROM is
summarized in the following algorithm:

(I) $\displaystyle\text{\emph{
Evolve}:}\qquad\boldsymbol{a}^{n}\quad\stackrel{{\scriptstyle
G-ROM~{}\eqref{eqn:g-rom-U}}}{{\longmapsto}}\quad\boldsymbol{w}^{n+1}$ (II)
$\displaystyle\text{\emph{
Filter:}}\qquad\left(M_{r}+\gamma\,\delta^{2}\,S_{r}\right)\overline{\boldsymbol{w}}^{n+1}=\boldsymbol{w}^{n+1}$
(III)
$\displaystyle\text{\emph{Relax:}}\qquad\boldsymbol{a}^{n+1}=(1-\chi)\,\boldsymbol{w}^{n+1}+\chi\,\overline{\boldsymbol{w}}^{n+1}\,,$

where $\chi\in[0,1]$ is a relaxation parameter. In Step (II), we use the ROM
differential filter (DF) with an explicit ROM lengthscale, $\delta$, which
represents the filtering radius. The DF acts as a spatial filter by
eliminating the small scales (i.e., high frequencies) from the input data
[27]. We note that, in Step (II), we modify the classical DF [27] by
introducing a new parameter $\gamma$. This new parameter $\gamma$ has a role
similar to that of the parameter $\alpha$ used in the ML-ROM (12): It controls
the amount of filtering used in the DF. Step (III) is a relaxation step in
which the EFR-ROM velocity approximation at the new time step is defined as a
convex combination of the approximations obtained in Step (I) and Step (II).
The relaxation parameter $\chi$ diminishes the magnitude of the numerical
diffusion [33, 30, 32] and increases the accuracy [34, 33]. The scaling
$\chi\sim\Delta t$, where $\Delta t$ is the time step size, is a popular
choice [33], but higher values have also been used (see, e.g., [34]). To our
knowledge, the EFR-ROM was first used in [14] without the relaxation step, and
in [15] with the relaxation step (see also [18] and references therein). EFR-
ROM was also investigated in [14, 15, 16, 17, 18].

##  3 ROM Lengthscales

In this section, we present two different ROM lengthscales: In Section 3.1, we
present a standard ROM lengthscale, denoted $\delta_{1}$, which is constructed
by using dimensional analysis arguments. In Section 3.2, we propose a novel
ROM lengthscale, denoted $\delta_{2}$, which is constructed by using energy
balance arguments. As explained in the introduction, both definitions aim at
expressing the ROM lengthscale as a function of the following two types of
input variables: (i) ROM variables (e.g., the ROM dimension, $r$, the total
number of ROM basis functions, $R$, the eigenvalues, $\lambda_{i}$, and the
ROM basis functions, $\boldsymbol{\varphi}_{i}$). (ii) FOM variables (e.g.,
the fine FOM mesh size, $h$, the FOM solution, $\boldsymbol{u}^{FOM}$, and the
computational domain characteristic lengthscale, $L$). Given these input FOM
and ROM variables, we then try to answer the following natural question: What
is the ROM lengthscale, $\delta$?

###  3.1 Standard dimensionality-based ROM Lengthscale $\delta_{1}$

In this section, we use dimensional analysis to construct the first ROM
lengthscale, $\delta_{1}$. To this end, we follow the approach used in Section
3.2 in [13], which, in turn, is based on the pioneering ML-ROM proposed in
[12] for a turbulent pipe flow.

To construct the ROM lengthscale $\delta_{1}$, we first define the
componentwise FOM velocity fluctuations, i.e., the unresolved component of the
velocity, which is computed by using FOM data:

$\displaystyle{u^{{}^{\prime}}_{i}}^{FOM}=\sum_{j=r+1}^{R}a_{j}^{FOM}\,\varphi_{j}^{i},\qquad
i=1,2,3,$ (13)

where $R$ is the total number of ROM modes, $\varphi_{j}^{i}$ are the
componentwise ROM basis functions, and
$a_{j}^{FOM}=(\boldsymbol{u}^{FOM},\boldsymbol{\varphi}_{j})$ are the ROM
coefficients computed by using FOM data. Using the componentwise FOM velocity
fluctuations ${u^{{}^{\prime}}_{1}}^{FOM},{u^{{}^{\prime}}_{2}}^{FOM}$, and
${u^{{}^{\prime}}_{3}}^{FOM}$ in the $x,y$, and $z$ directions, respectively,
we build the FOM velocity fluctuation vector field
${\boldsymbol{u}^{\prime}}^{FOM}=[{u^{{}^{\prime}}_{1}}^{FOM},{u^{{}^{\prime}}_{2}}^{FOM},{u^{{}^{\prime}}_{3}}^{FOM}]$.
Since ${\boldsymbol{u}^{\prime}}^{FOM}$ varies with time, we calculate the
time averaged value of ${\boldsymbol{u}^{\prime}}^{FOM}$, i.e.,

$\displaystyle\langle{\boldsymbol{u}^{\prime}}^{FOM}\rangle_{t}(\boldsymbol{x})=\frac{1}{M}\sum_{k=1}^{M}{\boldsymbol{u}^{\prime}}^{FOM}(\boldsymbol{x},t_{k})=\frac{1}{M}\sum_{k=1}^{M}\sum_{l=r+1}^{R}\biggl{(}\boldsymbol{u}^{FOM}(\cdot,t_{k}),\boldsymbol{\varphi}_{l}(\cdot)\biggr{)}\boldsymbol{\varphi}_{l}(\boldsymbol{x}),$
(14)

where $M$ is the number of snapshots.

To construct the ROM lengthscale $\delta_{1}$, we adapt equation (22) in [13]
to our computational setting (i.e., the turbulent channel flow in Section 4):

$\displaystyle\delta_{1}:=\left(\frac{\int_{0}^{L_{1}}\,\int_{0}^{L_{2}}\,\int_{0}^{L_{3}}\sum_{i=1}^{3}{u^{{}^{\prime}}_{i}}^{FOM}{u^{{}^{\prime}}_{i}}^{FOM}\,dx_{1}\,dx_{2}\,dx_{3}}{\int_{0}^{L_{1}}\,\int_{0}^{L_{2}}\,\int_{0}^{L_{3}}\sum_{i=1}^{3}\sum_{j=1}^{3}\frac{\partial{u^{{}^{\prime}}_{i}}^{FOM}}{\partial
x_{j}\hfill}\,\frac{\partial{u^{{}^{\prime}}_{i}}^{FOM}}{\partial
x_{j}\hfill}\,dx_{1}\,dx_{2}\,dx_{3}}\right)^{1/2},$ (15)

where $L_{1},L_{2}$, and $L_{3}$ are the streamwise, wall-normal, and spanwise
dimensions of the computational domain of the turbulent channel flow test
problem, respectively.

Note that a quick dimensional analysis shows that the quantity defined in (15)
has the units of a lengthscale:

$\displaystyle[\delta_{1}]=\left(\frac{\frac{m}{s}\,\frac{m}{s}\,m^{3}}{\frac{1}{s}\,\frac{1}{s}\,m^{3}}\right)^{1/2}=m\,.$
(16)

We note that the ROM lengthscale, $\delta_{1}$, defined in (15), depends on
the FOM velocity fluctuation vector field, ${\boldsymbol{u}^{\prime}}^{FOM}$.

We also note that an alternative lengthscale was defined in equation (23) in
[13]. Since this alternative lengthscale was not used in the numerical
investigation in [13] (because it was harder to implement), we do not consider
it in this study.

###  3.2 Novel Energy-Based ROM Lengthscale $\delta_{2}$

In this section, we use energy balancing arguments and propose a new ROM
lengthscale, $\delta_{2}$. Noticing that the ROM truncation level, $r$, has
the role of dividing the kinetic energy of the system, we can require that the
new ROM lengthscale, $\delta_{2}$, do the same. Specifically, we require that
the ratio of kinetic energy contained in the first $r$ ROM modes,
$\sum_{i=1}^{r}\lambda_{i}$, to the kinetic energy contained in the total
number of ROM modes, $\sum_{i=1}^{R}\lambda_{i}$, be equal to the ratio of the
kinetic energy that can be represented on an imaginary mesh of size
$\delta_{2}$, $KE(\delta_{2})$, to the kinetic energy that can be represented
on the FOM mesh, $KE(h)$:

$\displaystyle\frac{\sum_{i=1}^{r}\lambda_{i}}{\sum_{i=1}^{R}\lambda_{i}}=\frac{KE(\delta_{2})}{KE(h)}\,.$
(17)

###### Remark 3.1

We emphasize that the mesh of size $\delta_{2}$ is not used in the actual ROM
construction. Instead, this imaginary mesh is used to highlight the
fundamental difference between the physical lengthscales of the space
generated by the first $r$ ROM modes,
$\boldsymbol{X}^{r}=\text{span}\\{\boldsymbol{\varphi}_{1},\ldots,\boldsymbol{\varphi}_{r}\\}$,
and the physical lengthscales of the space generated by all $R$ ROM modes,
$\boldsymbol{X}^{R}=\text{span}\\{\boldsymbol{\varphi}_{1},\ldots,\boldsymbol{\varphi}_{r},\boldsymbol{\varphi}_{r+1},\ldots,\boldsymbol{\varphi}_{R}\\}$.
Indeed, since the FOM data was used to build all the ROM modes, all the
lengthscales of the space $\boldsymbol{X}^{R}$ can be represented on the FOM
mesh. However, the imaginary mesh of size $\delta_{2}$ can represent only the
lengthscales of the space $\boldsymbol{X}^{r}$; it cannot be expected to
represent the physical lengthscales of the space spanned by the higher index
ROM modes $\\{\boldsymbol{\varphi}_{r+1},\ldots,\boldsymbol{\varphi}_{R}\\}$,
which are generally associated with the small scales.

To compute the ratio $\frac{KE(\delta_{2})}{KE(h)}$ in (17), we transfer the
problem to the usual Fourier space. To this end, we first notice that
$\delta_{2}$ defines a cutoff wavenumber:

$\displaystyle k_{\delta_{2}}:=\frac{2\,\pi}{\delta_{2}}\,.$ (18)

We then notice that the kinetic energy in the system can be written in terms
of the energy spectrum, $E(\cdot)$:

$\displaystyle KE(k)=\int_{k_{0}}^{k}E(k^{\prime})\,dk^{\prime}\,,$ (19)

where $\displaystyle k_{0}=\frac{2\,\pi}{L}$ is the Fourier wavenumber that
corresponds to the computational domain characteristic lengthscale, $L$. In
the case of isotropic, homogeneous turbulence, we have the usual energy
spectrum given by Kolmogorov’s theory [28, 35]

$\displaystyle E(k)\sim C\,\varepsilon^{2/3}\,k^{-5/3}\,.$ (20)

Thus, the condition imposed in (17) can be written as

$\displaystyle\frac{\int_{k_{0}}^{k_{\delta_{2}}}E(k)\,dk}{\int_{k_{0}}^{k_{h}}E(k)\,dk}=\frac{\sum_{i=1}^{r}\lambda_{i}}{\sum_{i=1}^{R}\lambda_{i}}\stackrel{{\scriptstyle\text{notation}}}{{=}}\Lambda\,,$
(21)

where $\displaystyle k_{h}=\frac{2\,\pi}{h}$ is the highest Fourier wavenumber
that can be resolved on the given FOM meshsize, $h$. The LHS of (21) can be
evaluated by using (20):

$\displaystyle\int_{k_{0}}^{k_{\delta_{2}}}E(k)\,dk=C\,\varepsilon^{2/3}\,\int_{k_{0}}^{k_{\delta_{2}}}k^{-5/3}\,dk=C\,\varepsilon^{2/3}\,\frac{k_{\delta_{2}}^{-2/3}-k_{0}^{-2/3}}{-2/3}\,,$
(22)

and, similarly,

$\displaystyle\int_{k_{0}}^{k_{h}}E(k)\,dk=C\,\varepsilon^{2/3}\,\int_{k_{0}}^{k_{h}}k^{-5/3}\,dk=C\,\varepsilon^{2/3}\,\frac{k_{h}^{-2/3}-k_{0}^{-2/3}}{-2/3}.$
(23)

Plugging (22) and (23) back into (21), simplifying, and rearranging, we obtain

$\displaystyle
k_{\delta_{2}}^{-2/3}=\Lambda\,k_{h}^{-2/3}+\left(1-\Lambda\right)\,k_{0}^{-2/3}\,.$
(24)

Since $1\leq r\leq R$, $\Lambda$ satisfies the inequality $0<\Lambda\leq 1$.
Thus, (24) implies that $k_{\delta_{2}}^{-2/3}$ is a convex combination of
$k_{h}^{-2/3}$ and $k_{0}^{-2/3}$. Furthermore, as expected, $k_{\delta_{2}}$
satisfies the following asymptotic relations:

$\displaystyle\biggl{[}(r\rightarrow
R)\,\Longrightarrow\,(k_{\delta_{2}}\longrightarrow
k_{h})\biggr{]}\qquad\text{and}\qquad\biggl{[}(r\rightarrow
0)\,\Longrightarrow\,(k_{\delta_{2}}\longrightarrow k_{0})\biggr{]}\,.$ (25)

Using (24) together with (18), gives us a formula for $\delta_{2}$:

$\displaystyle\delta_{2}=\frac{2\,\pi}{k_{\delta_{2}}}=2\,\pi\,\left[\Lambda\,\left(\frac{2\,\pi}{h}\right)^{-2/3}+\left(1-\Lambda\right)\,\left(\frac{2\,\pi}{L}\right)^{-2/3}\right]^{3/2}=\left[\Lambda\,h^{2/3}+\left(1-\Lambda\right)\,L^{2/3}\right]^{3/2}\,.$
(26)

The new ROM lengthscale, $\delta_{2}$, defined in (26), depends on the FOM
mesh size, $h$, the ROM dimension, $r$, the total number of ROM basis
functions, $R$, the eigenvalues, $\lambda_{i}$, and the computational domain
characterisitic lengthscale, $L$.

######  Remark 3.2 (Asymptotic Behavior)

We note that the novel ROM lengthscale, $\delta_{2}$, satisfies the following
two natural asymptotic relations:

  * 1.

As $r$ approaches $R$, $\delta_{2}$ approaches $h$.

This simply says that, as expected, when the ROM dimension, $r$, approaches
the maximal ROM dimension, $R$, the ROM lengthscale, $\delta_{2}$, approaches
the minimal FOM lengthscale, $h$ (i.e., the spatial mesh size).

  * 2.

As $r$ approaches $1$, $\delta_{2}$ approaches $L$.

This simply says that, as expected, when the ROM dimension, $r$, approaches
the minimal ROM dimension, $1$, the ROM lengthscale, $\delta_{2}$, approaches
the maximal FOM lengthscale, $L$ (i.e., the dimension of the computational
domain).

##  4 Numerical Results

In this section, we perform a numerical investigation of the two lengthscales
discussed in Section 3: the standard dimensionality-based ROM lengthscale,
$\delta_{1}$, defined in (15), and the new energy-based ROM lengthscale,
$\delta_{2}$, defined in (26). To this end, we use two fundamentally different
ROMs: the ML-ROM presented in Section 2.2 and the EFR-ROM presented in Section
2.3. In each type of ROM, we use the two ROM lengthscales and compare the
results. To this end, in our numerical investigation, we test four types of
models: (i) ML-ROM1, which is the ML-ROM in which the velocity scale,
$U_{ML}$, is set equal to $\delta_{1}$; (ii) ML-ROM2, which is the ML-ROM in
which the velocity scale, $U_{ML}$, is set equal to $\delta_{2}$; (iii) EFR-
ROM1, which is the EFR-ROM in which the filter radius, $\delta$, is set equal
to $\delta_{1}$; and (iv) EFR-ROM2, which is the EFR-ROM in which the filter
radius, $\delta$, is set equal to $\delta_{2}$. Thus, to compare the two
lengthscales, we compare the four ROMs (i.e., ML-ROM1, ML-ROM2, EFR-ROM1, and
EFR-ROM2) in the numerical simulation of the 3D turbulent channel flow at
$Re_{\tau}=395$. We emphasize that the goal of this section is not to find the
best ML-ROMs and EFR-ROMs. Instead, we investigate whether the two
lengthscales are different and, if so, quantify their differences and how
these differences impact the ML-ROM and EFR-ROM results.

###  4.1 Numerical Setting

#### FOM

The computational domain is a rectangular box,
$\Omega=(-2\pi,2\pi)\times(0,2)\times(-2\pi/3,2\pi/3)$. We enforce no slip
boundary conditions on the walls at $y=0$ and $y=2$, and periodic boundary
conditions in the $x-$ and $z-$directions. We also use the forcing term
$\boldsymbol{f}=[1,0,0]^{\top}$ and the Reynolds number $Re_{\tau}=395$
($Re=13,750$). To generate the snapshots, we run an LES model using the
rNS-$\alpha$ scheme [36, 37] with the time step size $\Delta t=0.002$.

#### ROM

We collect a total of $5000$ snapshots from $t=60$ to $t=70$ (i.e., from a
time interval in which the flow is in the statistically steady regime) and use
the POD to generate the ROM basis. For illustrative purposes, we plot the
magnitude fields of the ROM basis functions
$\boldsymbol{\varphi}_{1},\boldsymbol{\varphi}_{25}$, and
$\boldsymbol{\varphi}_{50}$ in Fig. 1.

Figure 1:  The 1st, 25th, and 50th POD modes

For the ROM time discretization, we utilize the commonly used linearized BDF2
temporal discretization with a time step size $\Delta t=0.002$. As the ROM
initial conditions, we use the ROM projections of the LES approximations at
$t=60$ and $t=60.002$. For convenience, in our ROM simulations, $t=0$
corresponds to $t=60$ in the LES model.

To assess the ROMs’ performance, we use two different criteria: (i) the time
evolution of the kinetic energy, $E(t)$, and (ii) second-order statistics.

We define the ROM kinetic energy as follows:

$\displaystyle
E(t)=\frac{1}{2}\int_{\Omega}\bigl{(}u_{1}(\boldsymbol{x},t)^{2}+u_{2}(\boldsymbol{x},t)^{2}+u_{3}(\boldsymbol{x},t)^{2}\bigr{)}\,d\boldsymbol{x}\,,$
(27)

where $u_{1},u_{2}$, and $u_{3}$ are the components of the ROM velocity field
approximation.

Following [37], we consider the following two second-order statistics: (i) the
normalized root mean square (RMS) of the streamwise velocity component,
$U_{RMS}$:

$\displaystyle
U_{RMS}:=\frac{\biggl{|}\widetilde{\mathbb{R}}_{11}-\frac{1}{3}\sum_{j=1}^{3}\widetilde{\mathbb{R}}_{jj}\biggr{|}^{1/2}}{u_{1,\tau}},$
(28)

and (ii) the normalized streamwise-spanwise Reynolds stress tensor component,
$\mathbb{R}_{12}$:

$\displaystyle\mathbb{R}_{12}:=\frac{\widetilde{\mathbb{R}}_{12}}{u_{1,\tau}^{2}}.$
(29)

In these second-order statistics, the Reynolds stress tensor components are
calculated as follows:

$\displaystyle\widetilde{\mathbb{R}}_{ij}=\left\langle\left\langle
u_{i}u_{j}\right\rangle_{s}\right\rangle_{t}-\left\langle\left\langle
u_{i}\right\rangle_{s}\right\rangle_{t}\left\langle\left\langle
u_{j}\right\rangle_{s}\right\rangle_{t}\,,$ (30)

where $\langle\cdot\rangle_{s}$ denotes spatial averaging,
$\langle\cdot\rangle_{t}$ denotes time averaging, and $u_{i}$ are the
components of the given ROM or FOM velocity field approximations. The friction
velocity, $u_{1,\tau}$, which is used in (29), is calculated by using the
following formula:

$\displaystyle
u_{1,\tau}=\sqrt{\nu\,\frac{U_{\text{mean}}(y_{min})}{y_{min}}},$ (31)

where $\nu$ is the kinematic viscosity, $y_{min}$ is the minimum positive
$y$-value of the FOM mesh, and $U_{\text{mean}}=\left\langle\left\langle
u_{1}\right\rangle_{s}\right\rangle_{t}$ the average ROM or FOM velocity flow
profile.

###  4.2 Numerical Results: Lengthscale Comparison

In this section, we investigate the relative size of the two ROM lengthscales:
(i) the standard dimensionality-based lengthscale, $\delta_{1}$, defined in
(15), and (ii) the novel energy-based ROM lengthscale, $\delta_{2}$, defined
in (26). To calculate $\delta_{2}$ in equation (26), we define the FOM global
mesh size, $h$, as $h=\max_{K\in\mathcal{K}}r_{K}$, where the mesh
$\mathcal{K}$ is the set $\\{K\\}$ of tetrahedrons K and $r_{K}$ is the
inradius of the local tetrahedron, $K$. For the test problem used in our
numerical investigation, $h=1.1\times 10^{-1}$. Furthermore, we define the
maximal FOM lengthscale, $L$, as the vertical size of the computational
domain, i.e., $L=2$.

In Table 1, we list the $\delta_{1}$ and $\delta_{2}$ values for $r$ values
from $4$ to $50$. These results show that the two ROM lengthscales have very
different behaviors with respect to changes in $r$.

#### Magnitude Behavior

First, we notice that the $\delta_{2}$ magnitude is between one and two orders
of magnitude larger than the $\delta_{1}$ magnitude. Indeed, for the low $r$
values (i.e., $r=4,8,16$), $\delta_{2}$ is more than an order of magnitude
larger than $\delta_{1}$. For the large $r$ values (i.e., $r=32,40,50$),
$\delta_{2}$ is still an order of magnitude larger than $\delta_{1}$, but the
difference between the two lengthscales is smaller.

#### Asymptotic Behavior

Second, the asymptotic behavior of the ROM lengthscales with respect to $r$ is
fundamentally different: As $r$ increases, the $\delta_{1}$ magnitude remains
relatively unchanged. In contrast, as $r$ increases, the $\delta_{2}$
magnitude decreases by almost one order of magnitude. We emphasize that
$\delta_{2}$ has the natural asymptotic behavior expected from a ROM
lengthscale, as explained in Remark 3.2: (i) As $r$ increases toward its
maximal value, i.e., $R$ (which is $R=804$ in our numerical investigation),
$\delta_{2}$ approaches $h$. Indeed, the results in Table 1 show that for the
maximal $r$ value (i.e., for $r=50$), $\delta_{2}$ achieves the smallest
value, $\delta_{2}=4.32\times 10^{-1}$, which is the same order of magnitude
as the FOM mesh size, $h=1.1\times 10^{-1}$. (ii) As $r$ approaches its
minimal value (i.e., $1$), $\delta_{2}$ approaches $L$. Indeed, the results in
Table 1 show that for the minimal $r$ value (i.e., for $r=1$), $\delta_{2}$
achieves the largest value, $\delta_{2}=1.63\times 10^{0}$, which is the same
order of magnitude as the size of the computational domain, $L=2.0\times
10^{0}$.

$r$ | 4 | 8 | 16 | 32 | 40 | 50
---|---|---|---|---|---|---
$\delta_{1}$ | 4.64e-02 | 4.65e-02 | 4.68e-02 | 4.68e-02 | 4.66e-02 | 4.62e-02
$\delta_{2}$ | 1.63e00 | 1.41e+00 | 1.08e+00 | 6.84e-01 | 5.56e-01 | 4.32e-01
Table 1: ROM lengthscale values for different $r$ values.

###  4.3 Numerical Results: ML-ROM Investigation

In this section, we investigate the role played by the two ROM lengthscales,
$\delta_{1}$ and $\delta_{2}$, in the ML-ROM (12):

$\displaystyle\boldsymbol{\tau}=-\bigl{(}\alpha\,U_{ML}\,L_{ML}\bigr{)}\,S_{r}\,\boldsymbol{a}\,.$
(32)

We denote ML-ROM1 the ML-ROM in which $L_{ML}=\delta_{1}$ in (32) and ML-ROM2
the ML-ROM in which $L_{ML}=\delta_{2}$ in (32). To ensure a fair comparison
of ML-ROM1 and ML-ROM2, we use the same constant $\alpha$ and the same
velocity scale $U_{ML}$ (i.e., the time averaged streamwise velocity
component) in (32) for both models, and vary only the ROM lengthscale, i.e.,
$L_{ML}=\delta_{1}$ or $L_{ML}=\delta_{2}$. To vary the ROM lengthscale, we
vary the $r$ value in the definitions of $\delta_{1}$ and $\delta_{2}$. In
Figures 2–5 we plot the time evolution of the kinetic energy and the second-
order statistics of the ML-ROM1 and ML-ROM2 for different $r$ values and two
different $\alpha$ values: $\alpha=6\times 10^{-3}$ (Figures 2 and 4) and
$\alpha=6\times 10^{-4}$ (Figures 3and 5). (Results for more $\alpha$ values
are presented in the preliminary numerical investigation in [25].) As a
benchmark for the ROM results, we use the projection of the the FOM results on
the ROM basis (denoted as LES-proj in these plots).

Figure 2: Time evolution of the ML-ROM kinetic energy for $\alpha=6\times
10^{-3}$

Figure 3: Time evolution of the ML-ROM kinetic energy for $\alpha=6\times
10^{-4}$

(a) $r=4$

(b) $r=8$

(c) $r=16$

(d) $r=32$

(e) $r=40$

(f) $r=50$

Figure 4:  Second-order ML-ROM statistics for $\alpha=6\times 10^{-3}$

(a) $r=4$

(b) $r=8$

(c) $r=16$

(d) $r=32$

(e) $r=40$

(f) $r=50$

Figure 5:  Second-order ML-ROM statistics for $\alpha=6\times 10^{-4}$

#### Stability

Since the results in Table 1 show that $\delta_{2}$ is between one and two
orders of magnitude higher than $\delta_{1}$, we expect ML-ROM2 to yield more
stable results than ML-ROM1. Indeed, since we fix all the ML-ROM parameters
(i.e., the constant $\alpha$ and the velocity scale $U_{ML}$) and $\delta_{2}$
is significantly larger than $\delta_{1}$, we expect the ML-ROM2 artificial
viscosity to be higher than the ML-ROM1 artificial viscosity (and, thus, ML-
ROM2 to be more stable than ML-ROM1). This is clearly shown in the plots in
Figures 2–5, in which ML-ROM2 yields stable results for all $r$ values and for
both $\alpha$ values. In contrast, for the smallest $\alpha$ value,
$\alpha=6\times 10^{-4}$ (Figures 3 and 5), ML-ROM1 blows up for all $r$
values. Furthermore, for the largest $\alpha$ value, $\alpha=6\times 10^{-3}$
(Figures 2 and 4), ML-ROM1 blows up for the small $r$ values (i.e.,
$r=4,8,16$, and $32$). To quantify the stability of the two ML-ROMs, in Table
2, for different $r$ values, we list the threshold $\alpha_{0}$ value, i.e.,
the value that ensures that, if $\alpha>\alpha_{0}$, then the ML-ROM is
stable. These results show that, for each $r$ value, the threshold
$\alpha_{0}$ value is more than an order of magnitude lower for ML-ROM2 than
for ML-ROM1. Thus, we conclude that ML-ROM2 is more stable than ML-ROM1, which
is the same conclusion as that yielded by Figures 2–5.

| $r$ | 4 | 8 | 16 | 32 | 40 | 50
---|---|---|---|---|---|---|---
ML-ROM1 | $\alpha_{0}$ | 10e-3 | 9.8e-3 | 9.2e-3 | 8.5e-3 | 6.5e-3 | 6.2e-3
ML-ROM2 | $\alpha_{0}$ | 2.9e-4 | 3.4e-4 | 4.1e-4 | 6.7e-4 | 5.9e-4 | 7.5e-4
Table 2: ML-ROM threshold $\alpha_{0}$ values for different $r$ values.

#### Accuracy

Since $\delta_{2}$ is between one and two orders of magnitude higher than
$\delta_{1}$, we expect the ML-ROM accuracy to depend on the constant
$\alpha$. Indeed, the ML-ROM1 and ML-ROM2 plots in Figures 2–5 do not show a
clear winner: For the largest $\alpha$ value (i.e., $\alpha=6\times 10^{-3}$),
ML-ROM1 is more accurate than ML-ROM2 for large $r$ values (i.e., $r=40$ and
$50$) and less accurate for small $r$ values (i.e., $r=4,8$, and $16$). For
the smallest $\alpha$ value (i.e., $\alpha=6\times 10^{-4}$), ML-ROM2 is more
accurate than ML-ROM1 for all $r$ values (since ML-ROM1 simply blows up). We
also note that, for $\alpha=6\times 10^{-4}$, ML-ROM2 is quite accurate for
larger $r$ values (i.e., $r=32,40$, and $50$).

#### Parameter Sensitivity

Finally, we investigate the ML-ROM’s parameter sensitivity. Specifically, we
investigate which ROM lengthscale yields ML-ROMs that are less sensitive
(i.e., more robust) with respect to the ML-ROM’s $\alpha$ parameter. We
emphasize that the model sensitivity with respect to model parameters is a
well known issue that has hindered the development of closures and
stabilizations in CFD over the years [28, 27, 31]. Thus, finding robust (i.e.,
less sensitive) ML-ROMs that require as little parameter tuning as possible is
an important practical problem.

To quantify the ML-ROM’s parameter sensitivity, in Table 3, for different $r$
values, we list the optimal $\alpha$ value in ML-ROM, i.e., the $\alpha$ value
that ensures that the average ROM kinetic energy, $KE^{ROM}$, is the closest
to the average FOM (LES) kinetic energy, $KE^{LES}$. Specifically, we solve
the following optimization problem:

$\displaystyle\min_{\alpha}\left|\overline{KE}^{ROM}-\overline{KE}^{LES}\right|,$
(33)

where $\overline{KE}=\frac{1}{M}\sum_{k=1}^{M}KE(t_{k})$, and $M$ is the
number of snapshots.

| $r$ | 4 | 8 | 16 | 32 | 40 | 50
---|---|---|---|---|---|---|---
ML-ROM1 | $\alpha$ | 1.50e-2 | 1.38e-2 | 1.38e-2 | 1.28e-2 | 7.80e-3 | 6.51e-3
ML-ROM2 | $\alpha$ | 4.35e-4 | 4.42e-4 | 5.33e-4 | 6.70e-4 | 6.20e-4 | 7.50e-4
Table 3: ML-ROM optimal $\alpha$ values for different $r$ values.

The results in Table 3 show that the optimal ML-ROM $\alpha$ values are
sensitive with respect to changes in $r$. Indeed, as $r$ increases from $4$ to
$50$, $\alpha$ decreases by a factor of $2.3$ for ML-ROM1 and increases by a
factor of $1.7$ for ML-ROM2. Overall, the results in Table 3 show that the ML-
ROM’s parameter $\alpha$ sensitivity is higher for ML-ROM1 than for ML-ROM2.

###  4.4 Numerical Results: EFR-ROM Investigation

To further investigate the role played by the two ROM lengthscales,
$\delta_{1}$ and $\delta_{2}$, in reduced order modeling of turbulent flows,
in this section we consider the EFR-ROM presented in Section 2.3. We emphasize
that the EFR-ROM is completely different from the ML-ROM investigated in
Section 4.3: The EFR-ROM is a Reg-ROM based on numerical stabilization,
whereas the ML-ROM is a closure model. Thus, the EFR-ROM investigation in this
section could shed new light on the novel ROM lengthscale, $\delta_{2}$.

To ensure a fair comparison of the two ROM lengthscales, $\delta_{1}$ and
$\delta_{2}$, we fix the parameters $\gamma$ and $\chi$ in the EFR-ROM in
Section 2.3 and change only the filter radius, $\delta$, of the differential
filter used in Step (II) of the EFR-ROM algorithm. We denote the resulting
models as EFR-ROM1 (when $\delta=\delta_{1}$) and EFR-ROM2 (when
$\delta=\delta_{2}$). In our numerical investigation, we fix the $\chi$ value
for both EFR-ROM1 and EFR-ROM2 to $\chi=6\times 10^{-3}$, which is the minimum
value that yields a stable solution for $r=4$ EFR-ROM1. We also fix the
$\gamma$ value for both EFR-ROM1 and EFR-ROM2. In our numerical investigation,
we consider two $\gamma$ values. To vary the ROM lengthscale, we vary the $r$
value in the definitions of $\delta_{1}$ and $\delta_{2}$.

In Figures 6–9, we plot the time evolution of the kinetic energy and the
second-order statistics of the EFR-ROM1 and EFR-ROM2 for different $r$ values
and two different $\gamma$ values: $\gamma=8\times 10^{-2}$ (Figures 6 and 8)
and $\gamma=9\times 10^{-1}$ (Figures 7and 9). As a benchmark for the ROM
results, we use the projection of the the FOM results on the ROM basis
(denoted as LES-proj in these plots).

Figure 6:  Time evolution of the EFR-ROM kinetic energy for $\gamma=8\times
10^{-2}$

Figure 7:  Time evolution of the EFR-ROM kinetic energy for $\gamma=9\times
10^{-1}$

(a) $r=4$

(b) $r=8$

(c) $r=16$

(d) $r=32$

(e) $r=40$

(f) $r=50$

Figure 8:  Second-order EFR-ROM statistics for $\gamma=8\times 10^{-2}$

(a) $r=4$

(b) $r=8$

(c) $r=16$

(d) $r=32$

(e) $r=40$

(f) $r=50$

Figure 9:  Second-order EFR-ROM statistics for $\gamma=9\times 10^{-1}$

#### Stability

Since the results in Table 1 show that $\delta_{2}$ is between one and two
orders of magnitude higher than $\delta_{1}$, we expect EFR-ROM2 to yield more
stable results than EFR-ROM1. Indeed, since we fix all the EFR-ROM parameters
(i.e., $\chi$ and $\gamma$) and $\delta_{2}$ is significantly larger than
$\delta_{1}$, we expect the EFR-ROM2 filtering level to be higher than the
EFR-ROM1 filtering level (and, thus, EFR-ROM2 to be more stable than EFR-
ROM1). This is clearly shown in the plots in Figures 6–9, in which EFR-ROM2
yields stable results for all $r$ values and for both $\gamma$ values. In
contrast, for the smallest $\gamma$ value, $\gamma=8\times 10^{-2}$ (Figures 6
and 8), EFR-ROM1 blows up for all $r$ values. Furthermore, for the largest
$\gamma$ value, $\gamma=9\times 10^{-1}$ (Figures 7 and 9), EFR-ROM1 blows up
for the small $r$ values (i.e., $r=4$ and $r=8$). To quantify the stability of
the two EFR-ROMs, in Table 4, for different $r$ values, we list the threshold
$\gamma_{0}$ value, i.e., the value that ensures that, if $\gamma>\gamma_{0}$,
then the EFR-ROM is stable. These results show that, for each $r$ value, the
threshold $\gamma_{0}$ value is more than an order of magnitude lower for EFR-
ROM2 than for EFR-ROM1. Thus, we conclude that EFR-ROM2 is more stable than
EFR-ROM1, which is the same conclusion as that yielded by Figures 6–9.

| $r$ | 4 | 8 | 16 | 32 | 40 | 50
---|---|---|---|---|---|---|---
EFR-ROM1 | $\gamma_{0}$ | 1.0e0 | 1.0e0 | 9.0e-1 | 8.0e-1 | 7.0e-1 | 7.0e-1
EFR-ROM2 | $\gamma_{0}$ | 2.9e-2 | 3.3e-2 | 3.9e-2 | 5.5e-2 | 5.9e-2 | 7.5e-2
Table 4: EFR-ROM threshold $\gamma_{0}$ values for different $r$ values.

#### Accuracy

Since $\delta_{2}$ is between one and two orders of magnitude higher than
$\delta_{1}$, we expect the EFR-ROM accuracy to depend on the constant
$\gamma$. Indeed, the EFR-ROM1 and EFR-ROM2 plots in Figures 6–9 do not
display a clear winner: For the largest $\gamma$ value (i.e., $\gamma=9\times
10^{-1}$), EFR-ROM2 is more accurate than EFR-ROM1 for the small $r$ values
(i.e., $r=4$ and $8$) since EFR-ROM1 simply blows up. For the remaining $r$
values, EFR-ROM1 and EFR-ROM2 display the same accuracy level. For the
smallest $\gamma$ value (i.e., $\gamma=8\times 10^{-2}$), EFR-ROM2 is more
accurate than EFR-ROM1 for all $r$ values (since EFR-ROM1 simply blows up). We
also note that EFR-ROM2 is relatively accurate for the largest $r$ value
(i.e., $r=50$).

#### Parameter Sensitivity

To study the EFR-ROM’s parameter sensitivity, we investigate which ROM
lengthscale yields EFR-ROMs that are less sensitive (i.e., more robust) with
respect to the EFR-ROM’s $\gamma$ parameter. To quantify the ML-ROM’s
parameter sensitivity, in Table 5, for different $r$ values, we list the
optimal $\gamma$ value in EFR-ROM, i.e., the $\gamma$ value that ensures that
the average ROM kinetic energy, $KE^{ROM}$, is the closest to the everage FOM
(LES) kinetic energy, $KE^{LES}$. Specifically, we solve the following
optimization problem:

$\displaystyle\min_{\gamma}\left|\overline{KE}^{ROM}-\overline{KE}^{LES}\right|.$
(34)

| $r$ | 4 | 8 | 16 | 32 | 40 | 50
---|---|---|---|---|---|---|---
EFR-ROM1 | $\gamma$ | 1.01e0 | 9.98e0 | 9.03e-1 | 8.05e-1 | 8.44e-1 | 7.00e-1
EFR-ROM2 | $\gamma$ | 2.90e-2 | 3.30e-2 | 3.90e-2 | 5.50e-2 | 7.08e-2 | 7.50e-2
Table 5:  EFR-ROM optimal $\gamma$ values for $\chi=6\times 10^{-3}$ and
different $r$ values.

The results in Table 5 display a relatively low sensitivity of the optimal
EFR-ROM parameter $\gamma$ with respect to changes in $r$. Indeed, as $r$
varies, the order of magnitude of the optimal $\gamma$ remains the same for
both EFR-ROM1 and EFR-ROM2, although the latter is more sensitive than the
former.

One possible explanation for the relatively low sensitivity of the optimal
EFR-ROM parameter $\gamma$ is that the EFR-ROM parameter $\chi$, which
controls the amount of filtering in the relaxation step of the EFR-ROM
algorithm, is low ($\chi=6\times 10^{-3}$). Thus, only $0.6\%$ filtering is
applied at each time step of the EFR-ROM algorithm. Since only a low amount of
filtering is used, the effect of the filter radius (i.e., the ROM lengthscales
$\delta_{1}$ and $\delta_{2}$) is not as important in the EFR-ROM case as in
the ML-ROM case.

| $r$ | 4 | 8 | 16 | 32 | 40 | 50
---|---|---|---|---|---|---|---
EFR-ROM1 | $\gamma$ | 2.92e-1 | 6.25e-2 | 6.00e-2 | 4.84e-2 | 4.41e-2 | 3.61e-2
EFR-ROM2 | $\gamma$ | 2.36e-4 | 6.83e-5 | 1.12e-04 | 2.26e-4 | 3.10e-4 | 4.13e-4
Table 6:  EFR-ROM optimal $\gamma$ values for $\chi=6\times 10^{-2}$ and
different $r$ values.

To investigate whether a higher percentage of filtering yields a higher
sensitivity of the optimal EFR-ROM parameter $\gamma$, we increase the $\chi$
value. Specifically, we choose $\chi=6\times 10^{-2}$, i.e., $6\%$ filtering
at each EFR-ROM time step (Table 6). The results in Table 6 show that the
optimal EFR-ROM $\gamma$ value is very sensitive with respect to changes in
$r$. Indeed, as $r$ increases from $4$ to $50$, $\gamma$ decreases by almost
one order of magnitude. The optimal EFR-ROM2 $\gamma$ values are less
sensitive with respect to changes in $r$ than the optimal EFR-ROM1 $\gamma$
values. Indeed, although the EFR-ROM2 $\gamma$ values vary, their order of
magnitude generally stays constant.

Overall, the results in Tables 5–6 show that the EFR-ROM’s parameter $\gamma$
sensitivity is higher for EFR-ROM1 than for EFR-ROM2 when a high percentage of
filtering is used in the EFR-ROM algorithm. As expected, when a low percentage
of filtering is used, both EFR-ROM1 and EFR-ROM2 display a relatively low
$\gamma$ sensitivity.

Based on the results in Figures 6–9 and in Tables 5–6, we conclude that the
EFR-ROM investigation in this section yields qualitative results that are
similar to those yielded by the ML-ROM investigation in Section 4.3.

##  5 Conclusions

In this paper, we proposed a novel ROM lengthscale definition. This new ROM
lengthscale, denoted $\delta_{2}$, was constructed by using energy
distribution arguments. Specifically, we balanced the ROM and FOM energy
content with the energy content in the $\delta_{2}$ and $h$ scales,
respectively, where $h$ is the FOM mesh size. We emphasize that the novel ROM
lengthscale, $\delta_{2}$, is fundamentally different from the current ROM
lengthscales, which are built by using dimensional arguments.

We compared the new ROM lengthscale, $\delta_{2}$, with a standard dimensional
based ROM lengthscale, denoted $\delta_{1}$. To this end, we used these two
ROM lengthscales to build two mixing-length ROMs (ML-ROMs) and two evolve-
filter-relax ROMs (EFR-ROMs) in which all the other parameters were the same.
We investigated the four resulting ML-ROMs and EFR-ROMs in the numerical
simulation of the turbulent channel flow at $Re_{\tau}=395$.

The numerical results of our investigation yielded the following conclusions:

  1. 1.

The new energy-based ROM lengthscale, $\delta_{2}$, was signficantly (two
orders of magnitude) larger than the standard ROM lengthscale, $\delta_{1}$.
As a result, the ML-ROMs and EFR-ROMs based on the new ROM lengthscale were
significantly more stable than the ML-ROMs and EFR-ROMs based on the standard
ROM lengthscale.

  2. 2.

The new energy-based ROM lengthscale displayed the correct asymptotic behavior
with respect to the ROM dimension, whereas the standard dimensionality-based
ROM lengthscale did not.

  3. 3.

The ML-ROM parameters based on the new energy-based ROM lengthscale were less
sensitive (i.e., more robust) with respect to changes in the ROM dimension
($r$) than the ML-ROMs parameters based on the standard dimensionality-based
ROM lengthscale. The EFR-ROM parameters based on the new lengthscale were less
sensitive with respect to changes in $r$ than the EFR-ROMs parameters based on
the standard lengthscale when a significant percentage of filtering was
performed. For a low filtering percentage, as expected, the EFR-ROM parameters
based on the two lengthscales displayed a relatively low sensitivity. In this
setting, the standard lengthscale yielded less sensitive parameters than the
new lengthscale.

The numerical assessment of the new energy-based ROM lengthscale yielded
encouraging results. We plan to further investigate this ROM lengthscale in
the construction of other types of ROMs, e.g., large eddy simulation ROMs [13,
38]. We also plan to leverage the new energy based lengthscale to develop
scale-aware ROM strategies that are better suited for flow-specific
applications.

## Acknowledgments

The work of the first and fourth authors was supported by NSF through grant
DMS-2012253 and CDS&E-MSS-1953113. The third author gratefully acknowledges
the U.S. DOE Early Career Research Program support through grant DE-SC0019290
and the NSF support through grant DMS-2012255. Part of this work was funded
under the nuclear energy advanced modeling and simulation program.

## References

  * [1] E. Merzari, A. Obabko, P. Fischer, N. Halford, J. Walker, A. Siegel, Y. Yu, Large-scale large eddy simulation of nuclear reactor flows: Issues and perspectives, Nucl. Eng. Des. 312 (2017) 86–98.
  * [2] L. Fick, Y. Maday, A. T. Patera, T. Taddei, A stabilized POD model for turbulent flows over a range of Reynolds numbers: Optimal parameter sampling and constrained projection, J. Comput. Phys. 371 (2018) 214–243.
  * [3] Y. Liu, R. Hu, L. Zou, G. Hu, D. Nunez, SAM-ML: Machine learning enhanced system analysis module for thermal stratification analysis.
  * [4] L. Vergari, A. Cammi, S. Lorenzi, Reduced order modeling approach for parametrized thermal-hydraulics problems: inclusion of the energy equation in the POD-FV-ROM method, Progress in Nuclear Energy 118 (2020) 103071.
  * [5] L. Vergari, A. Cammi, S. Lorenzi, Reduced order modeling for coupled thermal-hydraulics and reactor physics problems, Progress in Nuclear Energy 140 (2021) 103899.
  * [6] E. Merzari, W. D. Pointer, P. Fischer, A POD-based solver for the advection-diffusion equation, in: Fluids Engineering Division Summer Meeting, Vol. 44403, 2011, pp. 1139–1147.
  * [7] P. Holmes, J. L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge, 1996.
  * [8] S. E. Ahmed, S. Pawar, O. San, A. Rasheed, T. Iliescu, B. R. Noack, On closures for reduced order models $-$ a spectrum of first-principle to machine-learned avenues, Phys. Fluids 33 (9) (2021) 091301.
  * [9] S. Volkwein, Proper orthogonal decomposition: Theory and reduced-order modelling, Lecture Notes, University of Konstanzhttp://www.math.uni-konstanz.de/numerik/personen/volkwein/teaching/POD-Book.pdf (2013).
  * [10] J. S. Hesthaven, G. Rozza, B. Stamm, Certified Reduced Basis Methods for Parametrized Partial Differential Equations, Springer, 2015.
  * [11] A. Quarteroni, A. Manzoni, F. Negri, Reduced Basis Methods for Partial Differential Equations: An Introduction, Vol. 92, Springer, 2015.
  * [12] N. Aubry, P. Holmes, J. L. Lumley, E. Stone, The dynamics of coherent structures in the wall region of a turbulent boundary layer, J. Fluid Mech. 192 (1988) 115–173.
  * [13] Z. Wang, I. Akhtar, J. Borggaard, T. Iliescu, Proper orthogonal decomposition closure models for turbulent flows: A numerical comparison, Comput. Meth. Appl. Mech. Eng. 237-240 (2012) 10–26.
  * [14] D. Wells, Z. Wang, X. Xie, T. Iliescu, An evolve-then-filter regularized reduced order model for convection-dominated flows, Int. J. Num. Meth. Fluids 84 (2017) 598––615.
  * [15] M. Gunzburger, T. Iliescu, M. Mohebujjaman, M. Schneier, An evolve-filter-relax stabilized reduced order stochastic collocation method for the time-dependent Navier-Stokes equations, SIAM-ASA J. Uncertain. 7 (4) (2019) 1162–1184.
  * [16] M. Girfoglio, A. Quaini, G. Rozza, A POD-Galerkin reduced order model for a LES filtering approach, J. Comput. Phys. 436 (2021) 110260.
  * [17] M. Girfoglio, A. Quaini, G. Rozza, Pressure stabilization strategies for a LES filtering reduced order model, Fluids 6 (9) (2021) 302.
  * [18] M. Strazzullo, M. Girfoglio, F. Ballarin, T. Iliescu, G. Rozza, Consistency of the full and reduced order models for evolve-filter-relax regularization of convection-dominated, marginally-resolved flows, Int. J. Num. Meth. Eng. 123 (14) (2022) 3148–3178.
  * [19] F. Sabetghadam, A. Jafarpour, $\alpha$ regularization of the POD-Galerkin dynamical systems of the Kuramoto–Sivashinsky equation, Appl. Math. Comput. 218 (10) (2012) 6012–6026.
  * [20] K. Kaneko, P.-H. Tsai, P. Fischer, Towards model order reduction for fluid-thermal analysis, Nucl. Eng. Des. 370 (2020) 110866.
  * [21] M. Gunzburger, T. Iliescu, M. Schneier, A Leray regularized ensemble-proper orthogonal decomposition method for parameterized convection-dominated flows, IMA J. Numer. Anal. 40 (2) (2020) 886–913.
  * [22] T. Iliescu, H. Liu, X. Xie, Regularized reduced order models for a stochastic Burgers equation, Int. J. Numer. Anal. Mod. 15 (4–5) (2018) 594–607.
  * [23] S. Aradag, S. Siegel, J. Seidel, K. Cohen, T. McLaughlin, Filtered POD-based low-dimensional modeling of the 3D turbulent flow behind a circular cylinder, Int. J. Num. Meth. Fluids 66 (1) (2011) 1–16.
  * [24] I. Farcas, R. Munipalli, K. E. Willcox, On filtering in non-intrusive data-driven reduced-order modeling, in: AIAA AVIATION 2022 Forum, 2022, p. 3487\.
  * [25] C. Mou, E. Merzari, O. San, T. Iliescu, A numerical investigation of the lengthscale in the mixing-length reduced order model of the turbulent channel flow, in: 19th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-19), Brussels, Belgium, 2022.
  * [26] S. L. Brunton, J. N. Kutz, Data-driven science and engineering: Machine learning, dynamical systems, and control, Cambridge University Press, 2019.
  * [27] L. C. Berselli, T. Iliescu, W. J. Layton, Mathematics of Large Eddy Simulation of Turbulent Flows, Scientific Computation, Springer-Verlag, Berlin, 2006.
  * [28] P. Sagaut, Large Eddy Simulation for Incompressible Flows, 3rd Edition, Scientific Computation, Springer-Verlag, Berlin, 2006.
  * [29] M. Couplet, P. Sagaut, C. Basdevant, Intermodal energy transfers in a proper orthogonal decomposition–Galerkin representation of a turbulent separated flow, J. Fluid Mech. 491 (2003) 275–284.
  * [30] P. F. Fischer, J. Mullen, Filter-based stabilization of spectral element methods, C. R. Acad. Sci. Paris Sér. I Math. 332 (3) (2001) 265–270.
  * [31] W. J. Layton, L. G. Rebholz, Approximate Deconvolution Models of Turbulence: Analysis, Phenomenology and Numerical Analysis, Vol. 2042, Springer Berlin Heidelberg, 2012.
  * [32] J. S. Mullen, P. F. Fischer, Filtering techniques for complex geometry fluid flows, Commun. Numer. Meth. Engng. 15 (1) (1999) 9–18.
  * [33] V. J. Ervin, W. J. Layton, M. Neda, Numerical analysis of filter-based stabilization for evolution equations, SIAM J. Numer. Anal. 50 (5) (2012) 2307–2335.
  * [34] L. Bertagna, A. Quaini, A. Veneziani, Deconvolution-based nonlinear filtering for incompressible flows at moderately large Reynolds numbers, Int. J. Num. Meth. Fluids 81 (8) (2016) 463–488.
  * [35] S. Pope, Turbulent flows, Cambridge University Press, Cambridge, 2000.
  * [36] L. Rebholz, C. Zerfas, K. Zhao, Global in time analysis and sensitivity analysis for the reduced NS-$\alpha$ model of incompressible flow, J. Math. Fluid Mech. 19 (3) (2017) 445–467.
  * [37] L. Rebholz, T. Kim, Y. Byon, On an accurate $\alpha$ model for coarse mesh turbulent channel flow simulation, Appl. Math. Model. 43 (2017) 139–154.
  * [38] X. Xie, D. Wells, Z. Wang, T. Iliescu, Approximate deconvolution reduced order modeling, Comput. Methods Appl. Mech. Engrg. 313 (2017) 512–534.


L PREPRINT - accepted by 21st IEEE Interregional NEWCAS Conference 2023
(NEWCAS 2023)3.5cm,1cm1.001.8

# Finite State Automata Design using 1T1R ReRAM Crossbar

Simranjeet Singh15, Omar Ghazal2, Chandan Kumar Jha3, Vikas Rana5, Rolf
Drechsler34,
Rishad Shafik2, Alex Yakovlev2, Sachin Patkar1, Farhad Merchant2 {simranjeet,
<EMAIL_ADDRESS>{si.singh<EMAIL_ADDRESS><EMAIL_ADDRESS>
{O.G.G.Awf2, rishad.shafik, alex.yakovlev<EMAIL_ADDRESS>1Indian Institute of Technology Bombay, India, 2Newcastle University, UK,
3University of Bremen, Germany,
4DFKI GmbH, Germany, 5Forschungszentrum Jülich GmbH, Germany

###### Abstract

Data movement costs constitute a significant bottleneck in modern machine
learning (ML) systems. When combined with the computational complexity of
algorithms, such as neural networks, designing hardware accelerators with low
energy footprint remains challenging. Finite state automata (FSA) constitute a
type of computation model used as a low-complexity learning unit in ML
systems. The implementation of FSA consists of a number of memory states.
However, FSA can be in one of the states at a given time. It switches to
another state based on the present state and input to the FSA. Due to its
natural synergy with memory, it is a promising candidate for in-memory
computing for reduced data movement costs. This work focuses on a novel FSA
implementation using resistive RAM (ReRAM) for state storage in series with a
CMOS transistor for biasing controls. We propose using multi-level ReRAM
technology capable of transitioning between states depending on bias pulse
amplitude and duration. We use an asynchronous control circuit for writing
each ReRAM-transistor cell for the on-demand switching of the FSA. We
investigate the impact of the device-to-device and cycle-to-cycle variations
on the cell and show that FSA transitions can be seamlessly achieved without
degradation of performance. Through extensive experimental evaluation, we
demonstrate the implementation of FSA on 1T1R ReRAM crossbar.

###### Index Terms:

FSA, Machine Learning, ReRAM, Memristors, In-Memory Computing

##  I Introduction

In-memory computing (IMC) using memristive devices has become popular in
elevating the von Neumann bottleneck by storing and processing data in memory,
especially for machine learning (ML) applications [1, 2, 3, 4, 5]. Memristive
devices connected in a crossbar structure allow the program to run in
parallel, making the IMC comparable with conventional computing in terms of
energy efficiency and performance [6]. Memristive devices, such as resistive
random access memory (ReRAM) [7], can be configured as multi-level cell [8],
where the device has multiple intermediate states between low resistive (LRS)
and high resistive state (HRS). Still, modern ML workloads require massive
storage resources and parallelism to accelerate. However, finite state
automata (FSA) capture the real-world constraint of finite memory to implement
learning applications [9]. The multi-level behavior of ReRAM can be used to
design FSA in memory.

FSA is an abstract machine that can be at exactly one of the finite number of
states at any given time. FSA changes its state from one to another, called a
transition, when a specific event occurs [10]. FSA has a significant edge on
applications that require low-latency or real-time processing, such as
automated verification [11], autonomous vehicles [12], and tsetlin machine
[13]. The energy efficiency, high density, and IMC properties of ReRAM devices
make them suitable for FSA implementation [14] [15]. Multiple states of ReRAM,
between LRS and HRS, can be mapped to the states of FSA. An FSA cell can be
represented as a single ReRAM device with a CMOS transistor in series (1T1R),
as shown in Fig. 1(a). The material stack of the memristive cell is shown in
Fig. 1(b), along with the I-V characteristics of the 1T1R cell in Fig. 1(c).
Fig. 1(c) shows the multiple states in SET (switching to LRS) and RESET
(switching to HRS) states, which have been utilized in this study to implement
the FSA on ReRAM devices [16].

However, transitions of FSA from one state to another with accurate detection
of the current FSA state under device variations remains challenging. In this
paper, we propose an architecture to implement the FSA on the ReRAM crossbar
for IMC. We show in the proposed architectures how a single 1T1R cell can be
used to design a six-state FSA. Next, we evaluate the proposed architecture in
terms of energy efficiency and performance. Moreover, we assess the
architecture under device variations such as device-to-device (D2D) and cycle-
to-cycle (C2C) [17]. To summarize the main contributions:

  * •

The integration of 1T1R ReRAM technology into FSA design by utilizing the
multi-level behavior of ReRAM. The gradual RESET method has been utilized to
achieve the multi-level behavior.

  * •

Investigation of the impact of D2D and C2C variations on state transitions and
detection.

  * •

Extensive evaluation of the efficiency of the proposed architecture in terms
of energy efficiency and latency.

Figure 1: 1T1R cell for FSA (a) cell structure, (b) device material stack, and
(c) multi-level characteristic of a device.

The remainder of the paper is organized as follows: Section II presents the
proposed architecture to design FSA on ReRAM. Section III presents
experimental results and validation for the design. Finally, Section IV
concludes the paper.

##  II Proposed Architecture

This section discusses the architecture to implement the FSA on the 1T1R ReRAM
crossbar, depicted in Fig. 2. At the core of the proposed architecture, 1T1R
cells connected in a crossbar structure called finite automaton (FA) are used.

Figure 2: Architecture to train the FSA using 1T1R cell ($FA_{mn}$), where ‘m’
and ‘n’ represent rows and columns in a crossbar, respectively, and ‘p’ is the
number of ADC bits, which is given as $\lceil\log_{2}(s)\rceil$, ‘s’ is the
number of states in a FA cell.

###  II-A FA using 1T1R cell

The FA cell is constructed using one memristor and one NMOS transistor in
series. The structure of a single FA cell is shown in Fig. 1.The multi-level
characteristics of the FA, which have been mapped to different states, are
shown in the I-V characteristics displayed in Fig. 1(c). In this study, the FA
has seven states (‘s’) from $S_{0}$ to $S_{6}$, starting from LRS to HRS.
$S_{0}$ has a minimum resistance of around $7.8$K$\Omega$ and $S_{6}$ has
maximum resistance of around $1.5$M$\Omega$. All other states are mapped into
the intermediate values between $S_{0}$ to $S_{6}$. Each FA in the crossbar
represents six states (excluding $S_{0}$), and they can be independently
programmed. However, multiple FAs can be combined to run a complex application
that needs more than six states. For this work, we will limit our study to the
working of independent FA in the crossbar. The state transitions of FA are
examined next.

###  II-B State transitions in FA

An FA has a finite number of states, seven in this study, and it changes the
state from one to another or the next state ($q_{i+1}$) based on the input
($i_{i}$) and current state ($q_{i}$) similar to a mealy machine. A pulse
generation module in Fig. 2 can generate the different pulse widths of a fixed
voltage amplitude. The control circuit selects the appropriate signal for the
next state based on the present state and input. The parameter to switch the
state from $S_{0}$ to any possible state is given in Table I. Next, the analog
demultiplexer (DeMUX) and bit-line encoder select a FA in the crossbar by
applying an ON voltage to the NMOS transistor and the required pulse signal at
the row of the crossbar.

For the functional correctness of the FSA transition, it is important to
identify the present state correctly. FSA can jump from the present state to
any other possible state in FA. So, it is expected from the FSA that it should
give the same current value for the same state transition. However, the
gradual RESET method limits the transitions of the states only in the forward
direction ($S_{1}\rightarrow S_{6}$). In order to change the state which is
less than the current state (backward direction), FA needs to switch to
$S_{0}$ (intermediate state) before switching to the next desired state. Also,
it is expected that states after switching can correctly be identified during
forward or backward direction switching. Therefore, an intermediate state is
added in every state transition in FA, which provides three main advantages;
(a) switching to any state in FA, (b) state retention while looping in the
same state, and (c) reducing the complexity of control circuitry.

Figure 3: State transitions in a FA cell (a) and (b) shows the use of FSA as
Krinsky learning automaton [18].

The state transition graph of single FA is shown in Fig. 3(a). It can be used
for learning applications such as Krinsky automaton [18]. The mapping of a
Krinsky learning automaton has been shown in Fig. 3(b). The unfavorable
response ($\beta=1$) gradually moves toward the boundary states separating the
two actions, behaving as binary states. In another response ($\beta=0$),
$S_{i}$ switches to $S_{1}$ and $S_{4}$ for $(1\leq i\leq 3)$ and $(4\leq
i\leq 6)$, respectively. The proposed architecture is highly flexible
regarding its control circuit and switching characteristics. It provides the
facility to transit from the current state to any next state via $S_{0}$ with
an adaptive STG control unit. The proposed approach can accommodate FSM with
more than 6 states by utilizing multiple 1T1R cells arranged to represent
different states and encoded into binary form, offering flexibility for
varying numbers of states.

TABLE I: State transition of 1T1R cell $\rm V_{fixed}$ = 1.8V, $\rm V_{SET}$=
$-2V$, $\rm V_{READ}$= $0.1V$
---
State |  $\rm P_{width}$ ($n$s) |  $I_{d}$ ($\mu$A) | Resistance (K$\Omega$)
S0 | 10$n$s at -2V | 12.8 | 7.8
S1 | 5$n$s | 12.6 | 8.0
S2 | 10$n$s | 1.6 | 95.2
S3 | 15$n$s | 0.56 | 196.1
S4 | 30$n$s | 0.3 | 342.5
S5 | 60$n$s | 0.2 | 588.2
S6 | 150$n$s | 0.07 | 1492.5

###  II-C Peripherals to control FSA

Various peripherals around the crossbar are required to implement FSA on the
ReRAM crossbar, as shown in Fig. 2. The control unit decided the transitions
in FA, which are functions of $x_{i}$ and $q_{i}$.

_Pulse generation module_ generates the voltage pulses with a required
duration for transitioning the state from one to another according to Table I.
Since state $S_{0}$ is an intermediate state, the small gap between $S_{0}$
and $S_{1}$ will not pose an issue in state estimation.

_Multiplexer (MUX), demultiplexer (DeMUX), and bit-line encoder_ are the
selection peripherals. For a given ($m\times n$) size crossbar, ($1\times m$)
sized DeMUX are attached to select a row of the crossbar to apply a pulse for
transition. Bit-line encoder enables the transistor of selected FA. At column,
($n\times 1$) MUX is connected to read the state of a FA.

_Current sense amplifier (CSA) and Analog-to-digital converter (ADC)_ are
sensing peripherals. The CSA converts the current to an amplified voltage,
which is further used by the ADC to detect the current state of FA. A common
CSA and ADC have been used in the proposed architecture, where a FA can be
read in each cycle. A ‘p’ bit ADC is required to correctly detect the ‘s’
number of states in FA.

_Control unit_ contains the algorithms to switch the states of FA in sequence.
The control unit manipulates the select lines of MUX and DeMUX to select an FA
and generate the required pulses for transitions via the pulse generation
module. It takes input from the digital interface ($x_{i}$) and ADC ($q_{i}$)
and calculates the required control signal to switch the state to $q_{i+1}$;
alternatively, $q_{i+1}(q_{i},x_{i})$, where $0<=|x_{i}|<=1$ and
$S_{1}<=q_{i}<=S_{6}$.

##  III Experimental Results

This section evaluates the proposed methodology in terms of energy efficiency
and area. First, we study the switching characteristics of a FA and the state
transitions from one to another. Next, we look at the impact of D2D and C2C
variations on state transitions.

Figure 4: Multi-states behavior using gradual RESET method.

###  III-A 1T1R cell characterization

1T1R cell used for this study is designed using a Pt/Ti/TiOx/HfO2/Pt material
stack memristive devices in series with a 45nm transistor. The material stack
used in the memristive device adheres to the characteristics of the
experimental devices. Fig. 1 shows the configuration of the 1T1R cell along
with the material stack of the device and I-V characteristics. Table II shows
the parameters used for the ReRAM device model. The device has multiple-state
characteristics, and the gradual RESET method has been used to achieve this
behavior. In the gradual RESET method, the device is initialized to an LRS
state by applying a positive voltage of a specific duration on the ohmic
electrode (OE). Next, the device is switched into HRS by applying RESET pulses
gradually. Fig. 4 shows the gradual RESET method, where the voltages across
the device and its corresponding resistance states have been plotted. The
gradual RESET method results in multi-level behavior of 1T1R cells.
Additionally, a transistor in series helps to control the current precisely.

TABLE II: Model parameters Symbol | Value | Symbol | Value
---|---|---|---
$\rm l_{cell}$ | 3 $nm$ | $\rm l_{det}$ | 4 $nm$
$\rm r_{det}$ | 20 $nm$ | $\rm N_{plug}$ | $20\times 10^{26}\ m^{-3}$
$a$ | 0.25 $nm$ | $\mu_{n}$ | $1\times 10^{-6}\ m^{2}/{Vs}$
$\epsilon$ | 17 $\epsilon_{0}$ | $\rm N_{disc,min}$ | $0.008\times 10^{26}\ m^{-3}$
$\epsilon\phi\beta$ | 5.5 $\epsilon_{0}$ | $\rm N_{disc,max}$ | $20\times 10^{26}\ m^{-3}$
$e\phi_{\beta n0}$ | 0.3 $eV$ | $e\phi_{\beta n}$ | 0.1 $eV$
$\Delta\rm W_{A}$ | 0.7 $eV$ | A | 0.00392 1/$\Omega$
$\rm R_{series}$ | 650 $\Omega$ | $\rm R_{o}$ | 719.244 $\Omega$
$\rm R_{th,line}$ | 90.47 $K\Omega$ | $\rm R_{th0}$ |  $1.572\times 10^{7}$ $\Omega$
Figure 5: State switching from $S_{0}$ to $S_{6}$. For the desired state, the
device is first switched to $S_{0}$ and then directly switched to the required
state by applying the appropriate pulse.

###  III-B Realization of state transitions

As the memristive devices change their states based on the value of integral
over time of the applied voltage, the width of the applied pulse can be used
for state transitions instead of gradually varying voltage. Table I shows the
switching of FA from $S_{0}$-$S_{6}$ with different pulse widths. The state
$S_{0}$ is the initial state or LRS state, which can be switched by applying a
10ns pulse of -2V voltage. The states from $S_{1}$ to $S_{6}$ (6 states) are
the actual states used for FA. Generally, the FA switches the state in the
forward ($S_{1}$-$S_{6}$) or backward ($S_{6}$-$S_{1}$) directions. However,
the gradual method only works in forward state switching. The backward state
switching is tackled by switching to the intermediate state ($S_{0}$) before
switching to the desired state. For example, the current state is $S{3}$, and
the next expected state is $S_{2}$. In this case, the device first is switched
to $S_{0}$ and then switched from $S_{0}$ to $S_{2}$ ($S_{3}\rightarrow
S_{0}\rightarrow S_{2}$). To reduce the complexity of the control circuit and
accurate state detection, the intermediate state has been used in forward, and
backward switching. Fig. 5 shows the switching of the states through the
intermediate state. Table I also shows the states’ current and resistance,
indicating enough margin between the state’s current/resistance for accurate
state detection.

###  III-C Impact of variations

The D2D and C2C variation in ReRAM devices can affect the switching behavior.
To simulate D2D variations, the random set of values for device parameters
such as radius, length, and minimum and maximum oxygen vacancy in the disc are
drawn from the experiment-verified Gaussian distribution [17]. These
variations were then independently applied to the available devices in the
crossbar. C2C variations are simulated by changing the variable parameters in
a period of a single cycle. It has been observed that states from $S_{0}$ to
$S_{3}$ (low states) have a larger impact of the variations compared to high
states ($S_{4}$ to $S_{6}$), which is around $\pm$50% change in the read
current for low states and $\pm$20% for high states. However, for low states,
the margin between the state is more than five times which enables the
accurate detection of the state even if the variations have a larger impact.
Moreover, switching through the intermediate state prevent error accumulation
over time and reduces the impact of variations.

Figure 6: STG of the control unit, (a) state transition cycle, and (b) reading the present state in FA.  TABLE III: 1T1R cell energy consumption State Switching |  Intermediate State |  Energy ($pJ$) |  State Switching |  Intermediate State |  Energy ($pJ$)
---|---|---|---|---|---
$S_{0}\rightarrow S_{1}$ | – | 1.74 | $S_{1}\rightarrow S_{2}$ | $S_{0}$ | 8.2
$S_{2}\rightarrow S_{3}$ | $S_{0}$ | 8.3 | $S_{3}\rightarrow S_{4}$ | $S_{0}$ | 8.5
$S_{4}\rightarrow S_{5}$ | $S_{0}$ | 8.8 | $S_{5}\rightarrow S_{6}$ | $S_{0}$ | 9.25
Average energy | 7.5pJ

###  III-D Control circuitry

An asynchronous digital controller has been designed using Workcraft [19] to
coordinate data flow between different components in the architecture. Faster
operation and lower power consumption are some of the advantages of
asynchronous circuits over global clocked-based circuits. The control unit’s
signal transition graphs (STG) [20] to perform state detection (Read cycle)
and state transitions (write cycle) are shown in Fig. 6. When there is a
request to read the state of FA, the control unit receives a data reading
strobe (_DR+_) from the digital interface environment, and a reading cycle
begins. It activates MUX and bit line encoder to select a device (MEN+). Next,
it starts reading the data for the $n_{th}$ FA (DN+) after receiving
acknowledgment from MUX (MACK+). Lastly, the data is read, and the next read
cycle is prepared by resetting DN-, MEN-, and MACK-. The final ACK signal is
sent out as an acknowledgment for the digital interface.

The control circuit for state transitions provides a facility to transition
the FA state from any present state to any other possible state. This
increases the flexibility of the proposed architecture to run any FSA
application. However, every state transition in FA is done via $S_{0}$ to
maintain functional correctness and reduce controller complexity. The control
circuitry handles this situation by generating an acknowledgment signal, which
includes the transition of $S_{0}$ and desired state ($S_{n}$). The transition
cycle is initiated whenever there is a request on a DW+ signal from the
environment. The first step for state transition includes switching to the
$S_{0}$ state. The FA is selected by enabling the row multiplexer and bit-line
encoder (MEN+). The FA changes the state from the $n_{th}$ state to $S_{0}$
and disables the row MUX before sending the acknowledgment for the final
transition. At this time of the cycle, FA is in $S_{0}$ and ready to be
switched to the desired state ($S_{n}$). Similar to $S_{0}$ switching, the
final transition starts by enabling the peripherals. Before the transition
cycle is finished, the signal DN- initiates the ACK+ that will be delivered to
the digital interface, resetting DW-, MEN-, and MACK-. The final
acknowledgment for the digital interface is ACK- signal, which indicates a
successful state transition.

###  III-E Energy and latency analysis

Each state transition in FA consumes different energy, which is given in Table
III. In every transition in FA, 7.5pJ energy is consumed on average. An FA has
to switch to the intermediate state before the desired transition, which
increases energy consumption. However, the intermediate state makes the state
switching robust against D2D and C2C variations. As the proposed architecture
used different pulse duration to state transitions, state $S_{6}$ takes 150ns
pulse. Hence the pulse generation module generates each pulse for a 150ns
period with varying widths. However, latency can be improved further by
increasing the voltage amplitude of the applied pulse, which increases energy
consumption. So, there is a trade-off between energy consumption and latency.

##  IV Conclusions

In this work, for the first time, we proposed the architecture to implement
the FSA using a 1T1R ReRAM crossbar. This paper offers insights into the scope
of FSA utilizing ReRAM and CMOS technology. We use the multi-level
characteristics of ReRAM, achieved using the gradual RESET method, to
implement FSA on the crossbar. We studied the impact of variation on state
transitions. Finally, we evaluated the proposed framework in terms of latency
and energy consumption. The results are encouraging and demonstrate the
potential for using ReRAM-based FSA designs. We will explore the prototyping
of the proposed designs and test the architecture with learning automaton
applications in the future.

## Acknowledgments

This work was supported in part by the Federal Ministry of Education and
Research (BMBF, Germany) in the project NEUROTEC II under Project 16ME0398K,
Project 16ME0399, German Research Foundation (DFG) within the Project PLiM (DR
287/35-1, DR 287/35-2) and through Dr. Suhas Pai Donation Fund at IIT Bombay.

## References

  * [1] G. Pedretti and D. Ielmini, “In-memory computing with resistive memory circuits: Status and outlook,” _Electronics (Switzerland)_ , vol. 10, 05 2021\.
  * [2] S. Yu, A. Soltan, R. Shafik, T. Bunnam, F. Xia, D. Balsamo, and A. Yakovlev, “Current-mode carry-free multiplier design using a memristor-transistor crossbar architecture,” in _2020 Design, Automation & Test in Europe Conference & Exhibition (DATE)_, 2020, pp. 638–641.
  * [3] C. K. Jha, P. L. Thangkhiew, K. Datta, and R. Drechsler, “Imagin: Library of imply and magic nor-based approximate adders for in-memory computing,” _IEEE Journal on Exploratory Solid-State Computational Devices and Circuits_ , vol. 8, no. 2, pp. 68–76, 2022.
  * [4] S. Froehlich and R. Drechsler, “Unlocking approximation for in-memory computing with cartesian genetic programming and computer algebra for arithmetic circuits,” _it - Information Technology_ , vol. 64, no. 3, pp. 99–107, 2022.
  * [5] F. Staudigl, F. Merchant, and R. Leupers, “A survey of neuromorphic computing-in-memory: Architectures, simulators, and security,” _IEEE Design & Test_, vol. 39, no. 2, pp. 90–99, 2022.
  * [6] A. Sebastian, M. Gallo, R. Khaddam-Aljameh, and E. Eleftheriou, “Memory devices and applications for in-memory computing,” _Nature Nanotechnology_ , vol. 15, 03 2020.
  * [7] R. Waser, R. Dittmann, G. Staikov, and K. Szot, “Redox‐based resistive switching memories – nanoionic mechanisms, prospects, and challenges,” _Advanced Materials_ , vol. 21, pp. 2632–2663, 07 2009.
  * [8] A. Siemon, D. Wouters, S. Hamdioui, and S. Menzel, “Memristive device modeling and circuit design exploration for computation-in-memory,” 05 2019, pp. 1–5.
  * [9] M. O. Rabin and D. Scott, “Finite automata and their decision problems,” _IBM Journal of Research and Development_ , vol. 3, no. 2, pp. 114–125, 1959\.
  * [10] A. Rezvanian, B. Moradabadi, M. Ghavipour, M. M. Daliri Khomami, and M. R. Meybodi, _Introduction to Learning Automata Models_. Cham: Springer International Publishing, 2019, pp. 1–49. [Online]. Available: https://doi.org/10.1007/978-3-030-10767-3_1
  * [11] K. Luo and Y. Liu, “Automatic verification of fsa strategies via counterexample-guided local search for invariants,” in _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , ser. IJCAI’19. AAAI Press, 2019, p. 1814–1821.
  * [12] M. Zhang, N. Li, A. Girard, and I. Kolmanovsky, “A finite state machine based automated driving controller and its stochastic optimization,” 10 2017, p. V002T07A002.
  * [13] O.-C. Granmo, “The Tsetlin Machine – A Game Theoretic Bandit Driven Approach to Optimal Pattern Recognition with Propositional Logic,” 2018. [Online]. Available: https://arxiv.org/abs/1804.01508
  * [14] A. Wheeldon, R. Shafik, T. Rahman, J. Lei, A. Yakovlev, and O.-C. Granmo, “Learning automata based energy-efficient AI hardware design for IoT applications,” _Philosophical Tran of the Royal Society A_ , vol. 378, no. 2182, p. 20190593, 2020.
  * [15] Y. Halawani, B. Mohammad, M. Abu Lebdeh, M. Al-Qutayri, and S. F. Al-Sarawi, “Reram-based in-memory computing for search engine and neural network applications,” _IEEE Journal on Emerging and Selected Topics in Circuits and Systems_ , vol. 9, no. 2, pp. 388–397, 2019.
  * [16] H. Aziza, S. Hamdioui, M. Fieback, M. Taouil, M. Moreau, P. Girard, A. Virazel, and K. Coulié, “Multi-level control of resistive ram (rram) using a write termination to achieve 4 bits/cell in high resistance state,” _Electronics_ , vol. 10, no. 18, 2021. [Online]. Available: https://www.mdpi.com/2079-9292/10/18/2222
  * [17] C. Bengel, A. Siemon, F. Cüppers, S. Hoffmann-Eifert, A. Hardtdegen, M. von Witzleben, L. Hellmich, R. Waser, and S. Menzel, “Variability-aware modeling of filamentary oxide-based bipolar resistive switching cells using SPICE level compact models,” _IEEE TSCAS I_ , vol. 67, no. 12, pp. 4618–4630, 2020\.
  * [18] R. Iraji, M. T. Manzuri-Shalmani, A. H. Jamalian, and H. Beigy, “Ija automaton: Expediency and $e$-optimality properties,” in _2006 5th IEEE International Conference on Cognitive Informatics_ , vol. 1, 2006, pp. 617–622.
  * [19] I. Poliakov, V. Khomenko, and A. Yakovlev, “Workcraft – a framework for interpreted graph models,” in _Applications and Theory of Petri Nets_ , G. Franceschinis and K. Wolf, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2009, pp. 333–342, (available from http://workcraft.org).
  * [20] L. Rosenblum and A. Yakovlev, “Signal Graphs: From Self-Timed to Timed Ones,” in _Int. Workshop on Timed Petri Nets_ , 1985, pp. 199–206.