maximum error bars observed, we derived an upper limit for the wiggling angle
of $\theta\leq 0.5^{\circ}$. We consider in the following that the redshifted
outflow PA is constant for all layers (velocities) at $\rm PA=295^{\circ}$,
corresponding to the PA derived here for the high-velocity emission.

Figure 21: Determination of the outflow axis PA. Top panels: Determination of
the high-velocity component radial center as a function of projected heights
above the disk for different outflow axis PA orientations. The black dotted
lines show an angle of $\pm 1^{\circ}$ with respect to the central axis.
Middle panel: 12CO PV diagram perpendicular to the outflow axis. $(V-V_{\rm
sys})$ unit is km s-1. Bottom panel: The black dashed profile shows the
emissivity integrated from the PV diagram between $(V-V_{\rm sys})=6.64$ km
s-1 to $(V-V_{\rm sys})=9.77$ km s-1, this domain is indicated by the red
dashed lines in the middle panel. The red curves show Gaussian fits used to
derive the radial positions of the two edge peaks. The vertical black dotted
lines in the bottom and middle panels indicate the positions of these peaks.
The average of these two radial positions give the radial center at this
height.

##  Appendix B Characterization of the cusps

Figure 22: Characterization of the cusps. Left panels: Individual 12CO channel
maps at different line-of-sight velocities illustrating the spatial evolution
of the cusps. Six cusps, labeled U0-U5, are clearly identified (U4 is only
visible on three channel maps). We identify the projected heights of the
bottom of the cusps (red circles) and their evolution in the channel maps
(white dashed lines). Right panel: Longitudinal PV diagram at $\delta
x=0^{\prime\prime}$ averaged over a slice of $\Delta Z=0.2^{\prime\prime}$.
The velocity and height locations of the cusps derived in the channel maps are
shown as red dots. The four orange crosses indicate the maximum on-axis
projected velocity (back side) of the top of the fit four bright elliptical
structures identified in transverse PV diagrams at the corresponding heights
(see Fig. 23).

##  Appendix C Ellipse fits in the transverse PV diagrams

In this section, we discuss an alternative method to derive the radial
component of the velocity, not constrained by our tomographic method. At a few
positions along the flow, transverse PV diagrams clearly show elliptical
structures nested inside the main emission (see Fig. 23). We assumed that each
of these ellipses traces one emitting layer of the outflow. If the outflow
velocity does not vary drastically with height, one layer of the outflow of
fixed radius will be projected as an inclined elliptical structure in the PV
diagram (Louvet et al. 2018). The velocity width of the ellipse at $x=0$ is
directly linked to the radial velocity component of the layer. We fit ellipses
by the naked eye to the structures observed. From these fits, we recovered the
radial velocity $V_{\rm R}$ as well as $V_{\rm Z}$ at a few specific positions
along the outflow.

Figure 23 shows the $V_{\rm Z}$ profile of the outflow reconstructed with the
tomographic technique and extrapolated until Z $\approx$ 2700 au. The red
arrows trace the poloidal vectors derived from individual elliptical fits. The
angle of the arrow is determined by the ratio between $V_{\rm Z}$ and $V_{\rm
R}$. The length and colors of the arrows correspond to $V_{\rm Z}$. The
$V_{\rm Z}$ values derived from the elliptical fits appear consistent with the
estimates from the tomographic study. The poloidal velocity direction is also
consistent with the conical lines of constant $V_{\rm Z}$ in the tomography.
This suggests that the flow is indeed aligned with these conical lines.

Figure 23: Constraints on the radial velocity component $V_{\rm R}$. Left
panels: Transverse PV diagrams averaged over a slice of width $\Delta
Z=0.2^{\prime\prime}$ at selected positions $\delta z$ along the flow. In red
are shown the ellipse fits. We show side by side the PV diagram with and
without the fit for more visibility. Right panel: $V_{\rm Z}$ tomography of
the outflow in the outflow referential, extrapolated beyond $Z>1200$ au. The
red arrows represent the velocity directions determined from the ellipse fits.
The arrow colors trace the $V_{\rm Z}$ values determined from the ellipse fit.
$R$ and $Z$ units are in au.

We studied the relation between the ellipses visible in the transverse PV
diagrams and the cusps visible in the channel maps. The on-axis maximal
velocity of the fit ellipses was compared to the cusp location in the on-axis
longitudinal PV diagram (see Fig. 22). The (R,Z) positions of the fit ellipses
were also compared to the cusp locations in the tomographic $V_{Z}$ map (see
Fig. 13). The ellipse located at $\delta z=1.1^{\prime\prime}$ comes from a
lower altitude region, not included in our cusp analysis because of crowding.
The ellipses at $\delta z=2.1^{\prime\prime}$ and $2.9^{\prime\prime}$ seem to
extend the U5 cusp, at higher and lower line-of-sight velocities respectively.
The ellipse at $\delta z=7.5^{\prime\prime}$ is not so clearly associated with
a cusp extension in Fig. 22. It could be the high velocity extension of the U3
cusp, or associated with a fainter cusp located between U2 and U3, not
included in our analysis. Moreover, Fig. 24 indicates that the tops of the
ellipses located at $\delta z=2.9^{\prime\prime}$ and $\delta
z=7.5^{\prime\prime}$ are possibly consistent with a cusp structure which was
not included in our study, due to their low S/N. Similarly, Fig. 24 shows that
the cusps U1 and U2 also coincide with internal structures in the transverse
PV diagram, but the stacking of ellipses in the PV diagram makes the
identification more confusing than in the channel maps.

The cusps could be precisely located only at moderate projected velocities
$<5.3$ km s-1, where they are sufficiently bright (see Fig. B.1), hence they
probe outer streamlines; in contrast elliptical structures in transverse PV
diagrams are best distinguished at high projected velocities $>5.3$ km
s-1where they have less overlap with each other (see Fig. C.1), hence most of
them probe inner faster flow streamlines. However, we see clear
correspondences between faint cusps in channel maps and the on-axis high-
velocity portion of some ellipses in transverse PV cuts, and vice-versa, which
demonstrates that they trace different portions of the same underlying
substructures, extending across the whole conical outflow. Therefore, assuming
a radial flow to deproject the cusp apparent positions seems fully justified.

Figure 24: Illustration of the correspondences between ellipses and cusps.
Left panels: transverse PV diagrams at $\delta z=2.9^{\prime\prime}$ and
$7.5^{\prime\prime}$ (from Fig 23), with red dots indicating the on-axis
projected velocities of cusps U1 and U2 at $\delta z=7.5^{\prime\prime}$. Each
of these two cusps corresponds in position and velocity to the maximum
velocity (back side) of an ellipse in the PV cut. Right panels: channel maps
at the same line-of-sight velocity as the orange crosses in the left panels,
showing that elliptical structures seem to correspond to fainter cusps in
position-position space, but falling outside of the velocity domain where cusp
characterization was possible.

##  Appendix D Biases in the tomographic reconstruction

We studied possible biases in the method presented in Sect. 3 to reconstruct
poloidal maps of $V_{\rm Z}$ and $j=R\times V_{\phi}$. We first studied the
bias introduced by projection effects considering one single layer of the
outflow using the analytical solutions from Eqs. 1,2,3. To study the impact of
beam convolution and multiple layers of the outflow, we also applied our
tomographic reconstruction method to the synthetic data cube presented in
Sect. 10 and estimate the difference between the reconstructed $V_{\rm Z}$ and
$j$ and the initial values of the model.

###  D.1 Single shell

Figure 25: Biases in the tomographic reconstruction for one single conical
shell. Top panels: Black curves show computed transverse PV diagrams for one
conical layer with constant $V_{\rm Z}$, $V_{\rm R}$, and $j$ seen at
different inclinations $i$ and with $\theta_{V}=\theta=17^{\circ}$. The red
dashed curve traces an elliptical fit. Bottom Panel: Computed relative bias in
$V_{\rm Z}$ (green curves) and $j$ (orange curves) for one conical layer as a
function of the inclination and for different combinations of $\theta$ and
$\theta_{v}$ : solid curves: $\theta=\theta_{v}=17^{\circ}$, dotted curves:
$\theta=17^{\circ}$ & $\theta_{v}=30^{\circ}$, and dashed curves:
$\theta=\theta_{v}=30^{\circ}$.

