# Differentiable physics-enabled closure modeling for Burgers’ turbulence

Varun Shankar  Carnegie Mellon University  Vedant Puri  Carnegie Mellon
University  Ramesh Balakrishnan  Argonne National Laboratory  Romit Maulik
Argonne National Laboratory<EMAIL_ADDRESS>Venkatasubramanian Viswanathan
Carnegie Mellon University<EMAIL_ADDRESS>
###### Abstract

Data-driven turbulence modeling is experiencing a surge in interest following
algorithmic and hardware developments in the data sciences. We discuss an
approach using the differentiable physics paradigm that combines known physics
with machine learning to develop closure models for Burgers’ turbulence. We
consider the 1D Burgers system as a prototypical test problem for modeling the
unresolved terms in advection-dominated turbulence problems. We train a series
of models that incorporate varying degrees of physical assumptions on an a
posteriori loss function to test the efficacy of models across a range of
system parameters, including viscosity, time, and grid resolution. We find
that constraining models with inductive biases in the form of partial
differential equations that contain known physics or existing closure
approaches produces highly data-efficient, accurate, and generalizable models,
outperforming state-of-the-art baselines. Addition of structure in the form of
physics information also brings a level of interpretability to the models,
potentially offering a stepping stone to the future of closure modeling.

turbulence, burgers, subgrid-stress modeling, differentiable physics, machine
learning, neural operators

††preprint: AIP/123-QED

##  I Introduction

Simulation of fluid flow for practical applications is characterised by high
Reynolds number turbulent flows over complex geometriesDeville, Fischer, and
Mund (2002); Karniadakis (1999). The temporal and spatial scales of motion in
such flows span several orders of magnitudeMoin and Mahesh (1998). The cost of
Direct Numerical Simulations (DNS), which attempt to capture all energy
containing lengthscales in the flow, grows superlinearly with the Reynolds
number of the flowKolmogorov (1941); Coleman and Sandberg (2010); Orszag
(1970), limiting the scope of fully resolving simulations to canonical
geometries at medium Reynolds numbersMoin and Mahesh (1998).

Turbulence closure models have been developed for cases where resolving and
evolving the entirety of velocity and pressure spectra is not possible.
Techniques such as Large Eddy Simulations (LES) and Reynolds Averaged Navier-
Stokes (RANS) allow for cheaper calculations of high Reynolds number flows by
resolving only a portion of the spectrum and ‘modeling’ the rest Pope (2000).
Put simply, the goal of turbulence modeling is to find a closed system that
can describe the time-evolution of observables, such as the resolved velocity
spectrum, via a ‘closure’ map that allows one to untangle the dependence of
said observables on unresolved variablesSagaut (2006); Adams (2007); Qi and
Harlim (2022).

The problem of turbulence closure modeling has been extensively studied for
over a century, dating back to Boussinesq’s eddy viscosity hypothesis in
1877Layton (2014). More recently, machine learning (ML) based methods have
utilized for modeling fluid flow problemsBrunton, Noack, and Koumoutsakos
(2020); Vinuesa and Brunton (2021); Duraisamy, Iaccarino, and Xiao (2019);
Kochkov  _et al._ (2021); Watt-Meyer  _et al._ ; Vlachas  _et al._ (2018);
Thuerey  _et al._ (2020); Portwood  _et al._ (2021); Wang  _et al._ (2020).
The success of machine learning is built upon the ability of deep neural
networks to approximate functions in high-dimensional spaces with a number of
parameters that does not scale exponentially with the dimension of the problem
spaceWeinan (2020). As the underlying problem in closure modeling is of
finding the functional mapping between high-dimensional vector spaces, deep
neural networks can be utilized in a supervised learning framework.

However, the performance of machine learning models to turbulence modeling has
largely been unsatisfactory. A probable reason is that while neural networks
are touted to break the curse of dimensionality, the cost of attaining an
accurate model may be larger. The cost of generating high-fidelity data for
the model does scale with dimensionalityPope (2000), and even if data is
readily available, machine learning models need to be trained on large number
of samples to accurately approximate even simple functions. In other words, an
optimal closure model may be present in the space of functions that a neural
network architecture can express, but there is no clear path towards realizing
it.

Some success has been seen, however, by augmenting neural network models with
physical priors. For example, Ling et al.Ling, Kurzawski, and Templeton
(2016); Milani, Ling, and Eaton (2021) improved upon conventional machine
learning models by embedding Galilean invariance in the prediction of Reynolds
stress tensor for Reynolds Averaged Navier-Stokes simulations. Shankar et al.
Shankar  _et al._ (2020) enhanced the output of a convolutional neural network
model on homogeneous isotropic turbulence by imposing the divergence-free
constraint on the model output. Both approaches can be understood as attempts
to narrow down the space of functions that a neural network architecture can
represent.

Rackauckas et al.Rackauckas  _et al._ (2020) further develops inductive biases
in data-driven learning by directly embedding trainable parameters in a
model’s governing differential equations. Training such parameters against a
loss function that depends on the solution to the governing system requires
passing gradient information through a differential equation solve. Therefore
a solver that is compatible with automatic-differentiation toolchains in
needed. Directly backpropagating gradients is computationally expensive for
large solves, so the adjoint optimization methodChen  _et al._ (2018); Farrell
_et al._ (2013); Errico (1997); Sigmund and Maute (2013) is utilized. CFD
codes such as SU2 have utilized the adjoint method for sensitivity analysis
and shape optimizationPalacios  _et al._ (2013).

This idea is formalized into the notion of Differentiable PhysicsRamsundar,
Krishnamurthy, and Viswanathan (2021), defined as the set of scientific
avenues emerging from the marriage of state-of-the-art partial differential
equation (PDE) solvers with differentiable programming. The promise of
differentiable physics is that machine learning models can be developed to
learn specific unknown or residual physics, leading to more interpretable
models than classical ML approaches. We envision the landscape of data-driven
models on an axis of increasing inductive biases in Figure 1.

Figure 1: Data-driven models can incorporate different inductive biases, from
locality in convolutional networks, symmetries in equivariant networks, to
differential equation forms in universal differential equations. Models on the
right are constrained by known physics, such as a dynamical system or
advective and dissipative terms, leading to interpretable models where the
learning problem has been restricted to unknown or residual physics.

In this work, we comprehensively test the ability of ML techniques to develop
turbulence closure models that are resolution independent, interpretable in
the language of flow physics, and generalizeable across a range of Reynolds
numbers. We consider the 1D Burgers system as the ideal test problem for
modeling the energy spectra observed in advection-dominated turbulence
problems, and narrow our focus to the task of subgrid-stress modeling. We test
a range of differentiable physics models, that vary from a black-box approach
to training coefficients in state-of-the-art turbulence models. Experiments
are conducted over a range of Reynolds numbers and grid resolutions to assess
the generalizability of the model, and over multiple resolution-independent
architectures. We find that physics-based inductive biases are critical to the
development of data-efficient and accurate ML closures, and that these ML
closures can outperform state-of-the-art closure baselines in a wide range of
flows.

##  II Background

###  II.1 Burgers Turbulence

Turbulent flows are characterised by a competition between viscous forces,
which damp out velocity fluctuations by converting kinetic energy into heat,
and inertial forces, such as gravitation or pressure gradients. These tend to
generate and preserve velocity and pressure fluctuationsAlexander Smits
(2013). In low Reynolds number flows, viscous forces dominate and the flow is
near perfectly damped, meaning that any arbitrary fluctuation in the velocity
field would be smoothed out in no timeThomas (1953). The velocity fields
adjusts almost instantaneously to any changes in the pressure gradient that
drives the flow. At high Reynolds numbers, however, viscous forces may not be
strong enough to dampen out velocity fluctuations, and even tiny disturbances
may cause the entire flow to destabilise into turbulenceRoberts (2020). Energy
dissipation primarily happens at the smaller scales as the rate of dissipation
scales with the wavenumber squared.

We restrict the scope of this work to the task of subgrid stress modeling in
advection-dominated turbulence problems. The unforced, viscous Burgers
problem, Eq. 1, an advection-diffusion type problemFrisch and Bec (2000) that
exhibits an energy cascade similar to the Navier-Stokes system, is the perfect
test bench. The Burgers system is governed by

$\displaystyle\partial_{t}u+u\partial_{x}u=\nu\partial_{xx}u,$ (1)

where $u(x,t)$, the solution to the Burgers system, represents a time-varying
velocity field and $\nu$ is the kinematic viscosity. The 1D viscous Burgers
system is a canonical PDE that is widely used as a simplified case study for
subgrid model development and analysis Falkovich and Sreenivasan (2006);
Frisch and Bec (2000); Love (1980); LaBryer, Attar, and Vedula (2015). The
nonlinear advective term coupled with viscous dissipation leads to similar
multiscale phenomena observed in more complex flow dynamics. Eq. 1 is
deterministic which allows us to directly compare trajectories rather than
statistical averages, speeding up model training. We insulate the problem from
arbitrary effects of boundary conditions by considering Eq. 1 on the entire
real line, $\mathbb{R}$, which is approximated by periodic boundary conditions
on an interval. The Burgers equation can also be written in its conservative
form

$\partial_{t}u+\frac{1}{2}\partial_{x}u^{2}=\nu\partial_{xx}u,$ (2)

which is typically more amenable to numerical calculation.

###  II.2 Large Eddy Simulations

As high Reynolds number flows have a wide energy spectrum, DNS of such flows
is often computationally intractable, and a turbulence model is needed. Large
Eddy Simulations (LES) resolves the flow variables starting from the intertial
range up to a cut-off length, $\Delta$, that depends upon the computation
gridLesieur and Metais (1996); Boris  _et al._ (1992). Critically, an LES
calculation ignores the smallest lengthscales in the flow, which are most
computationally expensive to resolve. Filtering is accomplished by

$\displaystyle u(\bm{x},t)$
$\displaystyle=\overline{u}(\bm{x},t)+u^{\prime}(\bm{x},t)$ (3)
$\displaystyle\overline{u}(\bm{x},t)$ $\displaystyle=G_{\Delta}\star
u(\bm{x},t),$

where $G_{\Delta}$ is a low-pass filtering operator. Quantities of interest
are interpolated onto a computational grid that can resolve flow features only
up to lengthscale $\Delta$. The grid has much lower resolution than what is
needed to fully resolve the flow, and evolving flow variables on this grid
lead to saving computational resources.

Filtering the Burgers equation in this manner leads to

$\displaystyle\partial_{t}\overline{u}+\frac{1}{2}\partial_{x}\overline{u}^{2}=\nu\partial_{xx}\overline{u}-\frac{1}{2}\partial_{x}\overline{u^{\prime
2}},$ (4)

the LES equations for the Burgers system. Evolving $\overline{u}$ per Eq. 4 is
not possible as one does not have access to the subgrid stress term,
${\eta}\mathrel{\mathop{:}}=\frac{1}{2}\overline{{u}^{\prime}{u}^{\prime}}$,
and ignoring it leads to an incorrect flow field and numerical instabilities.
For the Navier-Stokes equations, ignoring subgrid stresses leads to numerical
instabilities spurious high-frequency modesMullen and Fischer (1999). The
closure problem of LES is of modeling the effect of ${\eta}$ on the dynamics
of $\overline{u}$. This is called subgrid stress modeling. As one does not
have access to the unresolved quantities, subgrid stress models are formulated
in terms of the resolved quantities.

Kolmogorov hypothesized that the finer scales are largely similarKolmogorov
(1941), and are primarily responsible for energy dissipation. One approach to
closure modeling is to take out small amounts of energy from the high-
frequency modes in $\overline{{u}}$ by applying a mild low-pass filter at
every timestep. Fisher et al.Mullen and Fischer (1999) has applied this
methodology to the Navier-Stokes equations in the Computational Fluid Dynamics
(CFD) code NEK5000Paul F. Fischer and Kerkemeier (2008), and Layton et
al.Layton, Pruett, and Rebholz (2010) has prove existence of unique strong
solutions to the resulting continuum model for the same.

Another approach is to directly model the commutative error term, ${\eta}$ in
terms of $\overline{{u}}$ with a map $\eta(\overline{u})$, and evolve the flow
variables as per Eq. 4. Hypothesizing a functional form for additive closures
of this sort is highly problem dependent. A common hypothesis is the eddy
viscosity family of modelsFrisch and Kolmogorov (1995); Schmitt (2007) that
introduce an “effective” or “eddy” viscosity, $\nu_{t}$, to the momentum
equation for enhancing dissipation. This leads to a functional form of
${\eta}=\partial_{x}\nu_{t}\partial_{x}\overline{u}$, where $\nu_{t}$ could be
a fixed scalar based on known flow physics and computational grid, or be
defined by a transport equation. In § IV we consider several such models and
apply them to the Burgers problem.

