# Fast Macroscopic Forcing Method

Spencer H. Bryngelson<EMAIL_ADDRESS>Florian Schäfer
<EMAIL_ADDRESS>Jessie Liu  Ali Mani  School of Computational
Science & Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA
Department of Mechanical Engineering, Stanford University, Stanford,
California 94305, USA

###### Abstract

The macroscopic forcing method (MFM) of Mani and Park [1] and similar methods
for obtaining turbulence closure operators, such as the Green’s function-based
approach of Hamba [2], recover reduced solution operators from repeated direct
numerical simulations (DNS). MFM has already been used to successfully
quantify Reynolds-averaged Navier–Stokes (RANS)-like operators for homogeneous
isotropic turbulence and turbulent channel flows. Standard algorithms for MFM
force each coarse-scale degree of freedom (i.e., degree of freedom in the RANS
space) and conduct a corresponding fine-scale simulation (i.e., DNS), which is
expensive. We combine this method with an approach recently proposed by
Schäfer and Owhadi [3] to recover elliptic integral operators from a
polylogarithmic number of matrix–vector products. The resulting Fast MFM
introduced in this work applies sparse reconstruction to expose local features
in the closure operator and reconstructs this coarse-grained differential
operator in only a few matrix–vector products and correspondingly, a few MFM
simulations. For flows with significant nonlocality, the algorithm first
“peels” long-range effects with dense matrix–vector products to expose a more
local operator. We demonstrate the algorithm’s performance for scalar
transport in a laminar channel flow and momentum transport in a turbulent
channel flow. For these problems, we recover eddy diffusivity- and eddy
viscosity-like operators, respectively, at $1\%$ of the cost of computing the
exact operator via a brute-force approach for the laminar channel flow problem
and $13\%$ for the turbulent one. We observe that we can reconstruct these
operators with an increase in accuracy by about a factor of $100$ over
randomized low-rank methods. Applying these operators to compute the averaged
fields of interest has visually indistinguishable behavior from the exact
solution. Our results show that a similar number of simulations are required
to reconstruct the operators to the same accuracy under grid refinement. Thus,
the accuracy corresponds to the physics of the problem, not the numerics. We
glean that for problems in which the RANS space is reducible to one dimension,
eddy diffusivity and eddy viscosity operators can be reconstructed with
reasonable accuracy using only a few simulations, regardless of simulation
resolution or degrees of freedom.

###### keywords:

Multi-scale modeling, eddy diffusivity, turbulence modeling, operator
recovery, numerical homogenization

##  1 Introduction

††
Code available at: https://github.com/comp-physics/fast-mfm

Well-established equations describe even the most complicated flow physics.
Still, full-resolution simulations of them stretch computational resources.
Reduced-complexity surrogate models are a successful approach to reducing
these costs. Historically, physical insight and analytical techniques have
been used to develop these models, including the RANS closure models [4].
However, data-driven approaches are emerging as semi-automated tools to
accomplish the same task. Some approaches attempt to represent the time
evolution of the physical system via neural networks, a formidable task that
involves reducing the entire Navier–Stokes operator [5, 6]. An alternative
approach is to compute effective equations that act on spatial or temporal
averages. The governing equations are projected into the reduced or averaged
space, and a forcing function is applied to examine the effect of the
underlying fluctuations on the averaged behavior. Kraichnan [7] and Hamba [2]
examined Green’s function solutions (i.e., using Dirac-delta-function-type
forcing) to scalar and momentum transport equations to develop exact
expressions for closure operators. Similarly, the macroscopic forcing method
(MFM) of Mani and Park [1] quantifies closure operators exactly by examining
forcing and averaged responses, called input–output pairs. However, as a
linear-algebra-based technique, MFM does not require the use Dirac delta
functions as forcing basis functions, and others like polynomials [8] and
harmonic functions [9], can be used.

MFM has been successfully applied to close reacting flow equations [10, 11]
and analyze homogeneous isotropic turbulence [9] and turbulent channel flow
[12]. MFM is analogous to numerical homogenization, or the finite-dimensional
approximation of solution spaces of partial differential equations (PDEs)
[13]. These techniques amount to operator recovery or learning, where an
unknown operator is estimated from a set of input–output pairs obtained from
full-resolution simulations. These simulations are computationally expensive,
so there is a pressing need to reduce the number of samples required, which we
address in this work.

Using MFM, one constructs effective operators, or macroscopic operators,
acting on solution averages from full-resolution simulations, called direct
numerical simulations (DNS) of the governing or microscopic equations [1]. If
the macroscopic operators are linear, the MFM procedure is no different from
estimating a matrix from a limited number of matrix–vector products. The
number of microscopic simulations required to recover the macroscopic operator
exactly equals the number of macroscopic degrees of freedom, which can be
prohibitively large for many simulation problems, like high-Re turbulence.

One can partially address this problem by working in Fourier space [1] or
fitting a parametric model to approximate the eddy diffusivity operator [14,
12]. However, the former requires spatial homogeneity, and the latter’s
accuracy depends on the parametric model’s quality. Liu et al. [8] introduces
an improved model that uses the nonlocal eddy diffusivity operator’s moments
to approximate the full operator. While these are viable approaches and the
subject of ongoing work, the target of this work is to reconstruct the full
discretely-defined nonlocal eddy diffusivity operator, as opposed to
prescribing or modeling its shape.