In our tomographic study, we assumed that one shell will be projected as an
ellipse in transverse PV diagrams. We also assumed that the projected
velocities at the extrema radii allowed us to recover $V_{\rm Z}$ and
$V_{\phi}$. We determine in this section biases introduced by these two
assumptions.

We assumed a conical shell of radius R(Z) in the outflow referential and with
local opening angle $\theta$, such as $R(Z)=R_{0}+Z\tan{\theta}$. For the DG
Tau B inner conical outflow, the opening angle $\theta$ of the layers vary
between 17∘ and 12∘, increasing with decreasing velocities. The estimated
opening angle of the lower velocity emission contributing to the pedestal is
$\theta\simeq 30^{\circ}$. We defined the parameter $\theta_{v}$ corresponding
to the angle of the velocity vector with the Z axis in the poloidal plane
($V_{\rm R}=V_{\rm Z}\tan{\theta_{V}}$). For a flow parallel to the conical
surface, $\theta=\theta_{V}$. We also assumed a constant specific angular
momentum ($V_{\phi}=j/R(\phi)$) as well as constant $V_{\rm Z}$ and $V_{\rm
R}$ over the transverse slit width. Solving Eq. 2 with constant $\delta z$
gave:

$R(\phi)=\frac{\delta
z\tan{\theta}+R_{0}\sin{i}}{\sin{i}-\sin{\phi}\cos{i}\tan{\theta}}.$ (21)

In the case $i=\theta$ and $i<\theta$, the cut of the conical outflow will be
respectively a parabola and an hyperbola. Consequently, $R(\phi)$ will tend
toward infinite. This is not the case for the DG Tau B conical outflow. We
then implemented this solution into Eqs. 3,1 for multiple values of the
inclination. The resulting transverse PV diagrams are represented in Fig. 25
for $\theta_{V}=\theta=17^{\circ}$. The PV diagrams were generated using Eqs.
3, 2 with the solution from Eq. 21. The difference with an ellipse increases
with the inclination but is expected to be small for the DG Tau B case
($i=117^{\circ}$).

We assume in Sect. 3 that the radial edges of the ellipse correspond to
$\phi=0-\pi$. However, this is an approximation. The radial extrema of the
ellipse correspond to solutions of the equation: $\frac{\partial\delta
x}{\partial\phi}=0$. For a conical layer, the solution is obtained for:

$\sin{\phi}=\frac{\tan{\theta}}{\tan{i}}.$ (22)

This difference is small in our situation. However, this leads to bias in the
estimation of $R$, $Z$, $V_{\rm Z}$ and $V_{\phi}$. We determined the
estimated value of $V_{\rm est,Z}$ and $V_{\rm est,\phi}$ using Eqs. 6,7. We
then computed the relative differences with the real values $R_{\rm real}$,
$Z_{\rm real}$, $V_{\rm real,Z}$, and $V_{\rm real,\phi}$:

$\frac{R_{\rm est}-R_{\rm real}}{R_{\rm
real}}=\sqrt{1-\frac{\tan^{2}{\theta}}{\tan^{2}{i}}}-1$ (23) $\frac{Z_{\rm
est}-Z_{\rm real}}{Z_{\rm real}}=-\frac{\tan^{2}\theta}{\tan^{2}i}$ (24)
$\frac{V_{\rm est,Z}-V_{\rm real,Z}}{V_{\rm
real,Z}}=\tan{\theta}\times\tan{\theta_{V}}$ (25) $\frac{V_{\rm
est,\phi}-V_{\rm real,\phi}}{V_{\rm
real,\phi}}=\sqrt{1-\frac{\tan^{2}{\theta}}{\tan^{2}{i}}}-1.$ (26)

Hence, $R$, $Z$, and $j$ are systematically underestimated, while $V_{\rm Z}$
is overestimated. Figure 25 shows the relative biases due to inclination and
projection effects in the estimation of $V_{\rm Z}$ and $V_{\phi}$ for
different values of $\theta$ and $\theta_{v}$. In the conical outflow, where
$\theta_{V}=\theta\leq 17^{\circ}$, the inclination bias is expected to be
$\leq 10\%$ in $V_{\rm Z}$ and $Z$, $\leq$ 1 % in $j$ and $R$.

###  D.2 Multiple layers

Figure 26:  Biases in the tomographic reconstruction for a stacking of conical
shells. Left panels: The gray area shows the limits of the transverse PV
diagram predicted for a stacking of conical shells with a shearing factor
$f_{\rm sh}$ (see text). In black is shown the ellipse corresponding to a
single layer of maximal radius $R$. The blue and red dots correspond
respectively to the velocity of the ellipse at $\pm R$ and the maximal
velocity of the PV diagram at $\pm R$. The red dots are the ones used in our
method, while the blue dots correspond to the true measurements. The
differences in velocities will lead to a bias in our $j$ and $V_{\rm Z}$
estimation. The red line shows the slope of the PV diagram $\frac{\partial
V_{\rm los}(\pm R)}{\partial R}$, see text for more details. Right panel: Plot
of the relative bias in $V_{\rm Z}$ (green curves) and $j$ (orange curves) due
to the shearing aspect of the PV diagram as a function of the $x$ parameter
(see text). The solid, dashed, and dotted lines correspond to different
configurations of $\theta$ and $\theta_{v}$. The black histogram, labeled
C.O., shows the distribution of the $x$ parameter derived at different heights
Z and radii R in the conical outflow. The relative biases corresponding to the
a and b ellipses shown in the left panels are represented. The orange (resp.
green) colored contours trace the relative biases derived from applying the
tomographic method to the synthetic data cube. Contours outline 30 to 90 % of
the distributions in step of 15 %.

In the previous section, we estimated the bias in the estimation of $V_{\rm
Z}$, $V_{\rm\phi}$ due to projection effects for one single conical layer.
However, the DG Tau B transverse PV diagrams shows a clear shear-like velocity
structure suggesting a stacking of layers. We modeled this effect directly in
the transverse PV diagrams by stacking the elliptical projections for conical
layers of increasing radii at origin $R_{0}$ with the same opening angle
$\theta$. The velocity shear in $V_{\rm Z}$ of the conical layers was defined
by a shearing parameter $f_{\rm sh}=\frac{\partial V_{\rm Z}}{\partial R}$.
The specific angular momentum of each layer, $j$, was assumed constant with
$Z$ and vary between 20 and 90 au km s-1, increasing with increasing radius,
such as $j\times V_{\rm P}$ was kept constant to mimic the DG Tau B
observations.

Due to the stacking, the maximal velocity at $\pm R$ does not perfectly
describe the velocity of the ellipse of radius R. This effect is larger when
$f_{\rm sh}$ is small. Using a similar procedure as before, we derived the
relative difference between the estimation of $V_{Z}$ and $j$ with the
tomographic method and the input theoretical values. This was achieved for a
range of $f_{\rm sh}$ values between 0.2 to 20 km s-1 au-1 and with the three
($\theta$, $\theta_{V}$) configurations studied in the previous subsection.

In order to study efficiently this bias, we defined the a-dimensional
parameter $x=\frac{[V_{\rm los}(\pm R)]}{R}$ $\times$ ($\frac{\partial[V_{\rm
los}(\pm R)]}{\partial R}$ )-1, where we defined $[V_{\rm los}(\pm
R)]=1/2(V_{\rm los}(R)+V_{\rm los}(-R))$. This a-dimensional parameter can be
derived from the observations. We show in Fig. 26 the predicted relative
biases in $V_{\rm Z}$ and $j$ as a function of this parameter $x$. Biases
increase with increasing $x$ values illustrating the effect of the velocity
shear. We also show the distribution of observed x values computed at each
(R,Z) position in the conical outflow. The x values are concentrated around
$\simeq 0.5$, suggesting moderate biases in both $V_{\rm Z}$ and $j$ for the
conical outflow where $\theta_{V}=\theta\leq 17^{\circ}$ .

However, this modeling did not include the effect of beam smearing nor the
impact of our polynomial fitting method to describe the shape of the PV
diagram. In order to study these effects, we directly applied our tomographic
method to the synthetic data cube presented in Sect. 10 and determined the
relative differences between the computed and the input $V_{\rm Z}$ and $j$
values at each (R,Z) position along the conical flow. We show in Fig. 26 the
derived relative biases in $V_{\rm Z}$ and $j$ as a function of the $x$
parameter. Their distributions are broader and flatter than predicted,
especially in $j$, likely due to the combination of velocity shear and beam
smearing effects.