##  III Method

###  III.1 Numerical methods

We compute solutions to the Burgers equation via a pseudospectral method Boyd,
Marilyn, and Eliot (2000). The periodic 1D domain is discretized with a
regular grid of $N$ grid points. The periodic nature of the solutions allow
for expansion in a Fourier basis, as per

$u(x,t)=\sum_{k=-N/2}^{N/2}\hat{u}_{k}(t)e^{ikx}.$ (5)

The linear term of the equation can be computed pointwise in Fourier space

$\partial_{xx}\hat{u}(k,t)=-k^{2}\hat{u}(k,t),$ (6)

however the nonlinear term introduces interactions between all the wavemodes
and thus requires more care for efficient computation. Therefore, it is more
suitable to compute the quadratic term $u^{2}$ pointwise in real space at the
$N$ collocation points, then compute the spatial derivative in Fourier space.

The spatial discretization allows for a method of lines approach to time
evolution. The discrete Fourier modes from the spectral basis expansion
produces a system of $N$ ordinary differential equations (ODEs) that can be
solved with standard ODE integration schemes. We choose the Tsitouras adaptive
5th order scheme Tsitouras (2011).

###  III.2 Problem statement

Accurate simulation of the governing Burgers equation imposes certain
restrictions on the temporal and spatial discretization. As we use an adaptive
time-stepping scheme with fixed error tolerance, we focus on errors arising
from the spatial discretization. The spatial grid should be sufficiently dense
to resolve the smallest dissaptive length scales in the flow, typically much
smaller than the characteristic length scales of the flow.

We wish to approximate solutions to the Burgers equation on a considerably
coarser grid of $M$ points, by exploring a family of data-driven models
enabled by algorithmic differentiation of the temporal and spatial
discretization schemes. Each of the models is described by a PDE or set of
PDEs, where certain terms in the equations are computed by deep neural
networks. Thus, we enable significant interpretability of the models as
network outputs can be readily identified as components of traditional closure
modeling approaches.

###  III.3 Models

In this section, we provide a description of the various models that we use in
our experiments. First, we make clear our notational conventions.

$u(x,t)$ is the spatiotemporal velocity field we are modeling.
$\overline{u}(x,t)$ is the filtered velocity field on the LES grid. We drop
the parentheses for brevity. $\eta$ represents an unknown closure variable we
are modeling. The exact description of this variable differs across models. In
some models, an additional transport equation is introduced to govern the
dynamics of the closure variable. The $\overset{\bullet}{}$ and the $\nabla$
symbols correspond to the time derivative and spatial derivative of the field
respectively. A subscript 0 represents the initial condition of the field.

Neural networks are indicated with subscript $\theta$. These functions are
parameterized by a large number of weights with functional form described in
the subsequent section. Network inputs and outputs are written as:

$c,d=f_{\theta}(a,b;e),$ (7)

where $a$ and $b$ are explicit inputs to the network, and $c$ and $d$ are
outputs of the network with implicit dependence on variable $e$. For example,
$\eta=f_{\theta}(u;x)$ represents a network that takes the spatially varying
velocity field as input and produces a spatially varying $\eta$ field as
output.

The ground truth equation, or none model, is given by:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}$
(8)

Since there is no closure term or trainable networks, we use this as a
baseline comparison for our learned models.

We also compare our models to a couple baseline closure models Maulik and San
(2018), the constant Smagorinsky model (further denoted as smag-const), and
the dynamic Smagorinsky model (further denoted as smag-dyn). The constant
Smagorinsky modelSmagorinsky (1963) is given by:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$
(9) $\nu_{t}=(C_{s}\delta x)^{2}|\nabla\overline{u}|,$ (10)

where $\nu_{t}$ is the eddy viscosity, $C_{s}$ is the Smagorinsky constant,
and $\delta x$ is the grid spacing. While there are no trainable networks in
this model, $C_{s}$ can be optimized via hand tuning or gradient descent.

The dynamic Smagorinsky model Germano  _et al._ (1991); Lilly (1992) is
similar to the constant model, however the term $(C_{s}\delta x)^{2}$ is
determined from the state of the system. Studies have shown that $C_{s}$ is
not constant and may be variable with different flow characteristics
Smagorinsky (1993). The model uses a second test filter, denoted by $\sim$
with filter scale $\tilde{\delta}x$. We fix the filter ratio
$\kappa=\tilde{\delta}x/\delta x=2$. Specifically,

$(C_{s}\delta x)^{2}=\frac{\langle HM\rangle}{\langle M^{2}\rangle},$ (11)
$H=\nabla(\widetilde{\overline{u}}^{2}/2)-\nabla(\widetilde{\overline{u}^{2}}/2),$
(12)
$M=\kappa^{2}\nabla(|\nabla\widetilde{\overline{u}}|\nabla\widetilde{\overline{u}})-\nabla(\widetilde{|\nabla\overline{u}|\nabla\overline{u}}),$
(13)

where $\langle\rangle$ denotes averaging over the spatial domain.

The rest of the models in this section contain trainable networks in
functional forms that rely on varying degrees of assumptions typically used in
closure analysis. At one extreme, we have model resnet, which uses the
classical ResNet architecture He  _et al._ (2015) to predict the coarse-
grained flow field at the next time step:

$(\overline{u}_{t+\Delta t},\eta_{t+\Delta
t})=(\overline{u}_{t},\eta_{t})+f_{\theta}(\overline{u}_{t},\eta_{t},\nu;x)\Delta
t,$ (14) $\eta_{0}=g_{\theta}(\overline{u}_{0};x),$ (15)

where $\Delta t$ is a coarse time step. This model makes no assumptions about
the physical description of the system, including the fact that it may be
described as a PDE. The closure variable $\eta$ is used only to lift the model
into a higher-dimensional space governed by the neural network. The network
may choose to learn $\eta_{t}=0$, or it may leverage $\eta$ in the computation
of $\overline{u}_{t}$. Projecting or transforming an input into a higher-
dimensional subspace is a common practice in ML, where a large latent space of
hidden channels eases learning. This is a core component of augmented neural
ODEs Dupont, Doucet, and Teh (2019), which have shown to generalize better
than standard neural ODEs.

anode, an augmented neural ODE, is a small modification to resnet that assumes
an ODE inductive bias. The model is given by:

$\dot{\overline{u}},\dot{\eta}=f_{\theta}(\overline{u},\eta,\nu;x)$ (16)
$\eta_{0}=g_{\theta}(\overline{u}_{0};x),$ (17)

where $f_{\theta}$ explicitly describes the RHS of an ODE governing the
discretized flow field. The dynamics of $\overline{u}$ and $\eta$ are
determined purely by the network $f_{\theta}$ with no additional terms. This
involves only a small change in the code of resnet, changing the solver from
explicit Euler to an adaptive scheme.

Filtering the governing equations as in Eq. 4 introduces an additional closure
term to the equations. Model direct models the entirety of this term directly
as an additive correction:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\eta$
(18) $\eta=f_{\theta}(\nabla\overline{u},\nu;x).$ (19)

Model transport-I leverages the fact that the closure term can be written as
the divergence of an unknown subgrid stress. Here, $\eta$ represents this
stress analogue in 1D, which is additionally transported with another
dynamical equation:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla\eta$
(20)
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta$
(21) $\alpha,\beta=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu;x)$ (22)
$\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x).$ (23)

The coefficients of the dissipative and convective terms in the $\eta$
transport equation are given by neural networks.

Model transport-II uses the Boussinesq hypothesis that the subgrid stress has
a linear relationship with the strain in the flow. $\eta$ represents an eddy
viscosity analogue:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\eta\nabla\overline{u})$
(24)
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta$
(25) $\alpha,\beta=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu,\delta x;x)$
(26) $\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x).$ (27)

Model transport-I-p is nearly identical to transport-I, however the $\eta$
transport equation has an additional production term $\gamma$:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla\eta$
(28)
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta+\gamma$
(29) $\alpha,\beta,\gamma=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu;x)$
(30) $\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x).$ (31)

Model transport-II-p is equivalent to transport-II except for the additional
production term:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\eta\nabla\overline{u})$
(32)
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta+\gamma$
(33) $\alpha,\beta,\gamma=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu;x)$
(34) $\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x).$ (35)

The last model, smag-ml, attempts to improve on the standard constant
Smagorinsky model by learning a spatially varying Smagorinsky constant,
computed via neural network:

$\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$
(36) $\nu_{t}=(\eta\delta x)^{2}|\nabla\overline{u}|$ (37)
$\eta=f_{\theta}(\nabla\overline{u},\nu,\delta x;x).$ (38)

All of the models used in our analysis are tabulated in Table 5 of the
appendix.

###  III.4 Neural network architecture

To leverage our choice of pseudospectral discretization, we use Fourier Neural
Operators (FNO) as our network architecture Li  _et al._ (2020). All functions
with subscript $\theta$ in the previous section correspond to an FNO
architecture. FNOs combine pointwise in real space and pointwise in Fourier
space linear operations with pointwise nonlinearities in real space to enable
learning of a large suite of resolution-invariant function classes, including
linear and nonlinear operators.

Given an input feature $h^{n}\in\mathbb{R}^{N\times c_{in}}$, where $c_{in}$
is the number of channels or scalar fields, an FNO layer is defined as:

$h^{n+1}=\sigma(Bh^{n}+\mathcal{F}^{-1}(W_{k}\mathcal{F}(h^{n})_{k})),$ (39)

where $B\in\mathbb{R}^{c_{in}\times c_{out}}$ is a 2-dimensional tensor that
is contracted along $c_{in}$ identically pointwise in real space and
$W\in\mathbb{R}^{k_{max}\times c_{in}\times c_{out}}$ is a 3-dimensional
tensor that is contracted along $c_{in}$ pointwise in Fourier space for each
wavemode $k$ up to $k_{max}$. Therefore, the Fourier transform of the input is
truncated at $k_{max}$. The output of these operations are summed and a
pointwise nonlinearity $\sigma$ is applied. The parameters of the layer are
given by the tensors $W$ and $B$. Figure 2 illustrates the algorithm for
$\eta=f_{\theta}(\nabla\overline{u},\nu,\delta x;x)$ in the smag-ml model.

Figure 2: An example schematic of the FNO neural network architecture for the
function $f_{\theta}$ in smag-ml. At each time step, the network takes as
input $\nabla\overline{u}$, $\nu$, and $\delta x$, computing the $\eta$ field
that corresponds to a spatially varying subgrid term. Within the network,
there are two FNO layers that include a local linear bias from matrix $B$ and
a transformation of the Fourier modes $\hat{h}_{k}$ with matrix $W_{k}$. These
are summed and passed through a non-linearity to the next layer.

###  III.5 Optimization

Network parameters are optimized with respect to the loss function:

$L=\sum_{i=1}^{T}\sum_{j=1}^{M}\frac{(\overline{u}_{model}(x_{j},t_{i})-\overline{u}_{DNS}(x_{j},t_{i}))^{2}}{M*T}$
(40) $\hat{\overline{u}}_{DNS}=G\star\hat{u}_{DNS},$ (41)

where $u_{DNS}$ is the ground truth velocity field resolved on a grid of $N$
points, $\overline{u}_{DNS}$ is filtered using sharp low-pass filter $G$ with
cutoff $k=M/2$, $x_{j}$ correspond to the $M$ collocation points on the coarse
grid, and $t_{i}$ correspond to the $T$ timesteps in the discretized solution
trajectory.

Optimization is achieved with the gradient-based optimizer ADAM, shown to be
effective for deep learning applications Kingma and Ba (2015). Gradients are
computed with the automatic differentiation (AD) platform available in Pytorch
via reverse-mode autodifferentiation, or backpropagation Paszke  _et al._
(2019). We show a schematic of the training procedure in Figure 3.

Figure 3: A schematic of the training procedure for learning the parameters in
the data-driven models. An initial velocity field is evolved (offline)
according to the governing equations on a DNS-resolving grid and filtered to
generate ground truth data. During the optimization loop, the same field is
filtered and evolved with the data-driven model to compute an estimate for the
filtered field. The gradients are backpropagated through the ODE solve and
spatial discretization to update the parameters.