For many flows of practical interest, the nonlocal effects of closure terms
show diffusive behavior. Thus, work on operator recovery for elliptic PDEs is
closely related to MFM. [15] propose a “peeling” approach for recovering
hierarchical matrices from a polylogarithmic number of matrix–vector products,
although without rigorous bounds on the approximation error. Extensions of
this algorithm were proposed by [16, 17, 18]. In this setting,
eigendecompositions and randomized linear algebra have been used to recover
elliptic solution operators from matrix–vector products [19, 20, 21, 22].
Since the eigenvalues of elliptic solution operators follow a power law, these
methods require $\mathrm{poly}(1/\varepsilon)$ matrix–vector products to
obtain an $\varepsilon$ approximation of the operator. In contrast, Schäfer et
al. [23] showed that hidden sparsity of the solution operator results in an
$\varepsilon$ approximation from only $\mathrm{poly}(\log(1/\varepsilon))$
carefully crafted matrix–vector products. This speedup amounts to an
exponential reduction in the number of matrix–vector products. The sparsity
used by Schäfer and Owhadi [3] results from the locality of the partial
differential operator shared by local fluid models.

We use this approach to accelerate the MFM to create the Fast MFM. The Fast
MFM reveals the locality of the physical models to reduce the sample
complexity of standard MFM operator recovery. We apply Fast MFM to
inhomogeneous and turbulent problems and reconstruct the RANS closure
operators. Specifically, we consider passive scalar transport in a laminar 2D
channel flow following Mani and Park [1] and reconstruct the corresponding
eddy diffusivity operator, and momentum transport in a canonical turbulent 3D
channel flow at $\mbox{\text{Re}}_{\tau}=180$ following Hamba [24] and Park
and Mani [12] and reconstruct the corresponding eddy viscosity operator. These
examples display sufficient spatio-temporal richness in their dynamics to
argue that the Fast MFM can be applied more broadly. For example, one could
tackle the open closure problems associated with multiphase flows [25, 26, 27,
28], though we do not address such extensions here.

We briefly introduce MFM and similar approaches for recovering turbulence
closure operators in section 2. Section 3 details the mathematical foundations
of the sparse reconstruction procedure and section 3.5 applies it to MFM with
an extension to nonsymmetric operators, resulting in the Fast MFM. Results are
presented in section 4, focusing on the 2D and 3D problems analyzed by Mani
and Park [1] and Park and Mani [12]. Section 5 discusses the outlook of sparse
reconstruction methods like the one presented for other flow problems and PDEs
broadly.

##  2 Background on the Macroscopic Forcing Method (MFM)

###  2.1 The macroscopic forcing method

Given a set of linear microscopic equations,

$\displaystyle\mathcal{L}c=s,$ (1)

and an averaging operator $(s,c)\mapsto(\mathrm{bar}{s},\mathrm{bar}{c})$, the
macroscopic (averaged) operator $\overline{\mathcal{L}}$ is defined to
satisfy, for all microscopic solutions $c,s$ of (1),

$\displaystyle\overline{\mathcal{L}}\mathrm{bar}{c}=\mathrm{bar}{s}.$ (2)

Often, (1) are advection-diffusion equations for scalar transport or
linearizations of nonlinear PDEs such as Navier–Stokes equations, so the
averaging operation can be written as

$\displaystyle\mathrm{bar}{s}=\frac{1}{L_{2}\cdots
L_{N_{d}}}\int_{\Omega_{2}}\cdots\int_{\Omega_{N_{d}}}s(x_{1},\ldots
x_{N_{d}})\,\text{d}x_{2}\cdots\text{d}x_{N_{d}},$ (3)

where

$\displaystyle\Omega=\Omega_{1}\times\Omega_{2}\times\cdots\times\Omega_{N_{d}}$
(4)

is the physical (possibly spatio-temporal) domain and $L_{i}$ are the lengths
in each coordinate direction $x_{i}$, $i\in\\{1,\dots,N_{d}\\}$. In this
example, the averaged (2) is a univariate problem in the non-averaged
coordinate $x_{1}$. However, we point out that the techniques outlined in this
work apply to a wider range of possible averaging operations.

Figure 1:  Schematic of the MFM.

Using MFM, one can determine the exact linear operator
$\overline{\mathcal{L}}$ that acts on averages of flow statistics [1]. They
infer this operator by generating solution pairs to
$\overline{\mathcal{L}}\mathrm{bar}{c}=\mathrm{bar}{s}$, obtained from solving
the microscopic equations with forcing $\mathrm{bar}{s}$ and macroscopic
solution average $\mathrm{bar}{c}$.

Figure 1 show an example MFM procedure schematically for a two-dimensional
problem with coordinate directions $x_{1}$ and $x_{2}$. The relevant averaging
direction is $x_{2}$, with averaged “strips” indicating averaging. This
solution is forced by a field $s(x_{1},x_{2})$ as a Dirac delta function at a
specific $x_{1}$ coordinate equivalent to its averaged field
$\mathrm{bar}{s}(x_{1})$. The inverse solution operator $\mathcal{L}^{-1}$
solves the problem (1) for $c(x_{1},x_{2})$ given $\mathrm{bar}{s}(x_{1})$.
This is computationally equivalent to solving the full-resolution system (1)
or DNS. The averaged solution field $\mathrm{bar}{c}(x_{1})$ corresponds to a
column (“recover a column”) of a macroscopic solution operator
$\overline{\mathcal{L}}^{-1}$ under this averaging scheme. This procedure is
repeated for all non-averaged degrees of freedom.

In this example, the non-averaged coordinate is $x_{1}$, so each discretized
$x_{1,i}$ is locally forced (via a Dirac delta function) with
$s\equiv\mathrm{bar}{s}$. Completing the MFM procedure gives access to the
matrix representation $\overline{\mathcal{L}}$ via
$\bm{s}\mapsto\overline{\mathcal{L}}^{-1}\bm{s}$. Since evaluating this map
involves a high-resolution simulation, column-by-column construction of
$\overline{\mathcal{L}}^{-1}$ is intractable.