From this study, we estimate that in the conical part of the outflow (at
$(V-V_{\rm sys})\geq 2$ km s-1), the tomographic method suffers from a
relative bias $\leq 15$% in the estimation of $V_{\rm Z}$ and $\leq 20$% in
the estimation of $j$. The X values extend up to $\approx 1.3$ in the pedestal
region. With an opening angle $\theta$, $\theta_{v}>30^{\circ}$, a tomographic
study of the low velocity component would suffer from a relative bias larger
than 60% and 30% in the estimation of $V_{\rm Z}$ and $j$ respectively.

##  Appendix E Wind-driven shell analytical solutions

In this section, using Eqs. 1, 2,3 and 11, we derive an analytical solution
for predicted channel maps in the case of the generalized WDS model introduced
in Sect. 4. The WDS model is defined by three parameters: $C$, $\tau$ and
$\eta$, see Eqs. 11. The projection on the plane of the sky ($\delta x$,
$\delta z$) for the emissivity map at $(V-V_{\rm sys})=V_{\rm CM}$ could be
recovered by solving Eq. 3 with $V_{\rm los}=V_{\rm CM}$ and using Eqs. 11 for
$V_{\rm Z}$ and $V_{\rm R}$:

$\displaystyle\zeta$ $\displaystyle=$ $\displaystyle\frac{4V_{\rm CM}\tau
C\cos{i}}{\eta^{2}\sin^{2}{i}}$ (27) $\displaystyle R$ $\displaystyle=$
$\displaystyle-\eta\frac{\tan{i}}{2C}(\sin{\phi}\pm\sqrt{\sin^{2}{\phi}-\zeta})$
(28) $\displaystyle\delta x$ $\displaystyle=$ $\displaystyle R\cos{\phi}$ (29)
$\displaystyle\delta z$ $\displaystyle=$ $\displaystyle
CR^{2}\sin{i}-R\sin{\phi}\cos{i}.$ (30)

Equation 28 has no solution in the case $\zeta>1$ ($V_{\rm
CM}\cos{i}>\frac{\eta^{2}\sin^{2}{i}}{4\tau C}$). This corresponds to the case
where no emission is predicted at $V_{\rm los}=V_{\rm CM}$. Similarly, in the
situation $0\geq\zeta\geq 1$, only a fraction of $\phi$ will be projected such
as $\sin^{2}{\phi}>\zeta$. Developing Eq. 30, we obtained the following:

$\displaystyle\delta z$ $\displaystyle=$
$\displaystyle\eta\frac{\sin^{2}{\phi}\sin{i}}{2C}(\eta\tan^{2}{i}+1)-V_{CM}\tau\tan{i}$
(32)
$\displaystyle\pm\eta\frac{\sin{\phi}\sin{i}}{2C}(\eta\tan^{2}{i}+1)\sqrt{\sin^{2}{\phi}-\zeta}.$

We reformulated $\delta z$ as $z_{0}+\Delta Z\sin{\beta}$ with:

$\displaystyle z_{0}$ $\displaystyle=$
$\displaystyle\eta\frac{\sin{i}}{2C}(\eta\tan^{2}{i}+1)-V_{CM}\tau\tan{i}$
(33) $\displaystyle\Delta Z$ $\displaystyle=$
$\displaystyle\eta\frac{\sin{i}}{2C}(\eta\tan^{2}{i}+1)\sqrt{1-\zeta}$ (34)
$\displaystyle\sin{\beta}$ $\displaystyle=$
$\displaystyle\frac{-1+\sin^{2}{\phi}\pm\sin{\phi}\sqrt{\sin^{2}{\phi}-\zeta}}{\sqrt{1-\zeta}}.$
(35)

Similarly, manipulating Eq. 29, $\delta x$ could be reformulated as $\delta
x=\Delta X\cos{\beta}$, with the following:

$\Delta X=\mp\eta\frac{\tan{i}}{2C}\sqrt{1-\zeta}.$ (36)

Hence ($\delta x$, $\delta z$) trace an ellipse of center (0, $z_{0}$) and
aspect ratio:

$\Bigl{|}\frac{\Delta Z}{\Delta X}\Bigr{|}=(\eta\tan^{2}{i}+1)|\cos{i}|.$ (37)

Therefore, in the classical WDS models with radial velocity vectors
($\eta=1$), the aspect ratio of the ellipse on the channel maps only depends
on the inclination. We represent on the channel maps in Figs. 5, 7, and 8 the
two limiting ellipses computed at $V-\delta V/2$ and $V+\delta V/2$ to take
into account the width of the channel map.

##  Appendix F Global model

Figure 27: Schematic representation of the method used to create synthetic
data cube.

In order to comprehend the impact of projection or convolution effects on the
observations, we developed a locally axisymmetric code that allows us to
simulate optically thin observations of simple outflow models. We took
advantage of the axisymmetric hypothesis that permits to reduce the complexity
of the model by defining at each height the radius $R(Z)$ and the cylindrical
velocities $V_{\rm z}$, $V_{\rm r}$ and $V_{\phi}$ (See the schematic view of
Fig. 2). The dependency between the height and the radius or the velocities
vary with the model used. We then created at each height an emitting ring with
azimuth parametrized with $\phi$. As the morphology is axisymmetric, neither
the radius nor velocities depend on $\phi$. Under the assumption of optically
thin emission, the emissivity at position ($Z$,$\phi$) is proportional to the
elementary volume $dV=R(Z)d\phi dRdZ$. We added an additional variation of the
emissivity with the height and radius as a power-law with parameters $\alpha$
and $\beta$ respectively. Proper modeling of the emissivity would require the
temperature and chemistry to be solved, which is well beyond the scope of this
model. The positions of the outflow emission on the data cube ($\delta
x$,$\delta y$, $V_{\rm proj}$) were then defined by Eqs. 1,2,3.

We then created a data cube with the same spectral and spatial resolution than
our observations, and placed on each point ($\delta x$,$\delta y$, $V_{\rm
proj}$) the emissivity I. Under the assumption of optically thin emission, we
summed each emissivity corresponding to the same positions on the data cube.
We set a step size of 1∘ for $\phi$ and a fraction of the spatial pixel for Z.
We then convolved the data cube by a 2D Gaussian matching the spatial beam
characteristics in order to fully simulate the ALMA observations. The code,
written in Python 3, is publicly online777https://github.com/Alois-
deValon/Axoproj.

##  Appendix G Channel maps of disk-wind models

Figure 28: 12CO observed (first column) and synthetic disk-wind channel maps
at selected line-of-sight velocities. The second and third columns correspond
to the two disk-wind models with axis precession presented in Sect. 5.2. The
last column corresponds to a disk-wind model with axisymmetric density
enhancements, as presented in Sect. 5.3. $\delta x$ and $\delta z$ units are
arcseconds.


2023 1

ATIF et al. Towards Accelerating Particle-Resolved Direct Numerical Simulation
with Neural Operators

Mohammad Atif,

# Towards Accelerating Particle-Resolved Direct Numerical Simulation with
Neural Operators

Mohammad Atif  Vanessa López-Marrero  Tao Zhang  Abdullah Al Muti Sharfuddin
Kwangmin Yu  Jiaqi Yang  Fan Yang  Foluso Ladeinde  Yangang Liu  Meifeng Lin
Lingda Li  Computational Science Initiative, Brookhaven National Laboratory,
Upton, New York 11973, USA  IACS, Stony Brook University , New York, USA
Environmental and Climate Sciences Department, Brookhaven National Laboratory,
Upton, New York 11973, USA  Department of Mechanical Engineering, Stony Brook
University , New York, USA  Department of Mathematics, Emory University,
Atlanta, GA 30322, USA<EMAIL_ADDRESS>
(Date Month Year; Date Month Year; Date Month Year)

###### Abstract