Most of the operations in our solution algorithm consist of basic linear
operations and pointwise nonlinearities, operations that have formed the
foundation of classical multi-layer perceptron machine learning architectures.
Deep learning libraries have readily supported backpropagation through these
operations for many years and thus we do not discuss implementation details
further. However, we highlight two slightly more complex components of the
algorithm corresponding to aspects of the spatial and temporal discretization
and discuss their gradient computation.

####  III.5.1 Gradients of the Fourier transform

In reverse-mode gradient accumulation Speelpenning (1980), we wish to compute
vector-Jacobian products (VJPs) Paszke  _et al._ (2017); Griewank and Walther
(2008) of each function $f$ in the algorithm:

$\text{VJP}=\mathbf{v}^{T}f^{\prime}(\mathbf{x}),$ (42)

given the cotangent vector $\mathbf{v}$, input $\mathbf{x}$, and Jacobian
$f^{\prime}(\mathbf{x})$. The VJP can then be passed to the previous function
in the computational graph to compute the next VJP.

We leverage spectral basis expansion to compute spatial derivatives in the
flow and inside the FNO layers. This is implemented with a Discrete Fourier
Transform (DFT) operation. Since the DFT is ultimately a linear operation

$\mathcal{F}_{DFT}(\mathbf{x})=\mathbf{F}\mathbf{x},$ (43)

with DFT matrix $\mathbf{F}$, the Jacobian of the transform is simply the
matrix $\mathbf{F}$. To compute the VJP, we can equivalently perform matrix-
vector multiplication with the adjoint operator $\mathbf{F^{*}}$. Since
$\mathbf{F^{*}}=N\mathbf{F}^{-1}$, we can compute the VJP as an inverse
Fourier transform multiplied by the size of the signal. The benefit of this
approach comes from the ability to implement Fast Fourier Transform Frigo and
Johnson (2005) (FFT) algorithms during the backwards pass as well, instead of
constructing the entire DFT matrix explicitly.

####  III.5.2 Propagating gradients through time integration

Time evolution of the velocity field is accomplished by solving the ODE
resulting from spatial discretization of the field using an arbitrary
integration method. A naive implementation of AD would propagate gradients
through each of basic operations in the solver. Depending on the solver
tolerance, order, and trajectory length, this may result in a very large
number of operations, and all intermediate steps must be saved in memory for
AD to work. In our experiments, backpropagating through even moderately sized
trajectory lengths, networks, and batch sizes could exceed the 16GBs of memory
on our NVIDIA V100 GPU using naive autograd.

Therefore, a more memory efficient approach is desired. We use the adjoint
method in time, which has been proposed for use in machine learning
applications by Chen et. al Chen  _et al._ (2018). The forward ODE system is
given by the evolution of the discrete Fourier modes:

$\frac{d\hat{\mathbf{u}}_{k}(t)}{dt}=f(\hat{\mathbf{u}}_{k}(t),\theta),$ (44)

where $\theta$ represents any learnable parameters in the model. The loss can
be written as an integral:

$L=\int_{t_{0}}^{t_{T}}l(\hat{\mathbf{u}}_{k}(t))dt.$ (45)

The adjoint method produces a second ODE system that evolves the adjoint
variable $\lambda$ backwards in time, given by the equation:

$\frac{d\mathbf{\lambda}(t)}{dt}=-\lambda(t)^{T}\frac{\partial
f(\hat{\mathbf{u}}_{k}(t),\theta)}{\partial\hat{\mathbf{u}}_{k}}+\frac{\partial
l(\hat{\mathbf{u}}_{k}(t))}{\partial\hat{\mathbf{u}}_{k}},$ (46)

and an integral:

$\frac{dL}{d\theta}=-\int_{t_{T}}^{t_{0}}\lambda(t)^{T}\frac{\partial
f(\hat{\mathbf{u}}_{k}(t),\theta)}{\partial\theta}dt,$ (47)

to compute the parameter sensitivities. The RHS of the equations involve VJPs
that are efficiently evaluated with AD. The adjoint system can again be solved
with an arbitrary solver, potentially even different than the forward solver.
Specifically, we use the interpolated adjoint method, which reduces the memory
footprint by constructing an interpolation of the forward solution, which
appears in the adjoint equations. During the forward pass, the solution is
saved at a selected few timesteps and reused in the adjoint computation via
interpolation during the backwards pass.

##  IV Experiments and results

To assess the quality and practicality of our proposed models, we seek to
study model performance across a range of Burgers solutions.

###  IV.1 Viscosity generalization and temporal stability

In this experiment, we examine model generalizability with respect to initial
conditions of the velocity field and the viscosity parameter $\nu$. In
addition, we investigate stability of the models by extrapolating in time
relative to the training data.

The models are trained using the optimization procedure described in the
previous section on a ground truth dataset of 100 Burgers solutions. We vary
the initial conditions of the velocity field by randomizing the Fourier modes
of the initial field. The $k$th Fourier mode is denoted by $\hat{u}_{k}$. The
dataset is generated with the parameters tabulated in Table 1.

Parameter |  |  |  Range
---|---|---|---
$\nu$ |  |  |  $[10^{-4},10^{-3})$
$\hat{u}_{k}$ |  |  |  $\left([-0.5,0.5)+[-0.5,0.5)i\right)e^{-0.5k}$, $k=0..32$
$t$ |  |  |  $[0,3)$
$N$ |  |  |  4096
$M$ |  |  |  64
Table 1: Parameters and their distributions used to generate the training set
for the first experiment, including the viscosity $\nu$, the Fourier modes of
the initial velocity field $\hat{u}_{k}$, the timespan for integration $t$,
the DNS grid resolution $N$, and the coarse grid resolution $M$.

The neural networks in the trainable models contain 2 FNO layers with
$k_{max}=16$, 128 hidden channels, and appropriate input and output sizes. The
first layer uses a ReLU nonlinearity while the second uses the identity. The
models are trained for 500 epochs using the ADAM optimizer with learning rate
varying from $10^{-3}$ to $10^{-1}$ depending on the model. Generally, models
such as smag-ml that are more constrained by inductive biases required larger
learning rates.

Model performance is gauged using the mean-squared error of the predictions
relative to the ground truth solution on an unseen test set. We leverage two
test sets to assess generalizability across viscosities and stability in time.
The first test set consists of $25\times 25$ solutions – 25 discrete
viscosities in the range of $10^{-4.5}$ – $10^{-2}$ and 25 trajectories, or
different initial conditions, per viscosity. We compute the model error for
each prediction and plot the average ensemble error as a function of
viscosity.

The second test set contains 25 trajectories at $\nu=10^{-3}$. Here, we evolve
the trajectories between $t=[0,10)$, over 3 times longer than the training
set, to see how the models behave outside of the training regime. The average
error across all trajectories is plotted as a function of time.

Figure 4: Blue shaded regions indicate the section of the plot corresponding
to the training regime, while red shading indicates a test regime. (a) Plot of
the ensemble mean-squared error relative to the ground truth DNS for each
model as a function of viscosity. All models generally improve over the no
closure baseline except $anode$, but only $smagorinsky$ improves of the
constant Smagorinsky model. (b) Ensemble MSE for each model as a function of
time. The inclusion of a production term gives rise to anomalous behavior in
the extrapolation region.

Two non-trainable models none and smag-const provide some bounding benchmarks
for model comparison. In Figure 4(a), at the upper end, the none model with no
closure results in high error at low viscosities, where the dissipative length
scales are small, with progressively lower error at higher viscosities where
the dissipative length scales become larger. An effective closure model should
at the very least perform better than this, otherwise no additional physics
have been learned to improve over the baseline governing equations solved on a
coarse grid. At the lower end, the smag-const model is a well-known closure
model that has been shown to perform well in a wide range of viscosities. An
ideal model would produce lower loss than this baseline, indicating that the
neural closure functional form can show improvements over classical closure
models from theory.

Model performance is delineated well by model accuracy relative to the
baselines in Figure 4(a). Model anode has the highest loss, greater than the
none baseline, thus showing no improvement over the governing equations. Given
that the anode model contains no inherent physics or structural form to the
PDE, this result is unsurprising as the model must learn both the governing
dynamics and any closure approximations. Despite resnet having a similar
architecture to anode, resnet can achieve marginally better accuracy, although
it is still no better than the none baseline. Theoretically, anode should be
able to learn the same optima as resnet, since anode has the same functional
form as resnet, except solved using an adaptive time-stepper. However,
actually reaching those optima can be challenging during optimization due to
the ODE inductive bias.

Models transport-I and transport-I-p show marginal improvement over none, with
transport-I-p the better of the two. These models were particularly difficult
to train, producing very stiff ODEs that made integration difficult. This is
apparent in the fact that some of the viscosities are not shown as the
stiffness of those systems led to intractable solutions. We hypothesize that
this could be due to overfitting the training set, since all of the training
samples could be solved in a reasonable amount of time.

Model direct lies firmly between the two benchmarks. Accuracy is relatively
uniform with little dependency on viscosity. Thus, while direct consistently
improves over none at low viscosities, the advantage decreases with increasing
viscosity, becoming a disadvantage at the highest viscosities.

Models transport-II and transport-II-p have accuracies comparable to smag-
const and perform almost identically to each other within the viscosity test.
Generally, we see a consistent positive bias over the baseline closure
benchmark that diminishes only at the highest viscosities.

Model smag-ml is the only model with lower error than smag-const everywhere in
Figure 4(a). Given that the model leverages the functional form of the
classical Smagorinsky model and only adds a spatial dependency, again this
result is perhaps unsurprising. However, it is promising that the improvements
persist beyond the training region.

Overall, besides models none, transport-I, and transport-I-p, accuracy is
mostly independent of the viscosity within this range. While models have been
trained on viscosities between $10^{-4}$ and $10^{-3}$, they can be used on
viscosities well outside this range, with little to no erroneous behavior
outside the training regime.

Figure 4(b) shows model accuracy on the second test set. The blue region
indicates the time span the model was trained on, while the red region shows
the extrapolation range. The model hierarchy by accuracy is generally
preserved here, i.e. anode performs the worst, smag-ml the best, and
transport-II and transport-II-p are comparable to smag-const.

However, the stability test highlights the effect of the production term in
the models. The models largely have a monotonically decreasing slope in the
error over time, except for anode, transport-I-p, and transport-II-p, which
display aberrant behavior in the extrapolation range. The inclusion of the
production term results in improvements within the training time regime, and
in the case of transport-I and transport-I-p, the difference is significant.
Outside of this range, the production term can hinder the solution and lead to
growing errors. Deviations between transport-II and transport-II-p can be seen
at $t>5$, and error from transport-I-p grows to match transport-I at $t=10$.

Model |  |  |  Relative evaluation time
---|---|---|---
none |  |  |  1.0
anode |  |  |  2.54
resnet |  |  |  0.37
direct |  |  |  2.97
transport-I |  |  |  3.36
transport-I-p |  |  |  3.40
transport-II |  |  |  3.48
transport-II-p |  |  |  3.52
smag-ml |  |  |  3.24
smag-const |  |  |  1.36
Table 2: Table of evaluation times for the models on the second test set.
Times are normalized by the none model, which is the cheapest differential
equation-based model to compute.

Table 2 shows the model evaluation times relative to none. All of the models
except for resnet are more expensive than this baseline, which is expected
considering that the models require more operations and Fourier transforms to
compute than the governing dynamics. We found that the model evaluation times
have limited correlation with the number of operations needed at each time
step, and are more heavily influenced by the stiffness of the resulting ODE
system, given the adaptive integration scheme. Thus, resnet is much more
efficient than other models because it requires the fewest function
evaluations, coinciding with the fixed Euler time steps.

We also show training curves for this experiment in Figure 5, including the
train loss and training time as a function of epoch. We point out the train
loss of resnet in Figure 5(a), which is vastly different than the test loss in
Figure 4, different behavior from the other models tested. We believe the use
of differential equation inductive biases in the other models close the
learning gap between training and testing data distributions, promoting data-
efficient and generalizable models. Figure 5(b) clearly demonstrates the
effect of ODE stiffness. Since the number of operations per function
evaluation in each model is fixed, one might expect the training time to be
linear with epoch. However, as models like transport-I and direct evolve in
training, the curves become steeper, indicating more function evaluations
required for time integration.