###  2.2 The linear algebra of MFM

A linear algebraic perspective is useful for understanding MFM. To this end,
we denote as $\mathcal{L}$ the matrix representation of a discretized
advection-diffusion operator

$\displaystyle\mathcal{L}=\frac{\partial}{\partial
t}+u\cdot\nabla-\nabla\cdot(a\nabla),$ (5)

where coefficients $a$, $u$, are allowed to vary in space and time. The
inverse operator, $\mathcal{L}^{-1}$, takes a spatio-temporal forcing term $s$
as input and returns the spatio-temporal field $c=\mathcal{L}^{-1}s$ by
solving the PDE. Let $\bm{P}$ denote a projection onto coarse-scale features
of interest, for example, spatio-temporal averages, and $\bm{E}$ denote an
extension such that $\bm{P}\bm{E}=\mathcal{I}$, where $\mathcal{I}$ is the
identity matrix. In the example of fig. 1, rows of $\bm{P}$ correspond to
averages of $c$ in the $x_{2}$-direction of the domain, and the rows of
$\bm{E}$ extend $\mathrm{bar}{s}$ to $s$. The macroscopic operator can then be
expressed as [1]

$\displaystyle\overline{\mathcal{L}}=\left(\bm{P}\mathcal{L}^{-1}\bm{E}\right)^{-1},$
(6)

where $\overline{\mathcal{L}}$ and $\mathcal{L}^{-1}$ are now discretized.

Another perspective on MFM can be obtained considering
$\overline{\mathcal{L}}$ in discretized form. By using bases for its row and
column space that consider the row and column spaces of $\bm{P}$ and $\bm{E}$,
we obtain a $2\times 2$ block matrix. After eliminating the second block, the
macroscopic operator is obtained as the Schur complement of $\mathcal{L}$.

$\displaystyle\overline{\mathcal{L}}=\left(\left(\mathcal{L}^{-1}\right)_{1,1}\right)^{-1}=\mathcal{L}_{1,1}-\mathcal{L}_{1,2}\left(\mathcal{L}_{2,2}\right)^{-1}\mathcal{L}_{2,1}.$
(7)

Computing $\mathcal{L}^{-1}$ or
$\left(\mathcal{L}_{2,2}\right)^{-1}\mathcal{L}_{2,1}$ naïvely, column by
column, requires as many solutions of the full-scale problem as there are
coarse scale degrees of freedom.

###  2.3 Inverse MFM

As shown in fig. 2, $\overline{\mathcal{L}}$ is more local than its computed
inverse. Hoping to turn this locality into computational gains, Mani and Park
[1] propose an inverse MFM to directly compute matrix–vector products with
$\overline{\mathcal{L}}$ without first having to compute
$\overline{\mathcal{L}}^{-1}$ using the ordinary MFM. This procedure can be
interpreted as evaluating the right-hand side of (7) at the cost of solving a
system of equations in $\mathcal{L}_{2,2}$. If (1) is an evolution PDE, this
procedure can be interpreted as a control problem, where at each time step,
the microscopic portion of the forcing $s$ is chosen to maintain a target
average $\mathrm{bar}{c}$. The resulting averaged forcing $\mathrm{bar}{s}$ is
the same as $\overline{\mathcal{L}}\mathrm{bar}{c}$.

Figure 2:  Matrices $\mathcal{D}$, $\overline{\mathcal{L}}$, and
$\overline{\mathcal{L}}^{-1}$ as labeled for the 2D channel flow case with 50
non-averaged grid points (and so each is $50\times 50$) for illustration
purposes. The inverse operator matrix of (c) is nearly dense, though (a) and
(b) are more strongly banded. Similar behavior is observed for finer
discretizations, which correspond to larger matrices.  Figure 3:  The middle
row of the matrices of fig. 2 on a log scale. $\mathcal{D}$ and
$\overline{\mathcal{L}}$ have a similar degree of locality, with entries
decreasing in magnitude algebraically from the diagonal. The discretized
operator $\overline{\mathcal{L}}^{-1}$ is dense; along its diagonal entries
only decay modestly at the boundaries.

###  2.4 The eddy diffusivity operator

As done in Mani and Park [1], consider a problem in which the averaging
operation includes averaging over the entire temporal domain and all
directions except $x_{1}$. Using a Reynolds decomposition, the velocity and
scalar fields can be decomposed as

$\displaystyle u=\mathrm{bar}{u}+u^{\prime}\quad\text{and}\quad
c=\mathrm{bar}{c}+c^{\prime},$ (8)

where $(\cdot)^{\prime}$ denotes fluctuations about the mean. Substitution of
(8) into the advection–diffusion equation (5) with $\mathcal{L}c=0$ and
averaging results in the corresponding mean scalar equation, we have

$\displaystyle\frac{\partial}{\partial
x_{1}}\overline{u_{1}^{\prime}c^{\prime}}-a\frac{\partial^{2}\mathrm{bar}{c}}{\partial
x_{1}^{2}}=0,$ (9)

where the scalar flux, $\overline{u_{1}^{\prime}c^{\prime}}$, is unclosed, and
(9) can be written as $\overline{\mathcal{L}}\mathrm{bar}{c}=0$. The
macroscopic operator, $\overline{\mathcal{L}}$, contains both the closure for
the scalar flux term and the closed molecular diffusion term. Thus,
$\overline{\mathcal{L}}$ can be further decomposed as

$\displaystyle\overline{\mathcal{L}}=-\left(\frac{\partial}{\partial
x_{1}}\right)\left(\mathcal{D}+a\mathcal{I}\right)\left(\frac{\partial}{\partial
x_{1}}\right),$ (10)