[Abstract] We present our ongoing work aimed at accelerating a particle-
resolved direct numerical simulation model designed to study aerosol-cloud-
turbulence interactions. The dynamical model consists of two main components –
a set of fluid dynamics equations for air velocity, temperature, and humidity,
coupled with a set of equations for particle (i.e., cloud droplet) tracing.
Rather than attempting to replace the original numerical solution method in
its entirety with a machine learning (ML) method, we consider developing a
hybrid approach. We exploit the potential of neural operator learning to yield
fast and accurate surrogate models and, in this study, develop such surrogates
for the velocity and vorticity fields. We discuss results from numerical
experiments designed to assess the performance of ML architectures under
consideration as well as their suitability for capturing the behavior of
relevant dynamical systems.

††articletype: Research Article††journal: Journal††volume: 00

##  1 Introduction

Particle-resolved direct numerical simulations (PR-DNS), which resolve not
only the smallest turbulent eddies but also track the development and motion
of individual particles, are essential tools for studying aerosol-cloud-
turbulence interactions 1, 2. The associated mathematical models comprise
systems of partial differential equations (PDEs), including the Navier-Stokes
equations for fluid flow. In turbulent regimes, the numerical solution of the
Navier-Stokes equations comes at a high computational cost. To reduce the cost
of such high-fidelity direct numerical simulations (DNS), surrogate models (of
lower fidelity) are often employed. With judicious algorithmic choices, high-
fidelity models may be used in conjunction with surrogate models to reduce
computing costs, while retaining acceptable levels of accuracy in the
simulation results 3, 4.

The confluence of machine learning and computational fluid dynamics has led to
important progress in increasing the speed of DNS 5, turbulence modeling 6,
developing surrogate models 7, and learning patterns and PDEs from data 8, 9,
10. Among many recent advancements in machine learning (ML) and artificial
intelligence, deep neural networks for operator learning, which learn mappings
between function spaces 11, 12, 13, 14, have emerged as new scientific
computing tools with much potential to aid in the development of fast and
accurate surrogate models for dynamical systems 12, 15. Research with such
neural operator networks is currently very active in many application areas,
including climate sciences16, 17.

In the present study, we consider the use of such neural operator methods in
order to assess their performance as surrogates of a numerical PR-DNS model.
Due to the complexity of the PR-DNS model of interest to us, as a first step
we focus on a subset of equations used in the PR-DNS model rather than
attempting to replace the physical model in its entirety with a machine
learning surrogate. This will allow us to study the potential benefits of a
hybrid approach, where the neural operator surrogates are judiciously
incorporated within the numerical PR-DNS model.

We proceed to describe the PR-DNS model in Section 2. Our initial test case
and data set are presented in Section 3. The methodology and choice of neural
operator method appears in Section 4, followed by a discussion of the obtained
results in Section 5. Concluding remarks and future research directions are
given in Section 6.

##  2 Particle-Resolved DNS Model

A summary of the physical model ultimately of interest to us appears next.
Readers may refer to 1 for a thorough description of the model.

###  2.1 Fluid Dynamics Equations

The incompressible Boussinesq system

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\nabla)\mathbf{u}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho_{0}}\nabla
p+\nu\nabla^{2}\mathbf{u}+\mathbf{f}_{b}+\mathbf{f}_{e}$ (1)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=$ $\displaystyle 0$ (2)

is adopted to model the fluid velocity. In equations (1)–(2), $\mathbf{u}$ is
the fluid velocity field, $p$ denotes pressure, $\nu$ is the kinematic
viscosity, and $\rho_{0}$ is the density of dry air. The buoyancy force
$\mathbf{f}_{b}$ in (1) is given by

$\displaystyle\mathbf{f}_{b}$ $\displaystyle=$
$\displaystyle-\mathbf{g}\left[\frac{T-T_{0}}{T_{0}}+0.608(q_{v}-q_{v_{0}})-q_{c}\right],$
(3)

where $\mathbf{g}$ denotes the gravitational acceleration vector, $T$ is
temperature, $q_{v}$ is the water vapor mixing ratio, and $q_{c}$ is the
liquid water mixing ratio. $T_{0}$ and $q_{v_{0}}$ denote reference values. As
the buoyancy force $\mathbf{f}_{b}$ in (3) depends on the temperature $T$ and
the water vapor mixing ratio $q_{v}$, it couples the equations (1)–(2) for
fluid velocity with those below for temperature (4) and water vapor (5).
Finally, the external force $\mathbf{f}_{e}$ in (1) is introduced to maintain
a statistically stationary homogeneous turbulence 1.

The time evolution for the temperature $T$ and water vapor mixing ratio
$q_{v}$ is modeled by the equations

$\displaystyle\frac{\partial T}{\partial t}+(\mathbf{u}\cdot\nabla)T$
$\displaystyle=$ $\displaystyle\frac{L_{h}}{c_{{p}}}C_{d}+\mu_{T}\nabla^{2}T$
(4) $\displaystyle\frac{\partial q_{v}}{\partial
t}+(\mathbf{u}\cdot\nabla)q_{v}$ $\displaystyle=$ $\displaystyle-
C_{d}+\mu_{v}\nabla^{2}q_{v},$ (5)

where $L_{h}$ is the latent heat of water vapor condensation, $c_{p}$ is the
specific heat at constant pressure. $\mu_{T}$ and $\mu_{v}$ are the molecular
diffusivities for temperature and water vapor, respectively. The rate of
exchange between liquid and vapor, a.k.a. the condensation rate $C_{d}$, as
defined in equation (9), depends on droplet radii and thus couples the fluid
dynamics model with the equations for droplet growth and motion depicted in
Section 2.2 below.

###  2.2 Droplet Growth and Motion

The equations describing the condensation/evaporation of cloud droplets and
their motion are given by

$\displaystyle R_{i}(t)\,\frac{\mathrm{d}R_{i}(t)}{\mathrm{d}t}$
$\displaystyle=$ $\displaystyle A\cdot S(\mathbf{X}_{{i}},t)$ (6)
$\displaystyle\frac{\mathrm{d}\mathbf{X}_{{i}}(t)}{\mathrm{d}t}$
$\displaystyle=$ $\displaystyle\mathbf{V}_{{i}}(t)$ (7)
$\displaystyle\frac{\mathrm{d}\mathbf{V}_{{i}}(t)}{\mathrm{d}t}$
$\displaystyle=$
$\displaystyle\frac{1}{\tau_{p}}[\mathbf{u}(\mathbf{X}_{{i}},t)-\mathbf{V}_{{i}}(t)]+\mathbf{g}.$
(8)

Here, $R_{i}(t)$, $\mathbf{X}_{{i}}(t)$, and $\mathbf{V}_{{i}}(t)$ denote,
respectively, the radius, position coordinate, and velocity of the $i$-th
droplet; $A$ is a function of temperature and pressure; $S(\mathbf{X},t)$ is
the supersaturation; $\mathbf{g}$ is the gravitational acceleration vector;
and $\tau_{p}$ is a measure of the droplet inertial effect.

The condensation rate $C_{d}$ in equations (4)–(5) is defined as

$\displaystyle C_{d}(\mathbf{X},t)$ $\displaystyle=$
$\displaystyle\frac{4\pi\rho_{l}A}{\rho_{0}a^{3}}\,\sum_{i=1}^{n}S(\mathbf{X}_{{i}},t)R_{i}(t)\,,$
(9)

where $a$ is the size of a grid cell, $n$ is the number of droplets in the
grid cell, $\rho_{l}$ and $\rho_{0}$ are the densities of water and air, and
$R_{i}(t)$ is the radius of the $i$-th droplet. Note that the condensation
rate $C_{d}$ acts as a source (or sink) term for equations (4)–(5) and thus
couples the particle equations with the fluid dynamics model. Finally, the
full model (1)–(8) is augmented with suitable initial and boundary conditions.
Typical simulations of the PR-DNS model from Gao et. al (2018) track evolution
of approximately 10 million cloud particles in a physical domain of
$0.512m\times 0.512m\times 0.512m$ meshed onto a grid of $256\times 256\times
256$.

##  3 Data Set

As discussed in Section 2, the PR-DNS model captures the basic fields of fluid
dynamics like velocity, temperature, and pressure as well as the cloud physics
fields of vapor mixing ratio, supersaturation, and condensation rate. It also
tracks the location and size of cloud particles. As the first step to model
such a complex system, this paper aims to ascertain whether the Navier-Stokes
solver governing the velocity field can be replaced by a surrogate model. To
this effect, we select a simpler test case of doubly periodic shear layer
where the entire dynamics can be captured by a vorticity or a velocity field.