Figure 5: Training curves for the models in the first experiment. (a) MSE loss
of the models, normalized by none, as a function of epoch. There is
significant variability between the models in time to convergence and
magnitude of the overall change in loss. (b) Training time as a function of
epoch. For some models, the computational cost of evaluation remains constant,
leading to linear behavior. For others, the computational cost increases with
training iteration, demonstrating non-linear behavior with a positive second
derivative.

###  IV.2 Resolution invariance

In the previous experiment, we have shown that many of the models generalize
well over a wide range of viscosities and time. In this experiment, we wish to
examine our approach with regards to the coarse-grain scale. In the previous
experiment, we fixed $M=64$, here we will vary $M$ and thus the cutoff of the
low-pass filter, $k=M/2$.

Again, we generate a training dataset of 100 Burgers solutions. Since we
desire for all the cutoffs $M/2$ to lie in the inertial region of the flow
spectrum, we slightly modify the distribution of the initial condition, the
time span for integration, and the ground truth grid density $N$, realizing a
wider inertial range. We divide the set into 4 batches with
$M=64,128,256,512$. The dataset parameters are tabulated in Table 3.

Parameter |  |  |  Range
---|---|---|---
$\nu$ |  |  |  $[10^{-4},10^{-3})$
$\hat{u}_{k}$ |  |  |  $[-0.05,0.05)+[-0.05,0.05)i$, $k=0..32$
$t$ |  |  |  $[0,6)$
$N$ |  |  |  8192
$M$ |  |  |  64, 128, 256, 512
Table 3: Parameters and their distributions used to generate the training set
for the second experiment, including the viscosity $\nu$, the Fourier modes of
the initial velocity field $\hat{u}_{k}$, the timespan for integration $t$,
the DNS grid resolution $N$, and the coarse grid resolutions $M$.

We train the two best performing models from the previous experiment,
transport-II and smag-ml, as well as two modifications to the smag-ml model to
understand the impact of certain neural network architectural choices. These
models are trained for 200 epochs, which was sufficient for convergence, again
using the ADAM optimizer with learning rate varying from $10^{-3}$ to
$10^{-1}$ depending on the model.

Specifically, we add two models, smag-ml-small and smag-ml-local. smag-ml-
small reduces the number of hidden channels in the network to 2 from 128.
Since the FNO architecture requires each of the channels to be Fourier
transformed and inverse transformed at each timestep, the size of the hidden
layer has a large impact on the computational efficiency of the model.
Therefore, we wish to see how small the model can be made without sacrificing
accuracy. smag-ml-local switches the FNO layers for standard linear layers
resulting in a purely local closure correction. The FNO is a non-local
operator given the pointwise multiplication in Fourier space. Replacement with
a linear layer gives only pointwise multiplication by a linear operator in
real space. Thus, we can compare the efficacy of a non-local functional form
with a local functional form.

To test the models, we generate a new test set that includes variation in $M$.
Here, we look at a set of $13\times 25$ trajectories – 13 discrete $M$’s in
the range of $2^{6}$-$2^{9}$, including interpolation between the training
discretizations, and 25 samples per $M$. The $\nu$ and $\hat{u}_{k}$ values
are drawn from the same distributions as the training set, and $t=[0,20)$.

Figure 6: The ensemble mean-squared error of various models normalized by the
MSE of the constant Smagorinsky model, one of the baseline we hope to improve
over. The adaptive ML Smagorinsky models can reduce the error up to 40%, with
smag-ml-small demonstrating the best performance. smag-dyn can only achieve
lower error than smag-const at high resolutions, where the inertial region is
more fully resolved.

Figure 6 shows the MSE of the models normalized by the constant Smagorinsky
baseline. Since the training dataset encompasses variations in both $\nu$ and
$M$, we train the Smagorinsky constant via gradient descent on the same
training dataset as the models to obtain the optimal constant instead of hand-
tuning. The training procedure converged clearly to a value of $C_{s}=0.4$ as
shown in Figure 7. We also include a comparison to the non-trainable smag-dyn
baseline, given that most of the models in this experiment are Smagorinsky
model derivatives.

Figure 7: Convergence of the Smagorinsky constant $C_{s}$ to 0.4 as optimized
by gradient descent on the training dataset.

Similar to the previous experiment, Figure 6 shows that model transport-II
results in greater error than smag-const, while the adaptive Smagorinsky
approaches can improve over the constant model. Even when the ideal $C_{s}$ is
fit to the data, the spatiotemporally varying models outperform the baseline.
Within the adaptive trainable Smagorinsky models, the local correction can
improve over the constant baseline at all resolutions, though not as much as
the other two models. The two non-local models generally perform better than
all others, which fits with our understanding that turbulence is a non-local
phenomenon. However, surprisingly, the smag-ml-small model ultimately shows a
universal advantage over both the dynamic and constant Smagorinsky baselines.
Comparing the training and validation losses of the models provides some
insight. The overparameterized smag-ml model overfits the training set
resulting in worse performance on the test set. While this may be ameliorated
by a larger training dataset, we do not pursue this route and leave this study
for future work. Instead, we show that we can learn an effective and
computationally efficient correction to the baseline Smagorinsky model with
limited training data, due to our differentiable physics approach to end-to-
end training of ML models.

Model |  |  |  Relative evaluation time
---|---|---|---
none |  |  |  1.0
smag-const |  |  |  1.33
smag-dyn |  |  |  2.35
smag-ml |  |  |  3.65
smag-ml-small |  |  |  3.47
smag-ml-local |  |  |  1.77
transport-II |  |  |  3.97
Table 4: Table of evaluation times for the models used in the resolution
experiment. Times are normalized by the none model, which is the cheapest to
compute.

Table 4 lists the normalized evaluation times for the models in this
experiment. While the FNO-based models are unequivocally more expensive than
any of the baselines, smag-ml-local provides an interesting compromise between
computational cost and accuracy, achieving lower error than the baselines
everywhere except for the highest resolutions, where smag-dyn improves over
smag-const.

The training curves for the models are shown in Figure 8. From Figure 8(a),
convergence was achieved in 200 epochs and from Figure 8(b), these models did
not display the same kind of non-linear behavior in the training time. smag-
ml-small was much faster to train than smag-ml, despite the comparable
evaluation times.

Figure 8: Training curves for the models in the second experiment. (a) MSE
loss of the models, normalized by none, as a function of epoch. The inductive
biases of these models lead to fast convergence. (b) Training time as a
function of epoch. Although smag-ml and smag-ml-small have comparable final
evaluation costs, the differences become significant during training, in terms
of overall optimization time.

Lastly, we show some examples of predictions from model smag-ml-small compared
to smag-const. Figure 9 displays results from $M=64$, where the cutoff
wavenumber is close to the energy containing region of the spectra. Figure 10
presents results from $M=512$, where the cutoff wavenumber is on the other
side of the inertial region close to the dissipative region. The snapshots of
the velocity field in (a) and the Smagorinsky constant and eddy viscosity in
(b) are taken at $t=10$. (c) compares the DNS and model energy spectra at
$t=10$, where the cutoff wavenumber is indicated by the model spectrum. (d)
illustrates the MSE loss of the model smag-ml-small and smag-const over time,
including the relative loss of the model to the constant Smagorinsky baseline.

Figure 9: Statistics and snapshots comparing results from closure model smag-
ml-small to a baseline and DNS with $M=64$. (a) Sample snapshot of the
velocity field from DNS, the smag-const baseline, and smag-ml-small. (b)
Sample snapshot of the closure model eddy viscosity field and spatially
varying Smagorinsky constant from smag-ml-small. (c) Snapshot of the energy
spectra of the sample from DNS and smag-ml-small. The cutoff $M/2$ lies in the
very upper region of the inertial range. (d) MSE of the model and baseline
over time. Figure 10: Statistics and snapshots comparing results from closure
model smag-ml-small to a baseline and DNS with $M=512$. (a) Sample snapshot of
the velocity field from DNS, the smag-const baseline, and smag-ml-small. (b)
Sample snapshot of the closure model eddy viscosity field and spatially
varying Smagorinsky constant from smag-ml-small. (c) Snapshot of the energy
spectra of the sample from DNS and smag-ml-small. The cutoff $M/2$ lies very
close to the dissipative range. (d) MSE of the model and baseline over time.

##  V Conclusions

In this work, we systematically investigate the efficacy of various machine
learning based closure models for Burgers turbulence using an end-to-end
learning framework. We aim to develop a better understanding of how the
functional form of the data-driven closure impacts its ability to accurately
model subgrid scales, and seek to design a data-driven closure that is broadly
applicable and ideally outperforms baseline closure models. Our models are
discretization independent and thus defined by their PDE form and neural
network hyperparameters. Models are trained in an end-to-end fashion, meaning
the entire closure PDE is solved with a pseudo-spectral method at each
iteration of the optimization procedure. Differentiable programming closes the
gap between typical a priori learning and separate a posteriori testing that
is common in ML closure modeling today. Given a differentiable solution
algorithm, models can be trained via gradient descent directly on an a
posteriori loss function. This approach allows for great flexibility in model
design. Just as neural network hyperparameters can be easily tuned and
adjusted, the PDE form of the closure can be modified to incorporate various
degrees of known physics or physical assumptions. Once the high level PDE has
been defined, the system is solved and evaluated on the loss function, the MSE
of the velocity field relative to a ground truth DNS simulation.

We asses the models on a wide range of Burgers systems to determine their
generalizability. Given the spread of Reynolds numbers encountered in flow
problems and the need for a closure model that can accommodate this range, we
test our models’ performance on unseen Burgers systems with viscosities
spanning 2.5 orders of magnitude, including interpolation of the training
viscosities and extrapolation outside of this distribution. We also check the
stability of the models by extrapolation in time.

We find that models that incorporate more physical assumptions result in lower
errors on the test set than those without such inductive biases. For example,
use of the Boussinesq hypothesis in transport-II produces significantly better
approximations than transport-I, which only models a stress analogue directly.
The resnet and anode models have no assumptions of the PDE form, meaning they
must learn all the complex dynamics of the Burgers system with limited
training data, resulting in unsuccessful models. Increasing the size of the
training set may have allowed for improved learning, however we are
specifically interested in data-efficient models given the cost of generating
ground truth DNS simulations. At the other extreme, we leverage all the
assumptions of the baseline Smagorinsky model and only apply a small
modification to learn a spatially varying Smagorinsky constant. Here, we are
able to show universal improvement over the baseline Smagorinsky models at all
viscosities.

We take our analysis further by demonstrating the resolution invariance of our
closure approach. Subgrid closures model flow scales smaller than their cutoff
wavenumber, $k=M/2$ in our case. This cutoff should lie in the inertial region
of the energy spectra, but where exactly within the region it lies may be
variable depending on the coarse grid resolution and system parameters. A
universal closure model should be able to adapt to different cutoff
wavenumbers, which impacts the range of energy scales it must model. In our
second experiment, we examine the adaptive trainable Smagorinsky models in
more detail with regards to different coarse grid resolutions. We train the
model on four resolutions ranging from cutoffs very close to the dissipative
region to cutoffs close to the energy containing region. We show that the
model can continue to outperform the constant model at all resolutions,
including interpolation between the training resolutions. To the best of our
knowledge, this is the first use of a differentiable physics based closure on
varying uniform grids.

Finally, we show the impact of two neural network architectural choices – the
size of the hidden layer in the network and the non-locality of the FNO
functional form. It is typical in machine learning approaches to create very
large networks with 10s of thousands of parameters to create very expressive
networks that can act as universal approximators. However, we have included a
great deal of inductive bias in our model simply by leveraging the Smagorinsky
form, meaning our network only has to learn a small correction to see
improvement. We demonstrate that a small network of 520 parameters is enough
to learn an effective correction, even showing lower loss than the larger
model which can overfit the small training set. Lastly, we see the impact of
the non-local FNO by testing a model with linear layers. While indeed the
local adaptive model improves over the baseline, the advantage is not as
significant as the non-local model, although its computational efficiency is
greater.

Closure modeling remains a significant challenge in the fluid dynamics
community. Machine learning is enabling a new class of data-driven models, but
their performance on arbitrary flows can be difficult to characterize. We use
the Burgers system as a test bed for developing subgrid models in a
differentiable framework, training on an a posteriori loss and testing on a
large variety of Burgers systems. We show that a simple correction to existing
closure models using data-efficient neural networks can result in better
performance.

##  VI Acknowledgements