where $\mathcal{D}$ is the eddy diffusivity matrix. In continuous form, this
is equivalent to

$\displaystyle-\overline{u_{1}^{\prime}c^{\prime}}(x_{1})=\int_{y_{1}}\mathcal{D}(x_{1},y_{1})\frac{\partial\mathrm{bar}{c}}{\partial
x_{1}}{\bigg{|}}_{y_{1}}\text{d}y_{1},$ (11)

which generalizes to the nonlocal eddy diffusivity of Berkowicz and Prahm
[29]:

$\displaystyle-\overline{u_{\beta}^{\prime}c^{\prime}}(\bm{x})=\int_{\bm{y}}\mathcal{D}^{(\beta\alpha)}(\bm{x},\bm{y})\frac{\partial\mathrm{bar}{c}}{\partial
x_{\alpha}}{\bigg{|}}_{\bm{y}}\text{d}\bm{y},$ (12)

where $\alpha,\beta\in\\{1,\dots,N_{d}\\}$ are coordinate directions in the
macroscopic space.

As shown in fig. 2 and fig. 3, the eddy diffusivity matrix is significantly
more regular than $\overline{\mathcal{L}}$, making it a preferred target for
an operator recovery strategy. To leverage these properties in operator
recovery, Mani and Park [1] proposed a method for computing matrix–vector
products $\mathcal{D}\bm{a}$ as $\Sigma\overline{\mathcal{L}}\Sigma\bm{a}$,
where $\Sigma$ denotes the antiderivative and $\bm{a}$ is the vector state.
The objective of the present work is to recover $\mathcal{D}$, and thus
$\overline{\mathcal{L}}$, from as few matrix–vector products as possible, as
accurately as possible.

##  3 LU reconstruction of elliptic operators

###  3.1 Reconstructing elliptic operators from matrix–vector products

We use the LU variant of the _Cholesky reconstruction_ of Schäfer and Owhadi
[3] to construct the eddy diffusivity operator. Schäfer and Owhadi [3] prove
that the solution operators of divergence form elliptic partial differential
equations in dimension $N_{d}$ can be reconstructed to accuracy $\epsilon$
from only $\mathcal{O}\left(\log^{N_{d}+1}(\epsilon^{-1})\right)$ solutions
for carefully selected forcing terms. We briefly review this approach, which
forms the basis of this work.

###  3.2 Graph coloring

_Graph coloring_ allows one to reconstruct multiple columns of a sparse matrix
from a single matrix–vector product. The key idea is to identify groups of
columns with non-overlapping sparsity sets and use a right-hand side that only
activates those columns. As illustrated in fig. 4, the selected columns can be
read off from the resulting matrix–vector product. Similarly, graph coloring
can also reveal the leading columns of a sparse LU factorization. Once a row-
column pair of the LU factors are identified, it can be used to correct the
matrix–vector products to reveal later columns. This procedure, a variant of
which was first proposed by Lin et al. [15], is referred to as peeling (fig.
4).

###  3.3 Cholesky factors in wavelet basis

Figure 4: (a) Sparse recovery: Columns with non-overlapping but known sparsity
patterns (shown in color) can be recovered by a single matrix–vector product
via a carefully chosen vector. (b) Factorization: Cholesky factorizations with
leading-column sparsity patterns can be recovered similarly. (c) Peeling: If
denser columns of the factorization prevent recovery of sparser ones, identify
dense columns first and subtract their contribution to recover the sparser
ones. [This figure was adapted from Schäfer and Owhadi [3] with author
permission.]

It is well-known that the solution operators of elliptic PDEs are dense, owing
to the long-range interactions produced by diffusion. However, Schäfer et al.
[23] show that when represented in a multiresolution basis ordered _from
coarse to fine_ , solution operators of elliptic PDEs have almost sparse
Cholesky factors. This phenomenon is illustrated in fig. 5. The leading
columns of the Cholesky factors, corresponding to a coarse-scale basis
function with global support, are dense and therefore limit the efficiency of
graph coloring. However, they are few and can be identified efficiently and
removed via peeling. This procedure can be repeated to reveal progressively
finer columns. The growing number of basis functions on finer scales is
compensated for by their smaller support and, thus, increased gains due to
graph coloring. Thus, the number of matrix–vector products required is
approximately constant across levels. The resulting procedure is described in
algorithm 1. As described in Schäfer and Owhadi [3], the operation
$\texttt{scatter}_{c}$ takes in the vector obtained from a peeled
matrix–vector product computed in 5 and 6 and uses it to recover the columns
associated with the color $c$. In principle, further compression of the
resulting operator is possible using the techniques of Schäfer et al. [30].

1: $\mathbf{L}\leftarrow 0\times N$ empty matrix

2: $\mathbf{U}\leftarrow N\times 0$ empty matrix

3: $\mathbf{D}\leftarrow 0\times 0$ empty diagonal matrix

4: for $c\text{ color}$ do

5:
$\mathbf{L}_{\mathrm{new}}\leftarrow\texttt{scatter}_{c}\left(\mathbf{W}^{\top}\mathcal{D}\mathbf{M}_{:,c}-\mathbf{L}\mathbf{D}\mathbf{U}\mathbf{W}^{\top}\mathbf{M}_{:,c}\right)$

6:
$\mathbf{U}_{\mathrm{new}}\leftarrow\texttt{scatter}_{c}\left(\mathbf{W}^{\top}\mathcal{D}\mathbf{M}_{:,c}-\mathbf{U}^{\top}\mathbf{D}\mathbf{L}^{\top}\mathbf{W}^{\top}\mathbf{M}_{:,c}\right)$