The doubly periodic shear layer is a well-known test case of two dimensional
incompressible flow simulations. This test has been extensively employed to
assess the errors and stability of fluid dynamics models 18, 19. Here, an
incompressible fluid of kinematic viscosity $\nu$ is confined in a two
dimensional square domain of side $L=1$. Periodic boundary conditions are
applied in both $x$ and $y$ spatial directions. The velocity fields at the
initial time $t=0$ are given by

$\displaystyle u_{x}(y)$
$\displaystyle=\begin{cases}U_{0}\tanh[(4y-1)/w],\qquad y\leq 1/2\\\
U_{0}\tanh[(3-4y)/w],\qquad y>1/2\end{cases}$ (10) $\displaystyle u_{y}(x)$
$\displaystyle=U_{0}\delta\sin[2\pi(x+1/4)],$ (11)

where $u_{x},u_{y}$ are the velocities in $x,y$ directions respectively. The
parameters $w$ and $\delta$ control the width of the shear layer and the
magnitude of the perturbation. The Reynolds number ($\rm Re$) defined as ${\rm
Re}=U_{0}L\,/\,\nu$ represents the ratio between inertial and viscous forces.
The perturbation causes the shear layers to curl up into two vortices. We
modify the initial condition using random numbers to

$\displaystyle u_{x}=\begin{cases}U_{0}\,r_{1},\qquad y\leq 1/2,\\\
U_{0}\,r_{2},\qquad y>1/2,\end{cases}$ (12) $\displaystyle
u_{y}=\begin{cases}U_{0}\,r_{3},\qquad x\leq 1/2,\\\ U_{0}\,r_{4},\qquad
x>1/2,\end{cases}$ (13)

where $r_{1},r_{2},r_{3},r_{4}$ are uniformly distributed random numbers
between $[-1,1]$. This generates two opposite vortices that diffuse as time
progresses and convect around the periodic domain depending on their initial
locations.

The data set for the present study corresponds to Reynolds number $1000$. We
solve the Navier-Stokes equations for the velocity field using the entropic
lattice Boltzmann method 20, 21 on a grid of $64\times 64$ points. The flow
field is first allowed to evolve for $4\,t_{c}$ (where $t_{c}=L/U_{0}$ is the
convection time) so that the initial sharp discontinuities vanish. Thereafter,
time is reset to zero and the velocity (${\bf u}$) and vorticity
($\omega_{z}$) are sampled from time $t=0$ to $t=t_{c}$ in steps of
$0.02\,t_{c}$. Note that the vorticity field $\omega_{z}(x,y)$ is calculated
as the curl of the velocity using

$\displaystyle\omega_{z}(x,y)=\nabla\times{\bf u}(x,y),$ (14)

where $\nabla$ is the nabla operator and the velocity vector ${\bf u}$ has
components $u_{x},u_{y}$ for two-dimensional fluid flow. The dataset consists
of fields from 10,000 such simulations. The initial condition of the velocity
field [Eqns. (12),(13)] depend on uniformly distributed random numbers that
vary from one sample to another. Figure 1 visualizes the two dimensional
vorticity field of two such samples. The two oppositely rotating vortices are
visible in the figure. In general, in an incompressible viscous flow the
vortices stretch (absent for two dimensions) and turn under the influence of
each other’s field and diffuse as the time advances.

Figure 1: Visualization of vorticity field from two samples of a dataset of
10,000 simulations at ${\rm Re}=1000$. The red and blue vortices are rotating
clockwise and counter-clockwise respectively. The slices at different non-
dimensional times $t/t_{c}$ show the motion of the voritices for unit
convection time.

##  4 Methodology

In recent years, many studies have been conducted to build and evaluate deep
learning models as alternatives to traditional PDE solvers for a wide variety
of problems in transport phenomenon. They have been found particularly useful
in replacing components of multiscale problems with identifiable computational
bottlenecks or in modeling flows for which experimental data is available but
constitutive relationships are unknown. One of the most popular deep learning
approach for inverse problems is the physics-informed neural network (PINN)
22, which incorporates the governing PDEs in the loss function of the
minimization problem during training. PINNs are successful in solving inverse
problems for the transport coefficients they are trained on, however, they too
sometimes suffer from high computational expense and lack of generalizability.
These lacunae of PINNs are attributed to the fact that they learn mappings
between finite dimensional Euclidean spaces.

Neural operators 11, 12, 13, 14, which have gained popularity recently,
alleviate the limitations of PINNs. Particularly, the Fourier neural operators
(FNOs) 12, 14 learn mappings between infinite-dimensional function spaces. A
brief description of FNOs is provided in Appendix A. FNOs have become a
popular alternative to the conventional neural networks for a wide range of
physical applications like climate modeling 17, multiphase flows in porous
media 23, 24, and wave propagation 25. This popularity is due to their low
computational cost, relatively small errors, and the support of features like
zero-shot super-resolution for turbulent flows which are limited in other
machine learning methods. Li et al. 12 proposed two types of FNOs for
transient dynamics in two spatial dimensions:

  * •

The two dimensional FNO coupled with a recurrent neural network architecture
for recurrent propagation in time, which maps the input of a few consecutive
time steps to the next chronological time step.

  * •

The three dimensional FNO, which learns a mapping in two spatial and one
temporal dimensions. This maps an input time window to an entire output time
window in a single inference step.

In the present study, we choose the three dimensional FNO as we want to
predict fields over a long time interval without accumulating errors from step
by step temporal predictions. The Fourier series is truncated at a finite
number of lower modes while filtering the higher modes. The FNO model used
here has four Fourier layers where width represents the number of nodes in a
layer. The code is available on Github 12 from the original authors of FNO and
has also been implemented in NVIDIA’s Modulus package 26. We design a simple
study to understand the behavior and errors of surrogate models. As discussed
previously, the data set contains velocity and vorticity fields from 10,000
simulations of doubly periodic shear layer corresponding to Reynolds number
1000. We calculate the errors for different number of samples, input-output
training intervals, and time refinements. The ${\rm L}_{2}$ errors reported in
the following section are calculated at each time step in the prediction
interval and averaged over $N=500$ test samples. We note that the same test
set of size 500 is used for all numerical experiments that follow. The fields
are first normalized via the target data’s mean and standard deviation and
then the error is calculated as

$\displaystyle{\rm L}_{2}\,\,{\rm
Error}=\frac{1}{N}\sum_{N}\left[\sqrt{\frac{1}{n_{x}n_{y}}\sum_{x}\sum_{y}\left(\phi^{\rm
pred}(x,y)-\phi^{\rm true}(x,y)\right)^{2}}\right],$ (15)

where the field $\phi$ is vorticity or velocity, $n_{x},n_{y}$ are the number
of grid points in $x$ and $y$ spatial directions, and $N$ is the number of
test samples.

##  5 Numerical Experiments

The time stepping threshold of numerical schemes for solving partial
differential equations are well understood from their stability analysis.
However, an analogous threshold for data driven models can not be derived from
similar mathematics. The idea behind this section is (a) to experiment with
different time steps to understand their relation with the accumulated errors,
and (b) to estimate the amount of data that needs to be generated from PR-DNS
model which in itself is a computationally extensive task.

###  5.1 Required number of samples

As each simulation of PR-DNS takes substantial time and resources to run and
generates a large amount of output data, we would like to assess the smallest
number of samples required to train a sufficiently accurate surrogate model.
To that end, we train three FNO models with the vorticity field using $1000$,
$4000$, and $8000$ samples, respectively. All three models are trained with
$10$ time steps ($t=0$ to $t=0.18\,t_{c}$) as input and $40$ time steps
($t=0.20\,t_{c}$ to $t=0.98\,t_{c}$) as output. We calculate the errors on the
same $500$ test samples over the prediction time ($t=0.20\,t_{c}$ to
$t=0.98\,t_{c}$) and plot them in Fig. 2.

Figure 2: Performance of models trained on $1000$, $4000$, and $8000$ samples
as input: ${\rm L}_{2}$ error calculated on $500$ test samples over the
prediction time ($t=0.20\,t_{c}$ to $t=0.98\,t_{c}$)

It is observed from Fig. 2 that the errors of the model trained with $1000$
samples grows faster with respect to time than the other two models. The
errors for models trained with $4000$ and $8000$ samples do not have
significant difference in terms of prediction errors. Thus, in the subsequent
tests we will only evaluate models trained with $1000$ and $4000$ samples as
this will save significant training data generation costs for the PR-DNS
simulations.