This material is based upon work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE 1745016 awarded to VS. The
authors from CMU acknowledge the support from the Technologies for Safe and
Efficient Transportation University Transportation Center, and Mobility21, A
United States Department of Transportation National University Transportation
Center. This work was supported in part by Oracle Cloud credits and related
resources provided by the Oracle for Research program. RM is supported by the
U.S. Department of Energy, Office of Science, Office of Advanced Scientific
Computing Research, under Contract DE-AC02-06CH11357.

Table 5: Model descriptions Name |  Equations |  Time-stepper |  Layer type |  Hidden channels |  Notes
---|---|---|---|---|---
none |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}$ |  Tsit5 |  N/A |  N/A |
smag-const |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$,
$\nu_{t}=(C_{s}\delta x)^{2}|\nabla\overline{u}|$ |  Tsit5 |  N/A |  N/A |  This model contains no neural networks, but $C_{s}$ is a potentially trainable parameter.
smag-dyn |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$,
$\nu_{t}=(C_{s}\delta x)^{2}|\nabla\overline{u}|$,
$(C_{s}\delta x)^{2}=\frac{\langle HM\rangle}{\langle M^{2}\rangle}$,
$H=\nabla(\widetilde{\overline{u}}^{2}/2)-\nabla(\widetilde{\overline{u}^{2}}/2)$,
$M=\kappa^{2}\nabla(|\nabla\widetilde{\overline{u}}|\nabla\widetilde{\overline{u}})-\nabla(\widetilde{|\nabla\overline{u}|\nabla\overline{u}})$ |  Tsit5 |  N/A |  N/A |
resnet |  $(\overline{u}_{t+\Delta t},\eta_{t+\Delta t})=(\overline{u}_{t},\eta_{t})+f_{\theta}(\overline{u}_{t},\eta_{t},\nu;x)\Delta t$,
$\eta_{0}=g_{\theta}(\overline{u}_{0};x)$ |  Euler |  FNO |  128 |
anode |  $\dot{\overline{u}},\dot{\eta}=f_{\theta}(\overline{u},\eta,\nu;x)$,
$\eta_{0}=g_{\theta}(\overline{u}_{0};x)$ |  Tsit5 |  FNO |  128 |
direct |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\eta$,
$\eta=f_{\theta}(\nabla\overline{u},\nu;x)$ |  Tsit5 |  FNO |  128 |
transport-I |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla\eta$,
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta$,
$\alpha,\beta=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu;x)$,
$\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x)$ |  Tsit5 |  FNO |  128 |
transport-II |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\eta\nabla\overline{u})$,
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta$,
$\alpha,\beta=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu,\delta x;x)$,
$\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x)$ |  Tsit5 |  FNO |  128 |  The input $\delta x$ to $f_{\theta}$ is only used in the resolution study.
transport-I-p |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla\eta$,
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta+\gamma$,
$\alpha,\beta,\gamma=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu;x)$,
$\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x)$ |  Tsit5 |  FNO |  128 |
transport-II-p |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\eta\nabla\overline{u})$,
$\dot{\eta}=\alpha\cdot\nu\nabla^{2}\eta-\beta\cdot\overline{u}\nabla\eta+\gamma$,
$\alpha,\beta,\gamma=f_{\theta}(\nabla\overline{u},\nabla\eta,\nu;x)$,
$\eta_{0}=g_{\theta}(\overline{u}_{0},\nu;x)$ |  Tsit5 |  FNO |  128 |
smag-ml |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$,
$\nu_{t}=(\eta\delta x)^{2}|\nabla\overline{u}|$,
$\eta=f_{\theta}(\nabla\overline{u},\nu,\delta x;x)$ |  Tsit5 |  FNO |  128 |  The input $\delta x$ to $f_{\theta}$ is only used in the resolution study.
smag-ml-local |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$,
$\nu_{t}=(\eta\delta x)^{2}|\nabla\overline{u}|$,
$\eta=f_{\theta}(\nabla\overline{u},\nu,\delta x;x)$ |  Tsit5 |  Linear |  128 |
smag-ml-small |  $\dot{\overline{u}}=\nu\nabla^{2}\overline{u}-\overline{u}\nabla\overline{u}+\nabla(\nu_{t}\nabla\overline{u})$,
$\nu_{t}=(\eta\delta x)^{2}|\nabla\overline{u}|$,
$\eta=f_{\theta}(\nabla\overline{u},\nu,\delta x;x)$ |  Tsit5 |  FNO |  2 |

## References

  * Deville, Fischer, and Mund (2002) M. O. Deville, P. F. Fischer, and E. H. Mund,  _High-order methods for incompressible fluid flow_ , Vol. 9 (Cambridge university press, 2002).
  * Karniadakis (1999) G. E. Karniadakis, “Simulating turbulence in complex geometries,” Fluid Dynamics Research 24, 343 (1999).
  * Moin and Mahesh (1998) P. Moin and K. Mahesh, “Direct numerical simulation: A tool in turbulence research,” Annual Review of Fluid Mechanics 30, 539–578 (1998).
  * Kolmogorov (1941) A. Kolmogorov, “The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds’ Numbers,” Akademiia Nauk SSSR Doklady 30, 301–305 (1941).
  * Coleman and Sandberg (2010) G. Coleman and R. Sandberg, “A primer on direct numerical simulation of turbulence – methods, procedures and guidelines,” (2010).
  * Orszag (1970) S. A. Orszag, “Analytical theories of turbulence,” Journal of Fluid Mechanics 41, 363–386 (1970).
  * Pope (2000) S. B. Pope,  _Turbulent Flows_ (Cambridge University Press, 2000).
  * Sagaut (2006) P. Sagaut,  _Large eddy simulation for incompressible flows: an introduction_ (Springer Science & Business Media, 2006).
  * Adams (2007) N. Adams, “Mathematics of large eddy simulation of turbulent flows. by lc berselli, t. iliescu & wj layton. springer, 2006. 348 pp. isbn 987 3 540 26316 6. 74.85€,” Journal of Fluid Mechanics 582, 473–475 (2007).
  * Qi and Harlim (2022) D. Qi and J. Harlim, “Machine learning-based statistical closure models for turbulent dynamical systems,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 380, 20210205 (2022), https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2021.0205 .
  * Layton (2014) W. J. Layton, “The 1877 boussinesq conjecture : Turbulent fluctuations are dissipative on the mean flow,” (2014).
  * Brunton, Noack, and Koumoutsakos (2020) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine learning for fluid mechanics,” Annual Review of Fluid Mechanics 52, 477–508 (2020).
  * Vinuesa and Brunton (2021) R. Vinuesa and S. L. Brunton, “The potential of machine learning to enhance computational fluid dynamics,” arXiv preprint arXiv:2110.02085 (2021).
  * Duraisamy, Iaccarino, and Xiao (2019) K. Duraisamy, G. Iaccarino, and H. Xiao, “Turbulence modeling in the age of data,” Annual Review of Fluid Mechanics 51, 357–377 (2019), https://doi.org/10.1146/annurev-fluid-010518-040547 .
  * Kochkov  _et al._ (2021) D. Kochkov, J. A. Smith, A. Alieva, Q. Wang, M. P. Brenner, and S. Hoyer, “Machine learning–accelerated computational fluid dynamics,” Proceedings of the National Academy of Sciences 118 (2021).
  * (16) O. Watt-Meyer, N. D. Brenowitz, S. K. Clark, B. Henn, A. Kwa, J. McGibbon, W. A. Perkins, and C. S. Bretherton, “Correcting weather and climate models by machine learning nudged historical simulations,” Geophysical Research Letters , e2021GL092555.
  * Vlachas  _et al._ (2018) P. R. Vlachas, W. Byeon, Z. Y. Wan, T. P. Sapsis, and P. Koumoutsakos, “Data-driven forecasting of high-dimensional chaotic systems with long short-term memory networks,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474, 20170844 (2018).
  * Thuerey  _et al._ (2020) N. Thuerey, K. Weißenow, L. Prantl, and X. Hu, “Deep learning methods for Reynolds-Averaged Navier–Stokes Simulations of airfoil flows,” AIAA Journal 58, 25–36 (2020).
  * Portwood  _et al._ (2021) G. D. Portwood, B. T. Nadiga, J. A. Saenz, and D. Livescu, “Interpreting neural network models of residual scalar flux,” Journal of Fluid Mechanics 907, A23 (2021).
  * Wang  _et al._ (2020) R. Wang, K. Kashinath, M. Mustafa, A. Albert, and R. Yu, “Towards physics-informed deep learning for turbulent flow prediction,” in  _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ (Association for Computing Machinery, New York, NY, USA, 2020) p. 1457–1466.
  * Weinan (2020) E. Weinan, “Machine learning and computational mathematics,” Communications in Computational Physics 28, 1639–1670 (2020).
  * Ling, Kurzawski, and Templeton (2016) J. Ling, A. Kurzawski, and J. Templeton, “Reynolds averaged turbulence modelling using deep neural networks with embedded invariance,” Journal of Fluid Mechanics 807, 155–166 (2016).
  * Milani, Ling, and Eaton (2021) P. M. Milani, J. Ling, and J. K. Eaton, “Turbulent scalar flux in inclined jets in crossflow: counter gradient transport and deep learning modelling,” Journal of Fluid Mechanics 906, A27 (2021).
  * Shankar  _et al._ (2020) V. Shankar, G. Portwood, A. Mohan, P. Mitra, C. Rackauckas, L. Wilson, D. Schmidt, and V. Viswanathan, “Learning non-linear spatio-temporal dynamics with convolutional neural odes,” in  _Third Workshop on Machine Learning and the Physical Sciences (NeurIPS 2020)_ (2020).
  * Rackauckas  _et al._ (2020) C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, D. Skinner, A. Ramadhan, and A. Edelman, “Universal Differential Equations for Scientific Machine Learning,” arXiv e-prints , arXiv:2001.04385 (2020), arXiv:2001.04385 [cs.LG] .
  * Chen  _et al._ (2018) T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” CoRR abs/1806.07366 (2018), 1806.07366 .
  * Farrell  _et al._ (2013) P. E. Farrell, D. A. Ham, S. W. Funke, and M. E. Rognes, “Automated derivation of the adjoint of high-level transient finite element programs,” SIAM Journal on Scientific Computing 35, C369–C393 (2013).
  * Errico (1997) R. M. Errico, “What is an adjoint model?” Bulletin of the American Meteorological Society 78, 2577–2592 (1997).
  * Sigmund and Maute (2013) O. Sigmund and K. Maute, “Topology optimization approaches,” Structural and Multidisciplinary Optimization 48, 1031–1055 (2013).
  * Palacios  _et al._ (2013) F. Palacios, M. Colonno, A. Aranake, A. Campos, S. Copeland, T. Economon, A. Lonkar, T. Lukaczyk, T. Taylor, and J. Alonso, “Stanford university unstructured (su2): An open-source integrated computational environment for multi-physics simulation and design,” (2013).
  * Ramsundar, Krishnamurthy, and Viswanathan (2021) B. Ramsundar, D. Krishnamurthy, and V. Viswanathan, “Differentiable physics: A position piece,” (2021).
  * Alexander Smits (2013) I. M. Alexander Smits, “Wall-bounded turbulence,” Physics Today 66 (2013), 10.1063/PT.3.2114.
  * Thomas (1953) L. Thomas, “The stability of plane poiseuille flow,” Physical Review 91, 780 (1953).
  * Roberts (2020) M. S. Roberts, “Fluid instabilities and transition to turbulence,” in  _Computational Overview of Fluid Structure Interaction_, edited by K. Ghaedi, A. Alhusseny, A. Nasser, and N. Al-Zurf (IntechOpen, Rijeka, 2020) Chap. 2.
  * Frisch and Bec (2000) U. Frisch and J. Bec, “Burgulence,”  (2000), 10.48550/ARXIV.NLIN/0012033.
  * Falkovich and Sreenivasan (2006) G. Falkovich and K. R. Sreenivasan, “Lessons from hydrodynamic turbulence,” Tech. Rep. (Abdus Salam International Centre for Theoretical Physics, 2006).
  * Love (1980) M. Love, “Subgrid modelling studies with burgers’ equation,” Journal of Fluid Mechanics 100, 87–110 (1980).
  * LaBryer, Attar, and Vedula (2015) A. LaBryer, P. Attar, and P. Vedula, “A framework for large eddy simulation of burgers turbulence based upon spatial and temporal statistical information,” Physics of Fluids 27, 035116 (2015).
  * Lesieur and Metais (1996) M. Lesieur and O. Metais, “New trends in large-eddy simulations of turbulence,” Annual review of fluid mechanics 28, 45–82 (1996).
  * Boris  _et al._ (1992) J. P. Boris, F. F. Grinstein, E. S. Oran, and R. L. Kolbe, “New insights into large eddy simulation,” Fluid dynamics research 10, 199 (1992).
  * Mullen and Fischer (1999) J. S. Mullen and P. F. Fischer, “Filtering techniques for complex geometry fluid flows,” Communications in Numerical Methods in Engineering 15, 9–18 (1999).
  * Paul F. Fischer and Kerkemeier (2008) J. W. L. Paul F. Fischer and S. G. Kerkemeier, “nek5000 Web page,” (2008), http://nek5000.mcs.anl.gov.
  * Layton, Pruett, and Rebholz (2010) W. J. Layton, C. D. Pruett, and L. G. Rebholz, “Temporally regularized direct numerical simulation,” Applied mathematics and computation 216, 3728–3738 (2010).
  * Frisch and Kolmogorov (1995) U. Frisch and A. N. Kolmogorov,  _Turbulence: the legacy of AN Kolmogorov_ (Cambridge university press, 1995).
  * Schmitt (2007) F. G. Schmitt, “About boussinesq’s turbulent viscosity hypothesis: historical remarks and a direct evaluation of its validity,” Comptes Rendus Mécanique 335, 617–627 (2007), joseph Boussinesq, a Scientist of bygone days and present times.
  * Boyd, Marilyn, and Eliot (2000) J. Boyd, T. Marilyn, and P. Eliot, “Chebyshev and fourier spectral methods,” (2000).
  * Tsitouras (2011) C. Tsitouras, “Runge–kutta pairs of order 5(4) satisfying only the first column simplifying assumption,” Computers & Mathematics with Applications 62, 770–775 (2011).
  * Maulik and San (2018) R. Maulik and O. San, “Explicit and implicit LES closures for burgers turbulence,” Journal of Computational and Applied Mathematics 327, 12–40 (2018).
  * Smagorinsky (1963) J. Smagorinsky, “General circulation experiments with the primitive equations: I. the basic experiment,” Monthly weather review 91, 99–164 (1963).
  * Germano  _et al._ (1991) M. Germano, U. Piomelli, P. Moin, and W. H. Cabot, “A dynamic subgrid-scale eddy viscosity model,” Physics of Fluids A: Fluid Dynamics 3, 1760–1765 (1991).
  * Lilly (1992) D. K. Lilly, “A proposed modification of the germano subgrid-scale closure method,” Physics of Fluids A: Fluid Dynamics 4, 633–635 (1992).
  * Smagorinsky (1993) J. Smagorinsky, “Some historical remarks on the use of nonlinear viscosities,” Large eddy simulation of complex engineering and geophysical flows 3, 36 (1993).
  * He  _et al._ (2015) K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,”  (2015).
  * Dupont, Doucet, and Teh (2019) E. Dupont, A. Doucet, and Y. W. Teh, “Augmented neural odes,” in  _Advances in Neural Information Processing Systems_, Vol. 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019).
  * Li  _et al._ (2020) Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, “Fourier neural operator for parametric partial differential equations,” arXiv preprint arXiv:2010.08895 (2020).
  * Kingma and Ba (2015) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in  _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_, edited by Y. Bengio and Y. LeCun (2015).
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Z. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” CoRR abs/1912.01703 (2019), 1912.01703 .
  * Speelpenning (1980) B. Speelpenning, “Compiling fast partial derivatives of functions given by algorithms,”  (1980), 10.2172/5254402.
  * Paszke  _et al._ (2017) A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer, “Automatic differentiation in pytorch,” in  _NIPS 2017 Workshop on Autodiff_ (2017).
  * Griewank and Walther (2008) A. Griewank and A. Walther,  _Evaluating Derivatives_, 2nd ed. (Society for Industrial and Applied Mathematics, 2008) https://epubs.siam.org/doi/pdf/10.1137/1.9780898717761 .
  * Frigo and Johnson (2005) M. Frigo and S. G. Johnson, “The design and implementation of FFTW3,” Proceedings of the IEEE 93, 216–231 (2005), special issue on “Program Generation, Optimization, and Platform Adaptation”.