7:
$\mathbf{D}_{\mathrm{new}}\leftarrow\operatorname{diag}\left(\left(\mathbf{L}_{\mathrm{new}}+\mathbf{U}_{\mathrm{new}}\right)/2\right)^{-1}$

8:
$\mathbf{L},\mathbf{U},\mathbf{D}\leftarrow\texttt{hcat}\left(\mathbf{L},\mathbf{L}_{\mathrm{new}}\right),\texttt{vcat}\left(\mathbf{U},\mathbf{U}_{\mathrm{new}}\right),\texttt{dcat}\left(\mathbf{D},\mathbf{D}_{\mathrm{new}}\right)$
{Concatenate to prior solution}

9: end for

10: $\mathbf{L}\leftarrow\mathbf{L}\mathbf{D}$

11: return approximate LU factorization $\mathbf{L}\mathbf{U}$ of
$\mathbf{W}^{\top}\mathcal{D}\mathbf{W}$

Algorithm 1 LU reconstruction in wavelet basis given by $\mathbf{W}$ with
measurements given by $\mathbf{M}$.

###  3.4 Adaptation to Fast MFM

The eddy diffusivity operator is not a divergence-form elliptic solution
operator. In particular, it is not symmetric. Instead of the Cholesky recovery
of Schäfer and Owhadi [3], we use an LU recovery that recovers a sparse LU
factorization of the target matrix. Columns of $\mathrm{L}$ are recovered from
matrix–vector products, and rows of $\mathrm{U}$ are recovered from matrix-
transpose–vector products (transpose–vector products), which can be computed
by solving the adjoint equation of $\mathcal{L}$. No rigorous guarantees exist
for the accuracy of LU reconstruction applied to eddy diffusivity matrices.
However, Schäfer et al. [23] show a wide range of diffusion-like operators,
including those produced by fractional-order Matérn or Cauchy kernels, produce
sparse Cholesky factors, despite the lack of theory supporting this
observation.

|  |
---|---|---
(a) | (b) |  (c)
Figure 5:  (a) shows basis function on four different scales of a
multiresolution basis. (b) shows the decay pattern of the Cholesky
factorization of an elliptic Green’s function discretized in this basis. (c)
shows four columns of the Cholesky factor as spatial functions. [This figure
was reproduced from Schäfer and Owhadi [3] with author permission.]

###  3.5 Fast MFM on nonsymmetric operators

As remarked by Schäfer and Owhadi [3], a LU version of Cholesky reconstruction
that extends to nonsymmetric matrices requires not only matrix–vector products
but also transpose–vector products. Similar requirements arise in hierarchical
low-rank approaches [31, 15]. Schur complementation commutes with
transposition, in the sense that

$\displaystyle\left(\overline{\mathcal{L}}\right)^{\top}=\left(\left(\mathcal{L}^{-\top}\right)_{1,1}\right)^{-1}=\left(\mathcal{L}^{\top}\right)_{1,1}-\left(\mathcal{L}^{\top}\right)_{1,2}\left(\left(\mathcal{L}^{\top}\right)_{2,2}\right)^{-1}\left(\mathcal{L}^{\top}\right)_{2,1}.$
(13)

As a result, transpose–vector products with $\overline{\mathcal{L}}$ can be
obtained by applying (inverse or forward) MFM to the
$\overline{\mathcal{L}}^{\top}$. In the case of the discretized advection-
diffusion operator in (5), when the system is solved up to time $T$, the
transpose of $\overline{\mathcal{L}}$ can be obtained by replacing $v(x,t)$
with $-v(x,T-t)$, $a(x,t)$ with $a(x,T-t)$, and by using the solution at time
$T$ as the initial condition. Here $x$ is the spatial coordinate, and $t$
denotes time. The resulting PDE is often called the adjoint problem and
frequently arises in the computation of sensitivities of solutions of PDEs
with respect to their coefficients, boundary, and initial conditions. We
empirically validate our method using matrix–vector products and
transpose–vector products obtained from a full eddy diffusivity operator
constructed via brute force (column-by-column) MFM. We leave an adjoint-based
MFM that efficiently implements transpose–vector products as future work.

##  4 Results

###  4.1 Steady-state laminar channel flow

Consider a 2D domain representing a channel with left and right walls at
$x_{1}=\pm\pi$ with Dirichlet boundary condition $c=0$ and the top and bottom
walls with $x_{2}=0,2\pi$ with no flux condition $\partial c/\partial
x_{2}=0$. The scalar field $c(x_{1},x_{2})$ is governed by a steady
advection–diffusion equation with a uniform source term

$\displaystyle u_{1}\frac{\partial c}{\partial x_{1}}+u_{2}\frac{\partial
c}{\partial x_{2}}=0.05\frac{\partial^{2}c}{\partial
x_{1}^{2}}+\frac{\partial^{2}c}{\partial x_{2}^{2}}+1,$ (14)

where the unequal diffusion constants in the coordinate directions are an
outcome of directional nondimensionalization. The flow is incompressible and
satisfies no-penetration boundary conditions on the walls. The steady velocity
field is prescribed as

$\displaystyle u_{1}=(1+\cos(2x_{1}))\cos(2x_{2}),\quad
u_{2}=\sin(2x_{1})\sin(2x_{2}).$ (15)

The PDE is discretized on a uniform staggered mesh with $N_{1}$ and
$N_{2}=N_{1}/2$ grid points in the $x_{1}$ and $x_{2}$ coordinate directions.
Second-order accurate central differences are used. The advective fluxes at
the cell faces are computed via second-order interpolation and then
multiplication of the divergence-free velocity at the face centers. At the
cell faces $x_{1}=\pm\pi$, the fluxes are computed using ghost points that
enforce the specified Dirichlet boundary conditions, while at the cell faces
at the top and bottom boundaries, the no flux condition is naturally enforced.