###  5.2 Choosing the hyperparameters

Similar to most machine learning models, the Fourier neural operator’s
performance is dependent on the choice of the so-called hyperparameters. The
hyperparameters here consist of the number of mode, neural net width, learning
rate, batch size, scheduler step size, scheduler $\gamma$, etc. In this
section, we calculate the training mean squared error (MSE) with $4000$
training samples under a few batch sizes, widths, scheduler steps, and
scheduler $\gamma$s,learning rate, and number of modes. Fig. 4 shows that the
MSE with a width 10 is larger than those using other hyperparameter settings
(due to its use of smaller neural nets) and with four modes (due to inadequate
number of Fourier modes to represent a function). The training time for widths
20 and 40 are similar. As the width of 40 requires more floating point
operations, we select 20 as the default width for the remaining
investigations.

Figure 3: Mean squared error for models trained with 4000 input samples and
different hyperparameter settings. The order of hyperparameters from left to
right is width, batch size, scheduler step, scheduler $\gamma$, learning rate,
and number of modes respectively.

###  5.3 Varying input and output time windows

Each sample of doubly periodic shear layer consists of $51$ time snapshots
spreading uniformly over unit convection time in steps of $0.02$ convection
times. As the objective of the study is to accelerate Navier-Stokes solvers, a
model with maximum output and minimum input windows is more beneficial
computationally. In the meantime, a surrogate model with low prediction errors
for longer output times is preferential.

The data set contains 51 snapshots sampled over unit convection time, which
could be split into input and output windows for FNOs in many ways, however,
we chose (a) 10 input, 40 output, (b) 20 input, 30 output, and (c) 10 inputs,
10 outputs to understand the trends in their respective errors. Thus, we test
the error for three input-output time windows:

  * •

input $t=0$ to $t=0.18\,t_{c}$ (10 snapshots); output $t=0.20\,t_{c}$ to
$t=0.38\,t_{c}$ (10 snapshots)

  * •

input $t=0$ to $t=0.18\,t_{c}$ (10 snapshots); output $t=0.20\,t_{c}$ to
$t=0.98\,t_{c}$ (40 snapshots)

  * •

input $t=0$ to $t=0.38\,t_{c}$ (20 snapshots); output $t=0.40\,t_{c}$ to
$t=0.98\,t_{c}$ (30 snapshots).

This time refinement and input-output split with acceptable errors is expected
to depend on the flow conditions too. For Reynolds number $1000$ we find that
${\rm L}_{2}$ errors are the lowest for 10 inputs and 10 outputs (see Fig. 4).
The plot also shows that the surrogate models accumulate errors for longer
prediction times – a behaviour which is expected to restrict the longest
possible prediction times for surrogate models. This experiment reveals that
the errors grow if the input windows are much smaller than the output windows.
Thus a model with fewer input snapshots and more output snapshots is expected
to perform poorly in comparison to similar number of input output snapshots.

Next, we evaluate the errors for input-output time windows different from
those for which the model was trained on. We train the model for $t=0$ to
$t=0.18\,t_{c}$ as input and $t=0.20\,t_{c}$ to $t=0.38\,t_{c}$ as output.
Once the model is obtained, we change the input-output time windows to

  * •

input $t=0.20\,t_{c}$ to $t=0.38\,t_{c}$, output $t=0.40\,t_{c}$ to
$t=0.58\,t_{c}$

  * •

input $t=0.40\,t_{c}$ to $t=0.58\,t_{c}$, output $t=0.40\,t_{c}$ to
$t=0.58\,t_{c}$

  * •

input $t=0.40\,t_{c}$ to $t=0.58\,t_{c}$, output $t=0.60\,t_{c}$ to
$t=0.78\,t_{c}$.

Thus, we are essentially computing errors for input-output windows moved ahead
in time by $0.2\,t_{c}$. In Fig. 5, at the top we plot ${\rm L}_{2}$ errors
with the ground-truth as an input, whereas the bottom figure takes the
preceding $0.2\,t_{c}$ interval’s output as its input. It can be seen that
when the input is the ground truth the errors remain at relatively low levels.
However while using previous prediction results as input, errors grow more
rapidly. This is expected as inaccurate input amplifies output errors for both
numerical simulations and ML models.

Figure 4: ${\rm L}_{2}$ errors with different input and output time windows
for models trained with $4000$ samples. The numbers in the parentheses
indicate the number of snapshots in the corresponding data set.

Figure 5: (Top) ${\rm L}_{2}$ errors with the ground-truth as an input,
(bottom) ${\rm L}_{2}$ errors for the predictions with preceding $0.2\,t_{c}$
window as in input. The numbers in the parentheses indicate the number of
snapshots in the corresponding data set.

###  5.4 Varying the time refinement

The maximum time step per iteration for classical PDE solvers is typically
dependent on convergence and numerical stability. For instance, the Courant-
Levy-Friedrichs condition dictates the time step of explicit time integration
schemes. In contrast, the time steps of ML-based surrogate models is not bound
to the time step of PDE solvers employed to generate data. Thus in this
section we try to understand the maximum permissible time steps of FNO. We
compare ${\rm L}_{2}$ errors for two time refinements of $0.02\,t_{c}$ and
$0.04\,t_{c}$, with the input and output time window sizes both fixed to 10.
For the time refinement of $0.02\,t_{c}$, the input window contains 10 time
steps from $t=0$ to $t=0.18\,t_{c}$ and the output contains 10 time steps from
$t=0.20\,t_{c}$ to $t=0.38\,t_{c}$. For the time refinement of $0.04\,t_{c}$,
the input window contains 10 time steps from $t=0$ to $t=0.36\,t_{c}$ and the
output contains 10 time steps from $t=0.40\,t_{c}$ to $t=0.76\,t_{c}$. It is
obvious that the latter case is able to predict a longer time window thereby
saving substantial computational cost of high-fidelity PR-DNS. However, it is
important to ensure the errors remain within acceptable limits. Note that the
computational cost of numerical simulation is the same for both the cases
below ($0.4\,t_{c}$), however this study serves to answer which of the
following options is better in terms of accuracy:

  * •

Run PR-DNS from $t=0$ to $t=0.18\,t_{c}$ then use FNO to predict fields from
$t=0.20\,t_{c}$ to $t=0.38\,t_{c}$. This will be later integrated in a
workflow that runs PR-DNS again from $t=0.40\,t_{c}$ to $t=0.58\,t_{c}$ then
uses FNO to predict fields from $t=0.60\,t_{c}$ to $t=0.78\,t_{c}$.

  * •

Run PR-DNS from $t=0$ to $t=0.38\,t_{c}$ then use FNO to predict fields from
$t=0.40\,t_{c}$ to $t=0.78\,t_{c}$.

The computational cost of the first approach is slightly higher as one needs
to perform two FNO predictions, however, it should be noted that the cost of
prediction from FNO is only a small fraction of the entire workflow. From Fig.
6 it is seen that the first option accumulates an ${\rm L}_{2}$ error of
$0.018$ till $t=0.18\,t_{c}$ whereas than the second option accumulates an
${\rm L}_{2}$ error of $0.064$ till $t=0.78\,t_{c}$. If the first option is
integrated in a workflow that runs PR-DNS again from $t=0.40\,t_{c}$ to
$t=0.58\,t_{c}$ and performs an FNO prediction from $t=0.60\,t_{c}$ to
$t=0.78\,t_{c}$, the FNO prediction can be reasonably assumed to accumulate
similar ${\rm L}_{2}$ errors. Thus, the net error at $t=0.78\,t_{c}$ would be
around $0.036$ which is smaller than the first approach’s $0.064$. Therefore,
the second approach with smaller intervals but higher time refinement seems to
be the better one.

Figure 6: Varying the time refinement of input-output windows for different
sample sizes. The numbers in the parentheses indicate the number of snapshots
in the corresponding data set.

###  5.5 Vorticity versus velocity formulation

Figure 7: Comparison of ${\rm L}_{2}$ errors on 500 test samples: vorticity
data set versus velocity data set for 4000 training samples.

The previous experiments are conducted with vorticity field as input and
output data. This approach is advantageous for building surrogate models for
incompressible two-dimensional flows because

  * •