# Korenblum’s principle for Bergman spaces with radial weights

Iason Efraimidis  Departamento de Matemáticas, Universidad Autónoma de Madrid,
28049 Madrid, Spain<EMAIL_ADDRESS>,  Adrián Llinares  Department of
Mathematics and Mathematical Statistics, Umeå University, 90736 Umeå, Sweden.
Current address: Departamento de Análisis Matemático y Matemática Aplicada,
Universidad Complutense de Madrid, 28040 Madrid, Spain<EMAIL_ADDRESS>and
Dragan Vukotić  Departamento de Matemáticas, Universidad Autónoma de Madrid,
28049 Madrid, Spain<EMAIL_ADDRESS>Peter Duren (1935-2020) in memoriam

(Date: 05 April, 2024.)

###### Abstract.

We show that the Korenblum maximum (domination) principle is valid for
weighted Bergman spaces $A^{p}_{w}$ with arbitrary (non-negative and
integrable) radial weights $w$ in the case $1\leq p<\infty$. We also notice
that in every weighted Bergman space the supremum of all radii for which the
principle holds is strictly smaller than one. Under the mild additional
assumption $\liminf_{r\to 0^{+}}w(r)>0$, we show that the principle fails
whenever $0<p<1$.

###### Key words and phrases:

Weighted Bergman space, domination

###### 2010 Mathematics Subject Classification:

All authors are partially supported by PID2019-106870GB-I00 from MICINN,
Spain. The first author is supported by a María Zambrano contract, reference
number CA3/RSUE/2021-00386, from UAM and Ministerio de Universidades, Spain
(Plan de Recuperación, Transformación y Resiliencia). Second author’s work is
funded by the postdoctoral scholarship JCK22-0052 granted by The Kempe
Foundations

##  1\. Introduction

###  1.1. Korenblum domination (maximum) principle

Let $A^{p}$ denote the Bergman spaces of all analytic functions in the unit
disk ${\mathbb{D}}$ that are $p$-integrable with respect to the normalized
Lebesgue area measure $dA$. In the late 1980s, Korenblum conjectured that the
following principle should be valid: there exists a radius $c\in(0,1)$ such
that, for any two functions $f$, $g\in A^{2}$ with $|f(z)|\leq|g(z)|$ in
$A_{c}=\\{z\,:\,c\leq|z|<1\\}$, we have $\|f\|_{2}\leq\|g\|_{2}$. He proved
the statement in [8] under some additional assumptions on the zeros of $f$ and
$g$. It is quite easy to formulate a generalization of this question for any
positive $p$. Various authors produced related partial results on such
questions during the 1990s.

In 1999, Hayman [4] proved this conjecture for $p=2$, showing that
$c=1/25=0.04$ works (and hence so does any smaller positive value). Shortly
afterwards, Hinkkanen [6] extended the result for any $p\geq 1$ and
$c=0.15724$. Their methods of proof were somewhat similar. Hinkkanen relies on
one fixed estimate and therefore proved the result for the same radius,
independently of the value of $p$. Further relevant technical improvements
were obtained by Schuster [9] and also by Wang [10, 11, 12, 13] regarding the
largest possible radius $c$.

In 2018, Božin and Karapetrović [1] showed that Korenblum’s principle fails
for all $p$ with $0<p<1$ by exhibiting a universal counterexample. More
recently, Karapetrović [7] studied the problem for the more general mixed-norm
spaces $H(p,q,\alpha)$ which include the Bergman spaces with standard radial
weights as special cases when $p=q$. He was able to show that the principle
holds when $1\leq p\leq q<\infty$ and fails when $0<q<1$.

###  1.2. Description of our results

Karapetrović’s carefully written paper [7], together with other earlier works,
inspired our research and made us think of generalizing the results to Bergman
spaces with general radial weights. In fact, it turns out that the previous
results, both the positive and the negative ones, carry over with quite
similar proofs but also with some further simplifications to these general
cases.

In the present note we show that the Korenblum principle actually holds in any
weighted Bergman spaces $A^{p}_{w}$ with a non-negative and integrable radial
weight $w$ in the case $p\geq 1$ (Theorem 2). We also notice than in general,
if the principle holds, the supremum of all admissible radii is strictly
smaller than one (Proposition 1). Under the additional assumption
$\liminf_{r\to 0^{+}}w(r)>0$, we show that the principle fails whenever
$0<p<1$ (Theorem 3) and if such an assumption is omitted, it may hold even
then (Example 4).

We notice that proving these facts actually does not require changing many
details from the existing discussions in the literature, something that seems
to have remained unnoticed until now. All the necessary basic concepts are
reviewed and detailed proofs are given in the next section.

##  2\. Domination principle for weighted Bergman spaces with radial weights

###  2.1. Hardy and Bergman spaces

It is well known [2, Chapter 1] that for a subharmonic function $g$ in the
unit disk ${\mathbb{D}}$ its integral means

$\int_{0}^{2\pi}g(re^{i\theta})\frac{d\theta}{2\pi}\,,$

on the circle $\\{z\,:\,|z|=r\\}$, $0<r<1$, are increasing functions of the
radius $r$. In particular, for any function $f$ analytic in ${\mathbb{D}}$ and
$0<p<\infty$, the integral means of order $p$ of $f$:

$M_{p}(r;f)=\left(\int_{0}^{2\pi}|f(re^{i\theta})|^{p}\frac{d\theta}{2\pi}\right)^{1/p}$

are increasing functions of $r\in(0,1)$. The Hardy space $H^{p}$ is the set of
all analytic functions in ${\mathbb{D}}$ for which these means have finite
limits: $\|f\|_{H^{p}}=\lim_{r\to 1^{-}}M_{p}(r;f)<\infty$. This is not a true
norm if $0<p<1$ but the same notation is still used in this case.

Let $dA$ denote the normalized Lebesgue area measure:
$dA=\frac{1}{\pi}\,dx\,dy=\frac{1}{\pi}\,r\,dr\,d\theta$, where
$z=x+iy=re^{i\theta}\in{\mathbb{D}}$. For $p\in(0,\infty)$, the Bergman space
$A^{p}$ is defined as the set of all analytic functions $f$ in ${\mathbb{D}}$
that are $p$-integrable with respect to area measure:

$\|f\|_{p}=\left(\int_{\mathbb{D}}|f|^{p}\,dA\right)^{1/p}=\left(\int_{0}^{1}2rM_{p}^{p}(r;f)\,dr\right)^{1/p}<\infty\,.$

Again, the same comment regarding the case $0<p<1$ as above is in place. For
the theory of Bergman spaces, we refer the reader to [3] or [5].

###  2.2. Weighted Bergman spaces

In this section, the term weight will be used for a measurable (with respect
to $dA$) and non-negative function $w$ on ${\mathbb{D}}$. (Note that non-zero
constant functions belong to the space if and only if $w$ is Lebesgue
integrable.) To avoid trivial situations, we also assume that
$\int_{\mathbb{D}}w\,dA>0$; equivalently, $w$ is positive on some subset of
positive measure of the disk. For $0<p<\infty$, the weighted Bergman space
$A^{p}_{w}$ is the set of all holomorphic functions $f$ in ${\mathbb{D}}$ for
which

$\|f\|_{p,w}=\left(\int_{\mathbb{D}}|f|^{p}w\,dA\right)^{1/p}<\infty\,,$

We will say that the Korenblum domination principle holds for $A^{p}_{w}$ if
there exists $c\in(0,1)$ such that for arbitrary functions $f$, $g\in
A^{p}_{w}$ that satisfy $|f(z)|\leq|g(z)|$ whenever $c\leq|z|<1$, the
inequality $\|f\|_{p,w}\leq\|g\|_{p,w}$ holds. Any such value $c$ is said to
be an admissible radius. If no such $c$ exists, we will say that the Korenblum
principle fails for $A^{p}_{w}$.

Trivially, if $c$ is an admissible radius, any smaller positive value is also
admissible; hence, the set of all admissible radii is an interval. The
supremum of all admissible values is usually called the Korenblum radius and
is sometimes denoted by $c(p)$.