Figure 6:  Relative $L_{2}$ errors between the exact MFM operator
$\mathcal{D}$ and its reconstruction for the laminar flow configuration. In
(b), only the Fast MFM with automatically chosen parameterization is shown,
but for different resolutions $N_{1}$ and the number of matrix–vector products
is denoted by $n$ in the legend.

Figure 6 shows the macroscopic operator errors for the laminar flow
configuration. In (a), the mesh resolution in the $x_{1}$ (non-averaged)
coordinate is $N_{1}=2000$. Fast MFM errors are smaller than a truncated SVD
reconstruction of the same operator. The latter provides the optimal low-rank
approximation of $\mathcal{D}$ but requires access to the full operator and
is, therefore, not practical. A randomized low-rank representation is also
shown, which is available. The errors for this reconstruction are about $10$
times larger than the SVD. Compared to Fast MFM, these errors are also larger
as the number of matrix–vector products increases. The Fast MFM requires
choosing the distance between basis functions of the same color (see section
3.2) and the level at which the wavelet coefficients are truncated. The former
parameter dictates the cost-accuracy trade-off, and choosing a suitable
truncation can improve the method’s cost and stability. A sweep over a wide
range of parameters is shown in shaded markers. We use a heuristic to set
these parameters, which results in the non-shaded darker marks. Sometimes, the
heuristic still produces poor parameter choices resulting in larger Fast MFM
errors.

In fig. 6 (b), we perform a similar analysis but only show the Fast MFM
results for varying mesh sizes $N_{1}$. Errors decrease exponentially with the
number of matrix–vector products (corresponding to the number of DNSs) with
about the same fit coefficients regardless of $N_{1}$. This indicates that the
Fast MFM reconstruction is dependent on the physical locality of the operator,
not a numerical or discretized one. Thus, operator recovery for high-
resolution simulations has an out-sized benefit over traditional MFM.

Figure 7:  Reconstruction of $\mathrm{bar}{c}$ for the laminar, steady channel
flow problem. Results are shown for a randomized low-rank approximation, the
Boussinesq approximation, and the Fast MFM result. This visualization is 26
out of 2000 possible samples.

Figure 7 shows the application of the recovered operator of fig. 6 to compute
$\mathrm{bar}{c}$ using (10). The exact result is computed using $N_{1}=2000$
DNSs to recover $\mathrm{bar}{c}$. For fewer DNSs, $26$ out of the $2000$
total non-averaged degrees of freedom, the Fast MFM result matches the exact
result well, but the other methods do not. The Boussinesq approximation is
purely local [14]:

$\displaystyle-\overline{u_{1}^{\prime}c^{\prime}}=D_{\mathrm{Boussinesq}}\frac{\partial\mathrm{bar}{c}}{\partial
x_{1}},$ (16)

where

$\displaystyle
D_{\mathrm{Boussinesq}}(x_{1})=\int_{y_{1}}\mathcal{D}(x_{1},y_{1})\text{d}y_{1}.$
(17)

###  4.2 Turbulent channel flow

We next consider a fully-developed turbulent channel flow, reconstructing eddy
diffusivities, which, for momentum transport, are commonly referred to as eddy
viscosities. The incompressible Navier–Stokes equations are

$\displaystyle\frac{\partial u_{i}}{\partial t}+\frac{\partial
u_{j}u_{i}}{\partial x_{j}}=-\frac{\partial p}{\partial
x_{i}}+\frac{1}{\mbox{\text{Re}}}\frac{\partial^{2}u_{i}}{\partial
x_{j}\partial x_{j}}+r_{i},$ (18) $\displaystyle\frac{\partial u_{j}}{\partial
x_{j}}=0,$ (19)

where $\bm{r}$ is a body force, $p$ is the pressure, and $\bm{u}$ are
velocities. Following Mani and Park [1], the generalized momentum transport
equation associated with MFM for a computed $u_{i}$ field is

$\displaystyle\frac{\partial v_{i}}{\partial t}+\frac{\partial
u_{j}v_{i}}{\partial x_{j}}=-\frac{\partial q}{\partial
x_{i}}+\frac{1}{\mbox{\text{Re}}}\frac{\partial^{2}v_{i}}{\partial
x_{j}\partial x_{j}}+s_{i},$ (20) $\displaystyle\frac{\partial v_{j}}{\partial
x_{j}}=0,$ (21)

for a pressure-like term $q$ that ensures incompressibility.

Figure 8:  Diagram of the Fast MFM reconstruction procedure for the 3D
turbulent channel flow case.

We consider a case with $\mbox{\text{Re}}_{\tau}=u_{\tau}\delta/\nu=180$ where
$\delta$ is the channel half-width and $u_{\tau}$ is the friction velocity.
The mean flow is in the $x_{1}$ direction, the $x_{2}$ direction is wall-
normal, and the $x_{3}$ direction is the span with periodic boundaries. The
streamwise domain length is $2\pi$, and the spanwise length is $\pi$. The body
force, $\bm{r}$, is the mean pressure gradient in the periodic simulation and
is $\bm{r}=(1,0,0)$ in this nondimensionalized problem. The incompressible
Navier–Stokes equations are solved using direct numerical simulation with a
$144^{3}$ grid for $T=500\delta/u_{\tau}$ to ensure statistical convergence.
Simulation-result baselines, including the discretized $\mathcal{D}$ and
$\mathcal{L}$ matrices, for this case follow from Park and Mani [12] and are
used herein. Figure 8 shows our MFM procedure, averaging all independent
variables except for the wall-normal coordinate $x_{2}$. We thus recover the
Reynolds stresses as a function of $x_{2}$. The generalized eddy viscosity is
given by