One only has to deal with one variable i.e., the non-zero component of
vorticity vector ${\bm{\omega}}=(0,0,\omega_{z})$.

  * •

The continuity equation is identically satisfied in the vorticity formulation
as one can define a streamfunction $\psi$ such that
$\omega_{z}=-\nabla^{2}\psi$ which automatically produces a divergence-free
velocity field.

It should be noted that the full-fledged PR-DNS model solves the velocity
formulation of the three-dimensional Navier-Stokes equations. Thus, it is
important to assess the type and magnitude of errors when dealing with
velocity data sets.

In this section, we calculate ${\rm L}_{2}$ errors when the model training and
testing is performed with velocity fields and compare them with the vorticity
errors. It can be observed from Fig. 7 that the ${\rm L}_{2}$ error when using
the velocity data set is lower than that of vorticity. Further investigation
is required to find the reason for this behavior but it could be an artifact
of availability of twice the volume of training data, i.e., two velocity
fields $u_{x},u_{y}$ are available as as opposed to a single field of
$\omega_{z}$ in the vorticity formulation. However, on the flip side a
prediction from velocity data set does not guarantee a divergence-free
velocity field.

Figure 8: Comparison of vorticity fields of a test sample from training on
velocity versus vorticity data. The vorticity field for the velocity data set
is found by processing the velocity data.

We visualize the vorticity fields obtained from various experiments in Fig. 8.
The true vorticity field at the top is the target ground truth. The middle
plot is the vorticity field obtained from the model trained on vorticity data
set. The bottom plot is the vorticity field estimated by processing the
velocity data obtained from the model trained on the velocity data set. It can
be seen that the vortices in the target field convect away from the top right
corner without any appreciable structural difference. However, a loss in the
braids of the vortices is observed in the middle plot in addition to a noisy
vorticity field. The velocity data set performs better at preserving the
braids but shows noise levels similar to the middle plot. We find that the
errors in the continuity equation are of $\mathcal{O}(10^{-3})$ (see Fig. 9)
when using velocity data set. This will be addressed in the future by
incorporating the continuity equation in the loss function.

Figure 9: Visualization of the deviation from divergence-free condition for
fields obtained from training on velocity data set.

##  6 Conclusions and Future Directions

We have studied the performance of FNO for doubly periodic shear layer at a
relatively mild Reynolds number of 1000. We have found that the lowest
prediction errors are obtained when the input-output time intervals are
smaller than 0.2 convection times. We have also seen that the velocity data
set performs slightly better than the vorticity data, but, the velocity field
obtained is not entirely divergence-free. Unlike traditional numerical methods
at under-resolved grid refinements, FNOs do not give rise to any spurious
nonphysical vortices. Neglecting the model training cost, the velocity
prediction from FNO takes a fraction of the time compared to the numerical PDE
solver. For example, the inference from FNO takes 0.21 seconds on Nvidia V100
GPU compared to 0.65 seconds for LBM solver on AMD EPYC 7402P CPU for a mesh
of 64x64 over $0.2t_{c}$. Thus, by using FNO as the surrogate model for
velocity fields we save the cost of the numerical PDE solver which varies from
$70-90\%$ (depending on grid size and number of cores used) of the total cost.

The data from our full-fledged PR-DNS model is three dimensional and
corresponds to a Reynolds number of approximately 10,000. Thus, in the future
we plan to extend FNOs to predict three dimensional velocity fields at a
higher Reynolds number. This will require more data in terms of finer temporal
and spatial resolutions as well as number of samples. There are other neural
operators such as Clifford 27, Markov 28, and DeepONet 29 that could also be
used if they exhibit better performance or yield lower errors. As this study
focuses on only one component of the overall PR-DNS model (i.e., substitution
of velocity solver with FNO), it remains to be seen how the other components
of the model are impacted by it. However, we note that the errors in FNO’s
prediction are similar to numerical diffusion of Navier-Stokes solvers at
coarse grid resolution in that both manifest as an enhanced diffusion in the
vorticity field. This will certainly impact the other variables behavior and
will be investigated in future studies. Another outstanding question that
needs to be answered by future studies is the identification of components of
the PR-DNS model other than the Navier-Stokes equations that may benefit if
modeled via machine learning methods. Although in this work we have evaluated
the errors when replacing of the Navier-Stokes solver, our future
investigations will assess judicious replacement of thermodynamical governing
equations and/or particle propagation with ML-based surrogates as well.
Finally, incorporating governing PDEs in the loss function will lead to
smaller errors and is an avenue for further research.

*Acknowledgments The authors gratefully acknowledge financial support from the Laboratory Directed Research and Development program at Brookhaven National Laboratory. We also acknowledge administrative support from the Sustainable Research Pathways (SRP) Program for JY who was SRP student participant at Brookhaven National Lab.

*Financial disclosure

This work is supported by Laboratory Directed Research and Development program
at Brookhaven National Laboratory, which is sponsored by the US Department of
Energy, Office of Science, under Contract Number DE-SC0012704.

*Conflict of interest

The authors declare no potential conflict of interests.

## References

  * 1 Zheng Gao, Yangang Liu, Xiaolin Li, Chunsong Lu, J. Geophys. Res. Atmos. 123(4), 2194 (2018).
  * 2 Lulin Xue, Sudarsan Bera, Sisi Chen, Harish Choudhary, Shivsai Dixit, Wojciech W. Grabowski, Sandeep Jayakumar, Steven Krueger, Gayatri Kulkarni, Sonia Lasher-Trapp, Holly Mallinson, Thara Prabhakaran, Shin ichiro Shima, Bull. Am. Meteorol. Soc. 103(5), E1413 (2022), https://doi.org/10.1175/BAMS-D-22-0018.1.
  * 3 Benjamin Peherstorfer, Karen Willcox, Max Gunzburger, SIAM Review 60(3), 550 (2018), https://doi.org/10.1137/16M1082469.
  * 4 Nikos Spyropoulos, Dimitris Vlastos, George Papadakis, John M. Prospathopoulos, Vasilis A. Riziotis, J. Phys.: Conf. 2265(3) (2022), https://dx.doi.org/10.1088/1742-6596/2265/3/032008.
  * 5 Dmitrii Kochkov, Jamie A Smith, Ayya Alieva, Qing Wang, Michael P Brenner, Stephan Hoyer, Proc. Natl. Acad. Sci. 118(21), e2101784118 (2021).
  * 6 Karthik Duraisamy, Gianluca Iaccarino, Heng Xiao, Annu. Rev. Fluid Mech. 51, 357 (2019).
  * 7 J Nathan Kutz, Steven L Brunton, Bingni W Brunton, Joshua L Proctor, Dynamic mode decomposition: data-driven modeling of complex systems, (SIAM2016).
  * 8 Steven L Brunton, Joshua L Proctor, J Nathan Kutz, Proc. Natl. Acad. Sci. 113(15), 3932 (2016).
  * 9 Jean-Michel Tucny, Mihir Durve, Andrea Montessori, Sauro Succi, arXiv preprint arXiv:2305.06222 (2023).
  * 10 Tianle Wang, Bigeng Wang, Elisha Siddiqui Matekole, Mohammad Atif, in Smoky Mountains Computational Sciences and Engineering Conference, pp. 345–365.
  * 11 Tianping Chen, Hong Chen, IEEE Transactions on Neural Networks 6(4), 911 (1995).
  * 12 Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, arXiv:2010.08895 (2020).
  * 13 Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, George Em Karniadakis, Nat Mach Intell 3, 218 (2021), https://doi.org/10.1038/s42256-021-00302-5.
  * 14 Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, Journal of Machine Learning Research 24(89), 1 (2023), https://www.jmlr.org/papers/volume24/21-1524/21-1524.pdf.
  * 15 Lu Lu, Xuhui Meng, Shengze Cai, Zhiping Mao, Somdatta Goswami, Zhongqiang Zhang, George Em Karniadakis, Comput. Methods Appl. Mech. Eng. 393, 114778 (2022), https://doi.org/10.1016/j.cma.2022.114778.
  * 16 Peishi Jiang, Nis Meinert, Helga Jordão, Constantin Weisser, Simon Holgate, Alexander Lavin, Björn Lütjens, Dava Newman, Haruko Wainwright, Catherine Walker, Patrick Barnard, Digital Twin Earth - Coasts: Developing a fast and physics-informed surrogate model for coastal floods via neural operators (2021), arXiv:2110.07100.
  * 17 Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, Pedram Hassanzadeh, Karthik Kashinath, Animashree Anandkumar, FourCastNet: A Global Data-driven High-resolution Weather Model using Adaptive Fourier Neural Operators (2022).
  * 18 David L Brown, J. Comput. Phys. 122(1), 165 (1995).
  * 19 Mohammad Atif, Praveen Kumar Kolluru, Santosh Ansumali, Phys. Rev. E 106, 055307 (2022).
  * 20 Mohammad Atif, Praveen Kumar Kolluru, Chakradhar Thantanapally, Santosh Ansumali, Phys. Rev. Lett. 119, 240602 (2017).
  * 21 Seyed Ali Hosseini, Mohammad Atif, Santosh Ansumali, Iliya V Karlin, Comput. Fluids , 105884 (2023).
  * 22 Maziar Raissi, Paris Perdikaris, George E Karniadakis, J. Comput. Phys. 378, 686 (2019).
  * 23 Thomas J Grady II, Rishi Khan, Mathias Louboutin, Ziyi Yin, Philipp A Witte, Ranveer Chandra, Russell J Hewett, Felix J Herrmann, arXiv:2204.01205 (2022).
  * 24 Zhongyi Jiang, Min Zhu, Dongzhuo Li, Qiuzi Li, Yanhua O Yuan, Lu Lu, arXiv:2303.04778 (2023).
  * 25 Steven Guan, Ko-Tsung Hsu, Parag V Chitnis, arXiv:2108.09374 (2021).
  * 26 NVIDIA Modulus . https://docs.nvidia.com/deeplearning/modulus/index.html.
  * 27 Johannes Brandstetter, Rianne van den Berg, Max Welling, Jayesh K Gupta, arXiv preprint arXiv:2209.04934 (2022).
  * 28 Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, arXiv:2106.06898v2 (2021).
  * 29 Vivek Oommen, Khemraj Shukla, Somdatta Goswami, Rémi Dingreville, George Em Karniadakis, Npj Comput. Mater. 8(1), 190 (2022).