The exact value of $c(p)$ is not known for any $p\in(0,\infty)$, even for
unweighted Bergman spaces. Wang has improved various upper and lower bounds
from earlier papers in [10, 11, 12] and has also shown in [13] that
$\lim_{p\to\infty}c(p)=1$, thus answering a question posed in [6] by
Hinkkanen. It is not known whether $c(p)$ is an increasing or continuous
function of $p$. Such questions are certainly of interest but will not be
considered in the present note.

###  2.3. The Korenblum radius is smaller than one for weighted Bergman spaces

Considering the Korenblum principle and Korenblum radius also makes sense for
other spaces, such as Hardy spaces or the more general mixed norm spaces.
Every function $f\in H^{p}$ has radial limits almost everywhere and
$\|f\|_{H^{p}}$ can be computed in terms of the modulus of its radial limits
only. From here it is immediate that every $c\in(0,1)$ is an admissible
radius. Therefore the Korenblum radius is equal to $1$ for all $H^{p}$ spaces.

Our next result shows that in weighted Bergman spaces this situation is
impossible. Here we will still not assume the integrability of the weight but
we will consider only the cases when the space $A^{p}_{w}$ is non-trivial.

######  Proposition 1.

Let $0<p<\infty$ and let $w$ be a weight (not necessarily integrable) and
assume that $A^{p}_{w}\neq\\{0\\}$. If the Korenblum principle holds for such
space $A^{p}_{w}$, then the corresponding Korenblum radius must satisfy
$c(p)<1$.

###### Proof.

Assume that, on the contrary, there exists an increasing sequence of
admissible radii $(c_{n})_{n}$ such that $\lim_{n\to\infty}c_{n}=1$. Given
$f\in A^{p}_{w}\setminus\\{0\\}$, consider the functions $f_{n}$ given by

$f_{n}(z)=\frac{z}{c_{n}}f(z)\,.$

Clearly, $f_{n}\in A^{p}_{w}\setminus\\{0\\}$ and $|f(z)|\leq|f_{n}(z)|$
whenever $c_{n}\leq|z|<1$. Since $c_{n}$ is an admissible radius by
assumption, it follows that

$\int_{\mathbb{D}}|f(z)|^{p}w(z)\,dA(z)\leq\frac{1}{c_{n}^{p}}\int_{\mathbb{D}}|z|^{p}\,|f(z)|^{p}w(z)\,dA(z)\,.$

After taking the limit as $n\to\infty$, this implies that

$\int_{\mathbb{D}}(1-|z|^{p})|f(z)|^{p}w(z)\,dA(z)\leq 0\,,$

which is only possible if $|f(z)|^{p}w(z)=0$ almost everywhere in
${\mathbb{D}}$. Since $w(z)=0$ cannot be satisfied almost everywhere in
${\mathbb{D}}$ in view of our assumption on $w$, it follows that $f(z)=0$ on a
set of positive measure in ${\mathbb{D}}$, hence $f\equiv 0$, which is absurd.
This proves the claim. ∎

###  2.4. Radial weights

From now on, we only consider integrable (non-negative) weights; hence all
bounded analytic functions in ${\mathbb{D}}$ will belong to $A^{p}_{w}$.

A weight is said to be radial if $w(z)=w(|z|)$ for all $z\in{\mathbb{D}}$. Our
earlier assumptions in this case include the condition that $w(r)>0$ for all
$r$ in some subset of positive measure of $(0,1)$. In this case, the norm of a
function $f\in A^{p}_{w}$ can be computed as follows:

$\|f\|_{p,w}=\left(\int_{0}^{1}2rw(r)M_{p}^{p}(r;f)\,dr\right)^{1/p}\,,$

where $M_{p}(r;f)$ is as before.

It should be noted that $A^{p}_{w}$ may not be complete for some weights of
this type (for example, if $w(r)=0$ on $(R,1)$ for some $R\in(0,1)$, as the
reader can check after a certain amount of work, but this will not be relevant
at all for our purpose. In other words, our discussion is independent of the
completeness of the space.

###  2.5. The case $p\geq 1$

Our first result is positive, generalizing the earlier findings for $p\geq 1$.
We follow the methods of [6] and [9], explained very well in [7]. Notice that
in our considerations we only need the simplest limit of a certain integral,
without any need to compute or estimate certain values, unlike in earlier
papers. This is achieved at the expense of not giving any explicit bounds on
the value of $c$ for which the principle holds.

######  Theorem 2.

Let $1\leq p<\infty$ and let $w$ be a radial weight. Then the Korenblum
domination principle holds for the space $A^{p}_{w}$.

###### Proof.

Let $0<c<1$ (the appropriate value is to be defined later) and consider $f$,
$g\in A^{p}_{w}$ that satisfy $|f(z)|\leq|g(z)|$ whenever $c\leq|z|<1$.

Define $\omega(z)=\frac{f(z)}{g(z)}$ for $c<|z|<1$. Note that $\omega$ is
holomorphic in this annulus and without loss of generality we may assume that
$|\omega|<1$ there; otherwise $f$ and $g$ coincide and there is nothing left
to prove. Also, if $\omega$ is a constant, then $|f|$ is a constant multiple
of $|g|$ with the constant $<1$, and the statement of the Korenblum principle
follows trivially; thus, we may assume from now on that $\omega$ is not
constant.

If it is that case that

$\int_{\\{|z|<c\\}}\left(|f|^{p}-|g|^{p}\right)w\,dA\leq 0\,,$

then after integrating the positive function $\left(|g|^{p}-|f|^{p}\right)w$
over the annulus $\\{c<|z|<1\\}$, simple algebra readily yields the inequality
$\|f\|_{p,w}^{p}\leq\|g\|_{p,w}^{p}$, and the statement of the Korenblum
principle follows. Thus, from now on we may restrict our attention only to the
case when

$\int_{\\{|z|<c\\}}\left(|f|^{p}-|g|^{p}\right)w\,dA>0\,.$

For each $\rho$ in the interval $(c,1)$ there exists a point $\zeta_{\rho}$
such that

$|\zeta_{\rho}|=\rho\,,\quad|\omega(\zeta_{\rho})|=\max\\{|\omega(z)|\,:\,|z|=\rho\\}\,.$

Next, define

$\gamma(\rho)=\max\left\\{\frac{|\omega(z)-\omega(\zeta_{\rho})|}{1-|\omega(z)|^{2}}\,:\,|z|=\rho\right\\}\,.$

Note that $\gamma(\rho)>0$ for all $\rho\in(c,1)$ since $\omega$ is not
constant. It is important to stress that $\gamma$ depends on the choice of
functions $f$, $g$ that satisfy our assumptions.

We continue following Hinkkanen [6] and using some basic inequalities. To this
end, in the next step the assumption $p\geq 1$ is actually fundamental since
we will use the elementary double inequality

(1) $py^{p-1}(x-y)\leq x^{p}-y^{p}\leq px^{p-1}(x-y)\,,\quad x\,,y\geq
0\,,\quad p\geq 1\,.$

For $0<r<c<\rho<1$, using the fact that $|\omega(\zeta_{\rho})|<1$ for all $z$
with $|z|=\rho$ and the second inequality in (1), combined with the basic
triangle inequality, we obtain the following chain of estimates

$\displaystyle M_{p}^{p}(r,f)-M_{p}^{p}(r,g)$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}\left(|f(re^{i\theta})|^{p}-|g(re^{i\theta})|^{p}\right)\frac{d\theta}{2\pi}$
$\displaystyle\leq$
$\displaystyle\int_{0}^{2\pi}\left(|f(re^{i\theta})|^{p}-|\omega(\zeta_{\rho})g(re^{i\theta})|^{p}\right)\frac{d\theta}{2\pi}$
$\displaystyle\leq$
$\displaystyle\int_{0}^{2\pi}p|f(re^{i\theta})|^{p-1}\left(|f(re^{i\theta})|-|\omega(\zeta_{\rho})g(re^{i\theta})|\right)\frac{d\theta}{2\pi}$
$\displaystyle\leq$
$\displaystyle\int_{0}^{2\pi}p|f(re^{i\theta})|^{p-1}\left|f(re^{i\theta})-\omega(\zeta_{\rho})g(re^{i\theta})\right|\frac{d\theta}{2\pi}$
$\displaystyle\leq$ $\displaystyle\int_{0}^{2\pi}p|f(\rho
e^{i\theta})|^{p-1}\left|f(\rho e^{i\theta})-\omega(\zeta_{\rho})g(\rho
e^{i\theta})\right|\frac{d\theta}{2\pi}\,,$

where in the last step we have used the fact that the function
$|f|^{p-1}|f-\omega(\zeta_{\rho})g|$ is subharmonic. (To check this, just note
that its logarithm is clearly subharmonic since $p-1\geq 0$ and recall that
the exponential function of a subharmonic function is subharmonic.)

We continue the above chain of inequalities by using the definitions of
$\omega(z)$ and $\gamma(\rho)$ and the first inequality in (1):

$\displaystyle M_{p}^{p}(r,f)-M_{p}^{p}(r,g)$ $\displaystyle\leq$
$\displaystyle\int_{0}^{2\pi}p|f(\rho e^{i\theta})|^{p-1}\left(|g(\rho
e^{i\theta})|-|f(\rho e^{i\theta})|\right)\frac{\left|f(\rho
e^{i\theta})-\omega(\zeta_{\rho})g(\rho e^{i\theta})\right|}{|g(\rho
e^{i\theta})|-|f(\rho e^{i\theta})|}\frac{d\theta}{2\pi}$ $\displaystyle=$
$\displaystyle\int_{0}^{2\pi}p|f(\rho e^{i\theta})|^{p-1}\left(|g(\rho
e^{i\theta})|-|f(\rho e^{i\theta})|\right)\frac{|\omega(\rho
e^{i\theta})-\omega(\zeta_{\rho})|}{1-|\omega(\rho
e^{i\theta})|}\frac{d\theta}{2\pi}$ $\displaystyle\leq$ $\displaystyle
2\gamma(\rho)\int_{0}^{2\pi}\left(|g(\rho e^{i\theta})|^{p}-|f(\rho
e^{i\theta})|^{p}\right)\frac{d\theta}{2\pi}$ $\displaystyle=$ $\displaystyle
2\gamma(\rho)\left(M_{p}^{p}(\rho,g)-M_{p}^{p}(\rho,f)\right)\,.$

For a fixed $\rho$, integrate from $0$ to $c$ with respect to $2rw(r)\,dr$:

$\int_{0}^{c}\left(M_{p}^{p}(r,f)-M_{p}^{p}(r,g)\right)\,2rw(r)\,dr\leq
2\gamma(\rho)\left(M_{p}^{p}(\rho,g)-M_{p}^{p}(\rho,f)\right)\,\int_{0}^{c}2rw(r)\,dr\,,$

Next, divide both sides by $\gamma(\rho)$ and then integrate from $c$ to $1$
with respect to $\rho w(\rho)\,d\rho$ to get

(2) $\int_{c}^{1}\frac{\rho
w(\rho)}{\gamma(\rho)}\,d\rho\cdot\int_{\\{|z|<c\\}}\left(|f|^{p}-|g|^{p}\right)w\,dA\leq\int_{0}^{c}2\,r\,w(r)\,dr\cdot\int_{\\{c<|z|<1\\}}\left(|g|^{p}-|f|^{p}\right)w\,dA\,.$

In view of our assumption that the second factor on the left is positive, it
follows that

$\int_{c}^{1}\frac{\rho w(\rho)}{\gamma(\rho)}\,d\rho<\infty\,.$

If we can find a value $c\in(0,1)$ independent of $\gamma$ (that is,
independent of $f$ and $g$) and such that

(3) $\int_{0}^{c}2\,r\,w(r)\,dr\leq\int_{c}^{1}\frac{\rho
w(\rho)}{\gamma(\rho)}\,d\rho\,,$

then, after a cancelation, (2) will immediately yield

$\|f\|_{p,w}^{p}=\int_{\\{|z|<c\\}}|f|^{p}w\,dA+\int_{\\{c<|z|<1\\}}|f|^{p}w\,dA\leq\int_{\\{|z|<c\\}}|g|^{p}w\,dA+\int_{\\{c<|z|<1\\}}|g|^{p}w\,dA=\|g\|_{p,w}^{p}\,,$

which will prove our claim.

To complete the proof, it remains to show that (3) holds for some $c\in(0,1)$
independent of $\gamma$. On the one hand, the Monotone Convergence Theorem
implies that

$\lim_{c\to 0^{+}}\int_{0}^{c}2\,r\,w(r)\,dr=0\,.$