$\displaystyle-\overline{u_{\beta}^{\prime}u_{\gamma}^{\prime}}(\bm{x})=\int_{\bm{y}}\mathcal{D}^{(\beta\gamma)}(\bm{x},\bm{y})\frac{\partial\mathrm{bar}{u}_{1}}{\partial
x_{2}}{\bigg{|}}_{\bm{y}}\text{d}\bm{y}$ (22)

where $\mathcal{D}^{\beta\gamma}$ in our notation represents
$\mathcal{D}_{\beta\gamma 21}$ in traditional notation for the eddy viscosity
tensor [24].

Figure 9:  Relative errors in the recovered eddy viscosity kernels (a)
$\mathcal{D}^{21}$ and (b) $\mathcal{D}^{11}$ are shown for the
$\mbox{\text{Re}}_{\tau}=180$ turbulent channel flow configuration.

Figure 9 shows the errors in the recovered eddy viscosities. The errors are
computed as the difference in operator norms between the approximate and exact
solutions to the discretized problem. The exact solution is recovered via
brute force IMFM, which computes each non-averaged degree of freedom via
forcing each $s_{i}$ independently to recover all columns of
$\mathcal{D}^{\beta\gamma}$ for each $\beta$ and $\gamma$. In fig. 9, the
trends of (a) and (b) are similar, with the Fast MFM having smaller errors
than both the SVD, which is inaccessible in practice, and a randomized low-
rank approximation of it that is accessible. The differences in errors are
small for small numbers of matrix–vector produces as the peeling procedure
removes the long-range behaviors. For larger numbers of matrix–vector
products, the difference increases. Fast MFM has a factor of about 100 smaller
errors than the low-rank approximation for 100 matrix–vector products in both
(a) and (b).

Figure 10:  Reynolds stress reconstructions (a)
$\overline{u_{1}^{\prime}u_{1}^{\prime}}$ and (b)
$\overline{u_{2}^{\prime}u_{1}^{\prime}}$ for the turbulent channel flow
configuration. Results are shown for 20 out of 144 possible samples.

Figure 10 shows the Reynolds stress reconstructions for the turbulent channel
flow. While direct computation of Reynolds stresses needs only one DNS, we use
Reynolds stresses as a metric to assess the accuracy of the recovered eddy
viscosity operator, which inevitably requires multiple simulations. The Fast
MFM results (solid, thin line) are compared with the Boussinesq approximation
and a randomized low-rank procedure. Exact results are recovered via brute-
force MFM. For 20 simulations, the difference between the exact solution and
Fast MFM is not discernible. For the same number of simulations, the low-rank
procedure does not produce a reasonable approximation for either component.
The Boussinesq approximation is a good one for the transverse stress component
$\overline{u_{2}^{\prime}u_{1}^{\prime}}$ of fig. 10 (b) but does poorly with
the $\overline{u_{1}^{\prime}u_{1}^{\prime}}$ reconstruction in fig. 10 (a).

##  5 Conclusion

This work explores a linear algebra approach to reconstructing closure
operators. Fast MFM uses sparse recovery and peeling techniques, revealing
local behaviors that crafted forcings can simultaneously recover. Results show
that tens of simulations are required to reconstruct the eddy diffusivity
operator and averaged field to visual accuracy. This contrasts against brute-
force MFM, which forces each degree of freedom and is thus prohibitively
expensive; the Boussinesq approximation, which is shown to have inaccuracies
in some test cases; randomized low-rank approximations, which, while feasible,
have poor accuracy; and even SVD, which performs worse than Fast MFM and is
inaccessible in a simulation environment. While other ongoing work focuses on
modeling the nonlocal eddy diffusivity using partial differential equations
and limited information about the exact eddy diffusivity, the Fast MFM
procedure recovers full nonlocal eddy diffusivity operators at low sample
complexity. It is a stepping stone toward the long-term goal of sample-
efficient recovery of coarse-grained time integrators.

## Acknowledgements

This work used Bridges2 at the Pittsburgh Supercomputing Center through
allocation TG-PHY210084 (PI Spencer Bryngelson) from the Advanced
Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS)
program, which is supported by National Science Foundation grants #2138259,
#2138286, #2138307, #2137603, and #2138296. SHB also acknowledges the
resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge
National Laboratory, which is supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC05-00OR22725. SHB acknowledges
support from the Office of the Naval Research under grant N00014-22-1-2519 (PM
Dr. Julie Young). FS gratefully acknowledges support from the Office of Naval
Research under grant N00014-23-1-2545 (PM Dr. Reza Malek-Madani). JL was
supported by the Boeing Company. AM was supported by the Office of Naval
Research under grant N00013-20-1-2718. The authors gratefully acknowledge
Danah Park for providing the DNS and MFM data of the turbulent channel flow
simulations and Dana Lavacot for fruitful discussions.