*Supporting information

Additional supporting information may be found in the online version of the
article at the publisher’s website.

##  Appendix A Fourier Neural Operator (FNO) Learning

Informally, operator learning deals with the problem of learning mappings
between (infinite dimensional) function spaces. In particular, a neural
operator is a (suitably defined) neural network architecture designed to
approximate some operator. Within the context of differential equations, such
neural operators are sought as approximations, or surrogates, to the solution
operator of differential equations 11, 12, 13, 14. Following 14, we provide
next a brief description of Fourier Neural Operators (FNOs) for partial
differential equations (PDEs), as this is the class of neural operators we
adopt in the present study. For a thorough treatment, the reader is referred
to the original publications 12, 14.

Let $\mathcal{A}$ and $\mathcal{U}$ be two Banach spaces of functions defined,
respectively, on bounded domains $\Omega_{d_{a}}\subset\mathbb{R}^{d_{a}}$ and
$\Omega_{d_{u}}\subset\mathbb{R}^{d_{u}}$. Let

$\displaystyle\mathcal{G}:\mathcal{A}\rightarrow\mathcal{U}$ (16)

be a (possibly nonlinear) operator mapping functions $a\in\mathcal{A}$ to
functions $u\in\mathcal{U}$ and

$\displaystyle\mathcal{G}_{\theta}:\mathcal{A}\times\theta\rightarrow\mathcal{U}$
(17)

be a map parameterized by $\theta\in\mathbb{R}^{p}$ such that, for some
$\theta^{*}\in\mathbb{R}^{p}$, we have that

$\displaystyle\mathcal{G}_{\theta^{*}}\ \approx\ \mathcal{G}$ (18)

is a suitable approximation and thus $\mathcal{G}_{\theta^{*}}$ may serve as a
useful surrogate to the operator $\mathcal{G}$. In the case of PDEs, the space
$\mathcal{A}$ comprises “input” functions, such as initial conditions,
coefficients, or forcing terms, that the PDEs depend on, while the solution
operator $\mathcal{G}$ maps these inputs to the solution $u\in\mathcal{U}$
satisfying the PDEs with the given inputs. Since, in practice, solving PDEs
numerically is often a time-consuming task, fast and accurate surrogates are
frequently needed.

In 14, the authors proposed a neural operator architecture (to be described
shortly) and several classes of parameterizations, to serve as surrogates (18)
to operators of interest, such as solution operators of PDEs. They prove a
universal approximation theorem showing that their proposed neural operator,
parameterization classes of which include FNOs, can approximate any given
nonlinear continuous operator. Furthermore, the neural operator proposed in 14
is discretization-invariant, meaning that it shares the same model parameters
$\theta_{*}$ among different discretizations of the underlying function
spaces. This is advantageous since, upon training (i.e., learning the optimal
model parameters $\theta_{*}$), the resulting neural operators can then be
used with discretizations of the input space $\mathcal{A}$ and output space
$\mathcal{U}$ which are possibly different than the discretizations used when
training.

Figure 10: Architecture of a Fourier neural operator (FNO). (Figure adapted
from 14. Permission to reuse material made available under a CC-BY 4.0
license, see https://creativecommons.org/licenses/by/4.0/. Attribution
requirements are provided at https://jmlr.org/papers/v24/21-1524.html)

The architecture of a Fourier neural operator, FNO, is depicted in Figure 10.
Here, a FNO $\mathcal{G}_{\theta}:\mathcal{A}\rightarrow\mathcal{U}$ is
defined as a composition

$\displaystyle\mathcal{G}_{\theta}\equiv
Q\,\circ\,\sigma_{L}\,\circ\,\left(W_{L}+\mathcal{K}_{L}\right)\,\circ\,\sigma_{L-1}\,\circ\,\left(W_{L-1}+\mathcal{K}_{L-1}\right)$
(19) $\displaystyle\,\circ\,\
\ldots\,\circ\,\sigma_{1}\,\circ\,\left(W_{1}+\mathcal{K}_{1}\right)\,\circ\,P$

of linear operators $\\{W_{l}+\mathcal{K}_{l}\\}_{l=1}^{L}$ and nonlinear
activation functions $\\{\sigma_{l}\\}_{l=1}^{L}$ for some given number $L$ of
layers in the neural network, and where:

  1. 1.

First, the input $a\in\mathcal{A}$ is lifted to a higher dimensional
representation

$\displaystyle v_{0}$ $\displaystyle=$ $\displaystyle P\circ a$ (20)

by a map $P$, often parameterized by a shallow neural network $P_{\theta_{P}}$
with parameters $\theta_{P}$.

  2. 2.

The above lifting operation is followed by a number $L$ of iterative updates
$v_{0}\mapsto v_{1}\mapsto\ldots\mapsto v_{L}$, with

$\displaystyle v_{l}$ $\displaystyle=$
$\displaystyle\sigma_{l}\circ\left(W_{l}+\mathcal{K}_{l}\right)\circ
v_{l-1},\qquad l=1,\ldots,L,$ (21)

where, for each $l=1,\ldots,L$, $W_{l}$ is a linear operator parametrized by a
neural network $W_{\theta_{W}}$ with parameters $\theta_{W}$ and

$\displaystyle\mathcal{K}_{l}$ $\displaystyle\equiv$
$\displaystyle\mathcal{F}^{-1}\circ R\circ\mathcal{F},$ (22)

where $\mathcal{F}$ denotes the Fourier transform to frequency space and
$\mathcal{F}^{-1}$ its inverse (transforming to physical space) and $R$ is an
operator filtering frequencies, which is parameterized by a neural network
$R_{\theta_{R}}$ with parameters $\theta_{R}$.

  3. 3.

Finally, the output is

$\displaystyle u$ $\displaystyle=$ $\displaystyle Q\circ v_{L},$ (23)

where $Q$ is a map that projects $v_{L}$ to $u\in\mathcal{U}$, and is
parameterized by a neural network $Q_{\theta_{Q}}$ with parameters
$\theta_{Q}$.

The parameters of the FNO $\mathcal{G}_{\theta}$ from (19) are thus given by
$\theta=(\theta_{P},\theta_{W},\theta_{R},\theta_{Q})$.