On the other hand, if $0<c<\frac{1}{4}$ we can use fundamental estimates due
to Schuster. Specifically, in [9, p. 3528] the following function was defined:

$F(\rho,c)=\frac{2c}{\rho}\left(1+\frac{\rho^{2}}{c}\right)\cdot\frac{1-c^{12}}{1-c^{10}}\cdot\prod_{n=1}^{5}\frac{(1+\rho^{2}c^{2n-1})(1+\rho^{-2}c^{2n+1})(1+c^{2n})^{2}}{(1+\rho^{2}c^{2n-2})(1+\rho^{-2}c^{2n})(1+c^{2n-1})^{2}}$

and it was shown that

$\gamma(\rho)\leq\frac{F(\rho,c)}{\sqrt{1-F(\rho,c)^{2}}}=H(\rho,c)\,.$

Note that $H$ is an explicit function independent of $f$ and $g$ (hence of
$\gamma$). It also has the property that

$\lim_{c\to 0^{+}}H(\rho,c)=\frac{2\rho}{1-\rho^{2}}\,.$

Recalling that $w>0$ on a set of positive measure, Fatou’s Lemma then implies
that

$0<\frac{1}{2}\int_{0}^{1}w(\rho)(1-\rho^{2})\,d\rho\leq\liminf_{c\to
0^{+}}\int_{c}^{1}\frac{\rho w(\rho)}{H(\rho,c)}\,d\rho\leq\liminf_{c\to
0^{+}}\int_{c}^{1}\frac{\rho w(\rho)}{\gamma(\rho)}\,d\rho\,.$

Thus, there exists a small enough $c>0$ such that (3) holds independently of
$\gamma$.

This shows that the Korenblum principle holds for such a radius $c$ and
$A^{p}_{w}$. ∎

###  2.6. The case $0<p<1$

Next, we show that under a mild additional assumption on $w$, the Korenblum
principle fails for $A^{p}_{w}$ whenever $0<p<1$. It turns out that the same
example from [7] works, with appropriate modifications of the argument
employed there (and with some minor simplifications made in passing).

######  Theorem 3.

Let $0<p<1$ and let $w$ be a radial weight such that $\,\liminf_{r\to
0^{+}}w(r)>0$. Then the Korenblum domination principle fails for the space
$A^{p}_{w}$.

###### Proof.

Let $c\in(0,1)$ be arbitrary and choose $n$ large enough so that $n(1-p)>2$.
Consider the functions

$f(z)=\frac{c^{n}}{c^{n}+\varepsilon^{n}}\left(z^{n}+\varepsilon^{n}\right)\,,\quad
g(z)=z^{n}\,,$

where an appropriate small value of $\varepsilon\in(0,c)$ is to be chosen
later. Note that if $c\leq|z|<1$ then clearly

$|f(z)|\leq\frac{c^{n}}{c^{n}+\varepsilon^{n}}\left(|z|^{n}+\varepsilon^{n}\right)\leq\frac{c^{n}|z|^{n}+|z|^{n}\varepsilon^{n}}{c^{n}+\varepsilon^{n}}=|g(z)|\,.$

We need to show that $\|g\|_{p,w}<\|f\|_{p,w}$; that is,

(4)
$\frac{\left(c^{n}+\varepsilon^{n}\right)^{p}}{c^{np}}\int_{0}^{1}2rw(r)r^{np}\,dr<\int_{0}^{1}2rw(r)M_{p}^{p}(r;z^{n}+\varepsilon^{n})\,dr\,.$

In order to prove this, define

$I=\int_{0}^{\varepsilon}2rw(r)M_{p}^{p}(r;z^{n}+\varepsilon^{n})\,dr\,,\quad
J=\int_{\varepsilon}^{1}2rw(r)M_{p}^{p}(r;z^{n}+\varepsilon^{n})\,dr\,.$

By subharmonicity of the function $u(z)=|z^{n}+\varepsilon^{n}|^{p}$ and the
sub-mean value inequality, we have
$M_{p}^{p}(r;z^{n}+\varepsilon^{n})\geq\varepsilon^{np}$ for all $r\in(0,1)$,
hence

(5) $I\geq\varepsilon^{np}\int_{0}^{\varepsilon}2rw(r)\,dr=I^{\prime}\,.$

As was observed in [7], the inequality is actually strict; this can be seen
quickly by resorting to the Hardy-Stein identity. However, this is not really
needed for our purpose because in a chain of inequalities to be obtained later
there will be another strict inequality.

One can also see, as in [7, p. 501], that
$M_{p}^{p}(r;z^{n}+\varepsilon^{n})\geq r^{np}$, for $\varepsilon<r<1$. Again,
strict inequality can be proved but is not needed. Actually, this inequality
follows from the one observed earlier:
$M_{p}^{p}(r;z^{n}+\varepsilon^{n})\geq\varepsilon^{np}$. This is easily seen
as follows (by an obvious change of variable $\theta\mapsto-\theta$):

$M_{p}^{p}(r;z^{n}+\varepsilon^{n})=\int_{-\pi}^{\pi}|r^{n}e^{in\theta}+\varepsilon^{n}|^{p}\frac{d\theta}{2\pi}=\int_{-\pi}^{\pi}|r^{n}+\varepsilon^{n}e^{-in\theta}|^{p}\frac{d\theta}{2\pi}=\int_{-\pi}^{\pi}|r^{n}+\varepsilon^{n}e^{in\theta}|^{p}\frac{d\theta}{2\pi}\geq
r^{np}$

as before. Thus,

(6) $J\geq\int_{\varepsilon}^{1}2rw(r)r^{np}\,dr=J^{\prime}\,,$

In view of (5) and (6), by Bernoulli’s elementary inequality:

$\frac{\left(c^{n}+\varepsilon^{n}\right)^{p}}{c^{np}}=\left(1+\left(\frac{\varepsilon}{c}\right)^{n}\right)^{p}\leq
1+p\left(\frac{\varepsilon}{c}\right)^{n}\,,\quad 0<p<1\,,$

we see that, in order to obtain inequality (4), it suffices to show that

(7)
$I^{\prime}+J^{\prime}>\left(1+p\left(\frac{\varepsilon}{c}\right)^{n}\right)\int_{0}^{1}2rw(r)r^{np}\,dr\,.$

The last inequality is equivalent to

$\varepsilon^{np}\int_{0}^{\varepsilon}2rw(r)\,dr+\int_{\varepsilon}^{1}2rw(r)r^{np}\,dr>\left(1+p\left(\frac{\varepsilon}{c}\right)^{n}\right)\int_{0}^{1}2rw(r)r^{np}\,dr\,,$

which is the same as

(8)
$\varepsilon^{np}\int_{0}^{\varepsilon}2rw(r)\,dr-\int_{0}^{\varepsilon}2rw(r)r^{np}\,dr>p\left(\frac{\varepsilon}{c}\right)^{n}\int_{0}^{1}2rw(r)r^{np}\,dr\,.$

By putting together the two integrals on the left-hand side of (8) and
employing the obvious change of variable $r=\varepsilon\rho$, we see that (8)
amounts to

$\varepsilon^{np+2}\int_{0}^{1}(1-\rho^{np})2\rho
w(\varepsilon\rho)\,d\rho>p\left(\frac{\varepsilon}{c}\right)^{n}\int_{0}^{1}2rw(r)r^{np}\,dr\,.$

Changing back $\rho$ to $r$, this shows that the desired inequality (7) can be
rewritten as

(9) $\int_{0}^{1}(1-r^{np})2rw(\varepsilon
r)\,dr>\frac{p}{c^{n}}\varepsilon^{n(1-p)-2}\int_{0}^{1}2rw(r)r^{np}\,dr\,.$

In view of our initial assumption on $n$, the right-hand side of (9) tends to
zero as $\varepsilon\to 0^{+}$. Without using any assumption on the continuity
of $w$, one can show by elementary calculus that

$\liminf_{\varepsilon\to 0^{+}}w(\varepsilon r)=\liminf_{s\to 0^{+}}w(s)$

for every fixed $r\in(0,1)$. (Just consider the appropriate sequences to see
that each quantity is greater than or equal to the other.) Thus, taking the
$\liminf_{\varepsilon\to 0^{+}}$ of the left-hand side of (9), by Fatou’s
lemma we obtain that

$\liminf_{\varepsilon\to 0^{+}}\int_{0}^{1}(1-r^{np})2rw(\varepsilon
r)\,dr\geq\liminf_{s\to
0^{+}}w(s)\int_{0}^{1}(1-r^{np})2r\,dr=\frac{np}{np+2}\liminf_{s\to
0^{+}}w(s)>0\,.$

This shows that we can choose a small enough positive $\varepsilon$ so that
(9) and therefore (7) holds, which proves (4). ∎

Finally, it should be noted that, in order for the Korenblum principle to fail
for $0<p<1$, an assumption like $\,\liminf_{r\to 0^{+}}w(r)>0$ on the weight
cannot be dispensed with, as the following simple example shows.

######  Example 4.

Let

$w(r)=\left\\{\begin{array}[]{ll}0&\mbox{if \ $0\leq r<R$},\\\ 1&\mbox{if \
$R\leq r<1$}\,.\end{array}\right.$

Then $w$ clearly does not satisfy the condition $\,\liminf_{r\to
0^{+}}w(r)>0$, yet the Korenblum principle trivially holds for $A^{p}_{w}$ for
every $p>0$ with $c=R$ since

$\|f\|_{p,w}=\left(\int_{R}^{1}2rw(r)M_{p}^{p}(r;f)\,dr\right)^{1/p}\,.$

###  2.7. Conclusions

It should be pointed out that the main credit for our findings belongs to the
previous authors: to Hayman [4] and Hinkkanen [6] for devising the method used
in the proof of Theorem 2 and to Schuster [9] for making clear the estimates
used for all values of $p$, as well as to Božin and Karapetrović’s [1] for
producing a single and simple counterexample that works for all $p\in(0,1)$,
even in the general context of our Theorem 3. Also, Karapetrović very clear
exposition in [7] motivated this work.

Our contribution consists primarily in realizing that the estimates obtained
by the other authors actually work in the context of general Bergman spaces
with radial weights in place of the spaces with standard radial weights, and
also in noticing that all the estimates work without having to estimate any
values of the Beta function that come up in the computation for the standard
cases.

Our findings seem to show that the Korenblum principle is a purely classical
analysis statement that only requires standard estimates on integral means,
without resorting to any abstract context or “soft analysis” methods, and its
validity or failure are independent of the question of completeness of the
weighted Bergman space considered, which is another remarkable fact.

## References

  * [1] V. Božin, B. Karapetrović, Failure of Korenblum’s maximum principle in Bergman spaces with small exponents, Proc. Amer. Math. Soc. 146 (2018), 2577–-2584.
  * [2] P.L. Duren, Theory of $H^{p}$ Spaces, Academic Press, New York 1970.
  * [3] P. L. Duren and A. P. Schuster, _Bergman Spaces_ , Math. Surveys and Monographs Vol. 100, American Mathematical Society, Providence, Rhode Island 2004.
  * [4] W.K. Hayman, On a conjecture of Korenblum, Analysis 19 (1999), 195–205.
  * [5] H. Hedenmalm, B. Korenblum, K. Zhu, _Theory of Bergman Spaces_ (Graduate Texts in Mathematics, Vol. 199, Springer, New York, Berlin, etc. 2000.
  * [6] A. Hinkkanen, On a maximum principle for Bergman spaces, J. Anal. Math. 79 (1999), 335–344.
  * [7] B. Karapetrović, Korenblum maximum principle in mixed norm spaces, Arch. Math. 118 (2022), 497–-507.
  * [8] B. Korenblum, A maximum principle for the Bergman space, Publ. Mat. 35 (1991), 479–486.
  * [9] A. Schuster, The maximum principle for the Bergman space and the Möbius pseudodistance for the annulus, Proc. Amer. Math. Soc. 134 (2006), 3525–-3530.
  * [10] C. Wang, Refining the constant in a maximum principle for the Bergman space, Proc. Amer. Math. Soc. 132 (2003), 853–855.
  * [11] C. Wang, On Korenblum’s maximum principle, Proc. Amer. Math. Soc. 134 (2006), 2061–-2066.
  * [12] C. Wang, Domination in the Bergman space and Korenblum’s constant. Integral Equations Operator Theory 61 (2008), 423–-432.
  * [13] C. Wang, Some results on Korenblum’s maximum principle, J. Math. Anal. Appl. 373 (2011), 393–-398.