## References

  * Mani and Park [2021] A. Mani, D. Park,  Macroscopic Forcing Method: A tool for turbulence modeling and analysis of closures,  Physical Review Fluids 6 (2021) 054607.
  * Hamba [1995] F. Hamba,  An analysis of nonlocal scalar transport in the convective boundary layer using the green’s function,  Journal of Atmospheric Sciences 52 (1995) 1084–1095.
  * Schäfer and Owhadi [2023] F. Schäfer, H. Owhadi,  Sparse recovery of elliptic solvers from matrix–vector products,  arXiv:2110.05351 (2023).
  * Tennekes and Lumley [1972] H. Tennekes, J. L. Lumley, A first course in turbulence, MIT press, 1972.
  * Li et al. [2020] Z. Li, N. B. Kovachki, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, A. Anandkumar,  Fourier Neural Operator for parametric partial differential equations,  in: International Conference on Learning Representations, 2020, pp. 1–16.
  * Lu et al. [2021] L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis,  Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators,  Nature Machine Intelligence 3 (2021) 218–229.
  * Kraichnan [1987] R. H. Kraichnan,  Eddy viscosity and diffusivity: exact formulas and approximations,  Complex Systems 1 (1987) 805–820.
  * Liu et al. [2021] J. Liu, H. Williams, A. Mani,  A systematic approach for obtaining and modeling a nonlocal eddy diffusivity,  arXiv:2111.03914 (2021).
  * Shirian and Mani [2022] Y. Shirian, A. Mani,  Eddy diffusivity operator in homogeneous isotropic turbulence,  Physical Review Fluids 7 (2022) L052601.
  * Shende and Mani [2022a] O. B. Shende, A. Mani,  Closures for multicomponent reacting flows based on dispersion analysis,  Physical Review Fluids 7 (2022a) 093201.
  * Shende and Mani [2022b] O. B. Shende, A. Mani,  A nonlocal extension of dispersion analysis for closures in reactive flows,  arXiv:2201.10013 (2022b).
  * Park and Mani [2021] D. Park, A. Mani,  Direct calculation of the eddy viscosity operator in turbulent channel flow at $\mathrm{Re}_{\tau}=180$,  arXiv:2108.10898 (2021).
  * Altmann et al. [2021] R. Altmann, P. Henning, D. Peterseim,  Numerical homogenization beyond scale separation,  Acta Numerica 30 (2021) 1–86.
  * Hamba [2004] F. Hamba,  Nonlocal expression for scalar flux in turbulent shear flow,  Physics of Fluids 16 (2004) 1493–1508.
  * Lin et al. [2011] L. Lin, J. Lu, L. Ying,  Fast construction of hierarchical matrix representation from matrix–vector multiplication,  Journal of Computational Physics 230 (2011) 4071–4087.
  * Martinsson [2016] P.-G. Martinsson,  Compressing rank-structured matrices via randomized sampling,  SIAM Journal on Scientific Computing 38 (2016) A1959–A1986.
  * Levitt and Martinsson [2022a] J. Levitt, P.-G. Martinsson,  Randomized compression of rank-structured matrices accelerated with graph coloring,  arXiv:2205.03406 (2022a).
  * Levitt and Martinsson [2022b] J. Levitt, P.-G. Martinsson,  Linear-complexity black-box randomized compression of hierarchically block separable matrices,  arXiv:2205.02990 (2022b).
  * de Hoop et al. [2023] M. V. de Hoop, N. B. Kovachki, N. H. Nelsen, A. M. Stuart,  Convergence rates for learning linear operators from noisy data,  SIAM/ASA Journal on Uncertainty Quantification 11 (2023) 480–513.
  * Boullé and Townsend [2022] N. Boullé, A. Townsend,  Learning elliptic partial differential equations with randomized linear algebra,  Foundations of Computational Mathematics (2022) 1–31.
  * Stepaniants [2021] G. Stepaniants,  Learning partial differential equations in reproducing kernel Hilbert spaces,  arXiv:2108.11580 (2021).
  * Nelsen and Stuart [2021] N. H. Nelsen, A. M. Stuart,  The random feature model for input-output maps between Banach spaces,  SIAM Journal on Scientific Computing 43 (2021) A3212–A3243.
  * Schäfer et al. [2021] F. Schäfer, T. J. Sullivan, H. Owhadi,  Compression, inversion, and approximate PCA of dense kernel matrices at near-linear computational complexity,  SIAM Multiscale Modeling & Simulation 19 (2021) 688–730.
  * Hamba [2005] F. Hamba,  Nonlocal analysis of the Reynolds stress in turbulent shear flow,  Physics of Fluids 17 (2005) 115102.
  * Bryngelson et al. [2019] S. H. Bryngelson, K. Schmidmayer, T. Colonius,  A quantitative comparison of phase-averaged models for bubbly, cavitating flows,  International Journal of Multphase Flow 115 (2019) 137–143.
  * Bryngelson et al. [2020] S. H. Bryngelson, A. Charalampopoulos, T. P. Sapsis, T. Colonius,  A Gaussian moment method and its augmentation via LSTM recurrent neural networks for the statistics of cavitating bubble populations,  International Journal of Multiphase Flow 127 (2020) 103262.
  * Vié et al. [2016] A. Vié, H. Pouransari, R. Zamansky, A. Mani,  Particle-laden flows forced by the disperse phase: Comparison between Lagrangian and Eulerian simulations,  International Journal of Multiphase Flow 79 (2016) 144–158.
  * Ma et al. [2016] M. Ma, J. Lu, G. Tryggvason,  Using statistical learning to close two-fluid multiphase flow equations for bubbly flows in vertical channels,  International Journal of Multiphase Flow 85 (2016) 336–347.
  * Berkowicz and Prahm [1980] R. Berkowicz, L. P. Prahm,  On the spectral turbulent diffusivity theory for homogeneous turbulence,  Journal of Fluid Mechanics 100 (1980) 433–448.
  * Schäfer et al. [2021] F. Schäfer, M. Katzfuss, H. Owhadi,  Sparse Cholesky factorization by Kullback–Leibler minimization,  SIAM Journal on Scientific Computing 43 (2021) A2019–A2046.
  * Halko et al. [2011] N. Halko, P.-G. Martinsson, J. A. Tropp,  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions,  SIAM Review 53 (2011) 217–288.