#  Padé: A code for protoplanetary disk turbulence based on Padé differencing

Karim Shariff NASA Ames Research Center
Moffett Field, CA 94035, USA

###### Abstract

The Padé code has been developed to treat hydrodynamic turbulence in
protoplanetary disks. It solves the compressible equations of motion in
cylindrical coordinates. Derivatives are computed using non-diffusive and
conservative fourth-order Padé differencing, which has higher resolving power
compared to both dissipative shock-capturing schemes used in most astrophysics
codes, as well as non-diffusive central finite-difference schemes of the same
order. The fourth-order Runge-Kutta method is used for time stepping. A
previously reported error-corrected Fargo approach is used to reduce the time
step constraint imposed by rapid Keplerian advection. Artificial bulk
viscosity is used when shock-capturing is required. Tests for correctness and
scaling with respect to the number of processors are presented. Finally,
efforts to improve efficiency and accuracy are suggested.

protoplanetary disks, hydrodynamics, turbulence, methods:numerical

††journal: ApJ Supp.

##  1 Introduction

This paper presents a code called Padé, developed to simulate hydrodynamic
turbulence in protoplanetary disks. The code uses a Padé/compact finite-
difference scheme (Lele, 1992). Such schemes have spectral-like resolving
power. This means that they approach the ability of a spectral method to
compute derivatives exactly for all wave numbers that the mesh can support. In
particular, for advection problems they have zero diffusive error (as do all
central finite-difference schemes) and have small dispersion error across the
wavenumber range. Thus they have the ability to more accurately treat the
dynamics of the small scales supported by the mesh and produce energy spectra
that have a wider inertial (power-law) range for turbulent flows.

On the other hand, most astrophysical codes employ Godunov-type schemes that
were an elegant and mathematically supported breakthrough for capturing shock-
waves. Such schemes employ non-central upwinded finite-difference or flux
reconstruction schemes with a smart non-linear dissipation that is necessary
for capturing shocks, but which leads to excessive dissipation of vortical and
other smooth features such as density waves. To fix this issue there have been
attempts along many directions which include hybrid methods (Adams & Shariff,
1996; Pirozzoli, 2002), non-linear filtering (Yee & Sjögreen, 2018), and
vorticity-preserving schemes (Lerat et al., 2007; Seligman & Laughlin, 2017;
Seligman & Shariff, 2019). Since shocks have not been observed in
protoplanetary disk turbulence to date, we can side step the issue. For
treating shock waves that are not too strong, the Padé code provides an
optional artificial bulk viscosity treatment (Cook & Cabot, 2005; Mani et al.,
2009). This is not as elegant or oscillation free as Godunov methods, but is
designed to apply a diffusivity only where the divergence
$\nabla\cdot\mathbf{u}$ is strong.

The development of the Padé code was motivated by the discovery in the last
two decades of a number of mechanisms for the generation of hydrodynamic
turbulence in protoplanetary disks; see Lesur et al. (2022) for a
comprehensive review. These include the vertical shear instability (VSI),
convective over-stability (COS), and the zombie vortex instability (ZVI). Each
is most strongly amplified for a different range of
$\Omega\,t_{\mathrm{thermal}}$, the time scale for radiative relaxation of
temperature fluctuations back to the background (normalized by the orbital
frequency, $\Omega$). VSI is most strongly amplified when
$\Omega\,t_{\mathrm{thermal}}\ll 1$, COS when
$\Omega\,t_{\mathrm{thermal}}\sim 1$, and ZVI when
$\Omega\,t_{\mathrm{thermal}}\gg 1$. Turbulence can also be driven by the
magneto-rotational instability (MRI) in the radially inner and outer regions
of protoplanetary disks where ionization is sufficient.

As mentioned earlier, most astrophysical codes that are applied to
protoplanetary disks use dissipative shock-capturing methods. These codes
include Athena (Stone et al., 2008), Athena++ (Stone et al., 2020), Pluto
(Mignone et al., 2007), and Fargo3D (Benítez-Llambay & Masset, 2016). The
highest order scheme provided in Athena and Athena++ is the third order PPM
(Piecewise Parabolic Method). The highest order scheme in Pluto is a fifth
order WENO (Weighted Essentially Non-Oscillatory) finite difference shock-
capturing scheme denoted WENOZ_FD in the Pluto manual (Mignone, 2021, pg. 93).
Fargo3d employs a staggered mesh such that the scalar variables, density and
internal energy per unit volume, are cell-centered while vector quantities,
velocity and magnetic field, are located at the centers of cell faces. The
computation of fluxes at cell faces employs upwinded interpolation which is
necessarily dissipative. Time advancement uses operator splitting whereby
different sets of terms contributing to the time rate of change of flow
quantities are time-advanced separately, one after the other.

The one exception to the use of dissipative shock-capturing scheme is the
Pencil code (Brandenburg et al., 2021) which uses sixth-order central
differencing. Central schemes, including the Padé scheme, produce oscillations
at the Nyquist wavenumber of the mesh. To overcome this, Pencil uses a
dissipative fifth order upwind biased scheme, while Padé uses a filter with a
sharp cut-off (§3.3).

The rest of the paper presents the equations solved (§2), computational
schemes (§3), various tests (§4), and closing remarks (§5). Code availability
is described after the acknowledgements.

##  2 Equations discretized

###  2.1 Transport equations

The code solves the equations for mass and momentum (radial, angular, and
vertical) transport written in as close to flux-divergence form as possible in
cylindrical coordinates

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial z}(\rho
u_{z})+\frac{1}{r}\frac{\partial}{\partial r}(r\rho
u_{r})+\frac{1}{r}\frac{\partial}{\partial\phi}(\rho u_{\phi})=0,$ (1)
$\displaystyle\frac{\partial}{\partial t}(\rho u_{r})+\frac{\partial}{\partial
z}(\rho u_{r}u_{z})+\frac{1}{r}\frac{\partial}{\partial r}\left[r(p+\rho
u_{r}^{2})\right]+\frac{1}{r}\frac{\partial}{\partial\phi}(\rho
u_{\phi}u_{r})-\frac{\rho u_{\phi}^{2}}{r}=\frac{p}{r}+\rho
g_{r}+F^{\mathrm{v}}_{r},$ (2) $\displaystyle\frac{\partial}{\partial t}(\rho
u_{\phi}r)+\frac{\partial}{\partial z}(\rho
u_{\phi}ru_{z})+\frac{1}{r}\frac{\partial}{\partial r}(r\rho
u_{\phi}ru_{r})+\frac{\partial}{\partial\phi}(p+\rho
u_{\phi}^{2})=rF^{\mathrm{v}}_{\phi},$ (3)
$\displaystyle\frac{\partial}{\partial r}(\rho u_{z})+\frac{\partial}{\partial
z}(p+\rho u_{z}^{2})+\frac{1}{r}\frac{\partial}{\partial r}(r\rho
u_{z}u_{r})+\frac{1}{r}\frac{\partial}{\partial\phi}(\rho u_{z}u_{\phi})=\rho
g_{z}+F^{\mathrm{v}}_{z}.$ (4)

Here $\mathbf{F}^{\mathrm{v}}$ represents an optional viscous force given in
Appendix A. Implementation of characteristic boundary conditions requires
calculation of the flux Jacobian, and to make the flux look similar to that in
Cartesian coordinates, the term $-\partial p/\partial r$ that normally appears
on the right-hand-side of the radial momentum equation (2) was moved into the
radial advective flux on the left-hand-side. This is accomplished by adding
the quantity

$\frac{1}{r}\frac{\partial}{\partial
r}\left(rp\right)=\frac{1}{r}\left(r\frac{\partial p}{\partial
r}+p\right)=\frac{\partial p}{\partial r}+\frac{p}{r},$ (5)

to both sides of (2), which results in the source term $p/r$ on the right-
hand-side.

For the locally isothermal option, the pressure is computed as $p=\rho
c_{\mathrm{i}}^{2}(r)$ where $c_{\mathrm{i}}(r,z)$ is the local isothermal
sound speed. For the non-isothermal option we have

$p=(\gamma-1)e_{\mathrm{int}},$ (6)

where $e_{\mathrm{int}}$ is the internal energy (per unit volume). The
rationale for using the internal energy instead of its sum with the kinetic
energy is given in our work on the Fargo method (Shariff & Wray, 2018). The
transport equation for $e_{\mathrm{int}}$ is

$\frac{\partial e_{\mathrm{int}}}{\partial t}+\frac{\partial}{\partial
z}(e_{\mathrm{int}}u_{z})+\frac{1}{r}\frac{\partial}{\partial
r}(re_{\mathrm{int}}u_{r})+\frac{1}{r}\frac{\partial}{\partial\phi}(e_{\mathrm{int}}u_{\phi})=-p\nabla\cdot\vec{u}+Q^{\mathrm{v}}-\nabla\cdot\vec{q}_{\mathrm{cond}}.$
(7)

The first term on the right-hand-side of (7) is the pressure-dilatation term
which causes heating under compression. The dilatation is given by

$\nabla\cdot\vec{u}=\frac{\partial u_{z}}{\partial
z}+\frac{1}{r}\frac{\partial}{\partial r}(ru_{r})+\frac{1}{r}\frac{\partial
u_{\phi}}{\partial\phi}$ (8)

The last two terms $Q^{\mathrm{v}}-\nabla\cdot\vec{q}_{\mathrm{cond}}$ on the
right-hand-side of (7) represent viscous heating and conductive heat transfer,
respectively. They are activated only if viscous terms are activated; their
form is given in Appendix A.

The variables evolved at each grid point are $\vec{q}=(\rho,\rho u_{r},\rho
u_{\phi}r,\rho u_{z},e_{\mathrm{int}})$. When spatial derivatives are
approximated by Padé finite-differences (discussed below), one obtains a
system of ODEs for the time rate of change, $\partial_{t}\vec{q}$ of $\vec{q}$
at each grid point. This system is evolved using the fourth-order Runge-Kutta
method. Note that we do not employ operator splitting whereby different sets
of terms in $\partial_{t}\vec{q}$ are time evolved separately one after
another; this would produce a splitting error which we do not incur.

Currently, the main source of gravity in the code is from a mass $M$ at the
origin so that

$g_{r}=-\frac{GM}{R}\frac{r}{R},\hskip
7.11317ptg_{z}=-\frac{GM}{R}\frac{z}{R},\hskip
7.11317ptR\equiv\left(r^{2}+z^{2}\right)^{1/2},$ (9)

which are pre-computed. The code also allows other simple choices such as
uniform $g_{z}$ and the thin disk version of (9).

###  2.2 Fargo

The Fargo method was introduced by Masset (2000) to alleviate the time step
restriction resulting from fast Keplerian advection. Shariff & Wray (2018)
improved its accuracy by first noting that underlying the method is a
transformation of the azimuthal coordinate $\phi$

$\displaystyle\phi^{\prime}$
$\displaystyle=\phi-\overline{\Omega}(r)(t-t_{0}),$ (10) $\displaystyle
r^{\prime}$ $\displaystyle=r,$ (11) $\displaystyle t^{\prime}$
$\displaystyle=t,$ (12)

for the duration of a time step, $t_{0}<t\leq t_{0}+\Delta t$. Here
$\overline{\Omega}(r)$ is a prescribed rotation rate that one wishes to
subtract from determining the time step. At the end of the time step, one
brings the flow field back to original coordinates by performing a shift using
an FFT. The chain-rule for differentiation then implies that every $t$ and $r$
derivative in the transport equations carries an additional term:

$\displaystyle\frac{\partial}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
t^{\prime}}-\frac{\partial\phi^{\prime}}{\partial
t}\frac{\partial}{\partial\phi^{\prime}}=\frac{\partial}{\partial
t^{\prime}}-\overline{\Omega}(r)\frac{\partial}{\partial\phi^{\prime}},$ (13)
$\displaystyle\frac{\partial}{\partial r}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
r^{\prime}}-\frac{\partial\phi^{\prime}}{\partial
r}\frac{\partial}{\partial\phi^{\prime}}=\frac{\partial}{\partial
r^{\prime}}-(t-t_{0})\frac{\partial\overline{\Omega}}{\partial
r}\frac{\partial}{\partial\phi^{\prime}},$ (14)
$\displaystyle\frac{\partial}{\partial\phi}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial\phi^{\prime}}.$ (15)

The second term in $\partial/\partial t$ serves to remove
$\overline{\Omega}(r)$ from the azimuthal advection velocity, therefore, it no
longer influences the time step. To see how this works, consider the
transformed mass transport equation (1):

$\frac{\partial\rho}{\partial
t^{\prime}}+\frac{1}{r}\frac{\partial}{\partial\phi}\left(\rho
u^{\prime}_{\phi}\right)+\frac{1}{r}\left[\frac{\partial}{\partial
r^{\prime}}+\chi\right]\left(r\rho u_{r}\right)+\frac{\partial}{\partial
z}\left(\rho u_{z}\right)=0,$ (16)

where

$u^{\prime}_{\phi}=u_{\phi}-\overline{\Omega}(r)r$ (17)

is a shifted velocity that results in a less restrictive time step. In (16),
the symbol $\chi$ denotes the operator

$\chi=-(t-t_{0})\frac{\partial\overline{\Omega}}{\partial
r}\frac{\partial}{\partial\phi}.$ (18)

It arises from the second term on the right-hand-side of the last member of
(14), and must be added to every $r$ derivative in the transport equations;
the code refers to it as the ‘extra operator.’ The extra operator $\chi$ is
missing in the original implementation (Masset, 2000): its neglect results in
an error of $\mathcal{O}\left(\Delta t\right)$. The original implementation
also “integerizes” the prescribed $\overline{\Omega}(r)$ to allow a shift of
the flow field by an integer number of grid intervals in $\phi$ at the end of
a time step. The integer shift jumps at certain radial locations and this
results in additional error at those locations.

This code provides options to completely or partially revert to Masset’s
original algorithm if desired. In particular, integer shifting can be selected
as opposed to the more accurate real-valued shifting. Also, inclusion of the
extra operator $\chi$ can be suppressed. At the end of each time step, the
flow field is shifted back to original coordinates. Real-valued shifting is
performed using an FFT from the FFTW library (Frigo & Johnson, 2005).

Following publication of Shariff & Wray (2018), we were contacted by Prof.
Pablo Benítez Llambay (Universidad Adolfo Ibáñez, Chile) who is a main
developer of Fargo3d. He pointed out that Fargo3d performs advection by
operator splitting, and therefore there is no coordinate transformation
requiring chain-rule terms. For instance, consider the continuity equation
with only $\phi$ advection for simplicity:

$\frac{\partial\rho}{\partial
t}+\frac{1}{r}\frac{\partial}{\partial\phi}\left(\rho u_{\phi}\right)=0,$ (19)

with the decomposition

$u_{\phi}=\overline{\Omega}(r)r+u_{\phi}^{\prime}.$ (20)

In the operator splitting approach, $\rho$ is first partially evolved
according to advection by the residual velocity $u_{\phi}^{\prime}$,

$\frac{1}{r}\frac{\partial}{\partial\phi}\left(\rho u_{\phi}^{\prime}\right),$
(21)

whose time step restriction is not severe. Next, the term
$\overline{\Omega}(r)\partial_{\phi}\rho$, representing advection by Keplerian
flow, is treated in Fargo3d by real-valued shifting (Colin McNally, Queen Mary
Univ. of London, Private Communication). Thus, it is true that no coordinate
transformation is implied. However, operator splitting has
$\mathcal{O}\left(\Delta t^{2}\right)$ error which is of the same order as the
error made by dropping the extra chain-rule terms. Thus, it would appear that
for higher order time integration schemes, chain-rule terms are necessary for
consistency.

###  2.3 Artificial pressure for shock-capturing

Padé schemes, like spectral schemes, are not designed to capture shocks.
However, one may encounter shocks in protoplanetary disks. Examples include
the infall accretion shock (Neufeld & Hollenbach, 1994), and bow shocks due to
solid bodies moving supersonically relative to the gas. For these reasons, the
code allows for an optional treatment of shocks using artificial bulk
viscosity (Cook & Cabot, 2005; Mani et al., 2009). This method results in an
artificial pressure, $p_{\mathrm{art}}$, which is then added to the physical
pressure. The actual calculation of $p_{\mathrm{art}}$ is

$p_{\mathrm{art}}=-\beta_{\Delta}\nabla\cdot\vec{u},$ (22)

where $\beta_{\Delta}$ is the artificial bulk viscosity. Note that artificial
pressure is positive in regions of compression (dilatation
$\nabla\cdot\vec{u}<0$) and negative in regions of expansion. The artificial
bulk viscosity is made sensitive to the dilatation as follows:

$\beta_{\Delta}=C_{\mathrm{ap}}\rho\ell_{\mathrm{grid}}^{2}|\nabla\cdot\vec{u}|,$
(23)

were $C_{\mathrm{ap}}$ is a user specified coefficient and the grid size
squared is

$\ell_{\mathrm{grid}}^{2}=\begin{cases}\Delta z^{2},&\text{\ in 1D};\\\
r\Delta\phi\Delta r,&\text{\ in the planar case};\\\ (r\Delta\phi\Delta
r\Delta z)^{2/3},&\text{\ in 3D}.\end{cases}$ (24)

Since the absolute value function is not smooth, the Padé filter (§3.3) with
$\epsilon_{\mathrm{filter}}=0.2$ is applied to $\beta_{\Delta}$.

The artificial pressure term imposes a time step constraint appropriate for a
viscous (second derivative) term. The time step must satisfy

$(\lambda_{\mathrm{ap}})_{\mathrm{max}}\Delta t<\mathrm{CFL},$ (25)

where $(\lambda_{\mathrm{ap}})_{\mathrm{max}}$ is the maximum eigenvalue of
the bulk viscosity operator and $\mathrm{CFL}$ is the Courant-Friedrichs-Lewy
limit specified by the user. For the fourth-order Runge-Kutta method,
$\mathrm{CFL}<\sqrt{8}$ for stability. The maximum eigenvalue is estimated as

$(\lambda_{\mathrm{ap}})_{\mathrm{max}}=\left(\frac{\pi^{2}\beta_{\Delta}}{\rho\ell_{\mathrm{grid}}^{2}}\right)_{\mathrm{max}},$
(26)

where a factor of $\pi$ comes from assuming that each numerical derivative has
spectral accuracy. Recall that $\beta_{\Delta}\propto\ell_{\mathrm{grid}}^{2}$
so the grid size actually drops out in (26). For the shock-tube and density
wave test cases we report below, this eigenvalue was a factor of 1.08 and 1.61
larger, respectively, than the eigenvalue for the Euler terms.

##  3 Numerics

###  3.1 Padé differentiation

The motivation for Padé differencing (also referred to as compact
differencing) is its high resolving power (§3.2); see Lele (1992) for a
comprehensive presentation and various extensions. The idea and initial
development of such schemes is due to the Czech born astronomer and numerical
analyst Zdeněk Kopal around 1959; see the bibliographical notes (item IX-C) in
his book (Kopal, 1961). The basic idea, developed in Chapter 9 of the book, is
to write the exact first derivative operator as an exact function of the
central difference operator. This operator function is first expanded in a
truncated Taylor polynomial, which results in a conventional difference
scheme. The key idea, however, is to next obtain a rational polynomial
approximation (known as a Padé approximant) to the Taylor series. It is known
that Padé111Henri Padé was a French mathematician who, for his doctoral
thesis, studied (c. 1890) the approximation of functions by rational
polynomials, now known as Padé approximants. approximants to ordinary
functions have a greater range of accuracy and radius of convergence than a
Taylor series. In the present case, one obtains difference schemes with better
resolving power. The simplest scheme, and the one we use both in the interior
of non-periodic domains and for periodic domains is

$\alpha f^{\prime}_{j-1}+f^{\prime}_{j}+\alpha
f^{\prime}_{j+1}=\frac{a}{2h}\left(f_{j+1}-f_{j-1}\right),$ (27)

where $h$ is the uniform grid spacing, $\alpha=1/4$ and $a=3/2$. Equations
(27) constitute a tridiagonal system of equations along each line of data in
the mesh. The system is solved efficiently using the Thomas algorithm, which
is simply Gaussian elimination applied to a tridiagonal matrix. Since Gaussian
elimination is recursive, namely, operating on row $j$ depends on the result
of row $j-1$, the memory cache cannot be preloaded with the required data.
Similarly, any available vectorization hardware cannot be engaged by the
compiler. To overcome this, we follow the standard practice of having each
step of the Thomas algorithm inner-loop over a bundle of independent
inversions for different grid lines of the mesh. Equation (27) is fourth-order
accurate, i.e., its truncation error is $\mathcal{O}\left(h^{4}\right)$. To
allow for non-uniform meshes, numerical differentiation is performed with
respect to the continuous grid index variable $\xi$ (such that $\xi_{j}=j$ and
$h=1$), and the chain rule is used, e.g.,

$\frac{\partial f}{\partial z}=\frac{\partial
f}{\partial\xi}\left(\frac{dz}{d\xi}\right)^{-1}.$ (28)

Kopal’s operator calculus is abstruse and unwieldy, but once the basic form of
Padé schemes is recognized, an easier approach to develop them is to simply
write down a specific form with a desired grid-point stencil, and obtain the
unknown coefficients by setting the Taylor series error to zero at various
orders. This is the approach followed in Lele (1992) who developed a number of
extensions, including resolving power optimization, boundary schemes, higher
derivatives, conservation, and filtering.

For robustness, it is desirable that a code maintain positivity of density and
internal energy. Negative values can arise in strongly evacuated regions for
low-order schemes, and near very strong discontinuities for higher order
schemes; see, for example, Hu et al. (2013) who present a simple method for
ensuring positivity for finite-difference schemes that can be written as a
difference of numerical fluxes. The present method does not guarantee
positivity, however, the code has not encountered difficulties for problems of
subsonic turbulence for which it is intended. In the future, it may be
possible to implement the method of Hu et al. (2013) using the “reconstruction
by primitive function” trick (Harten et al., 1987; Shu & Osher, 1989), which
is also discussed in Merriman (2003). In this method, one obtains numerical
fluxes by differentiating the primitive function; this differentiation would
be performed using the Padé scheme.

###  3.2 Resolving power of Padé differencing

The advantage of Padé schemes is, first of all, their compactness: for (27),
for instance, fourth-order accuracy is obtained with a stencil width of three
rather than four in the case of standard central differencing. More
importantly, they have better resolving power than standard central
differencing. This means that for the same formal order of accuracy, they
provide an accurate derivative up to higher wave numbers. This is quantified
by the so-called effective wave number analysis: substitute

$f_{j}=\exp(\mathrm{i}kx_{j})\hskip 7.11317pt\mathrm{and}\hskip
7.11317ptf^{\prime}_{j}=\mathrm{i}k_{\mathrm{eff}}(k)\exp(\mathrm{i}kx_{j}),\hskip
7.11317pt0\leq k\leq\pi/h,$ (29)

into (27) to obtain the effective wavenumber $k_{\mathrm{eff}}(k)$.

Figure 1: (a) Effective wavenumber for various difference schemes. (b)
Transfer function, $T(k)$, for the Padé filter. The grid spacing is $h$. The
limit $kh=\pi$ corresponds to the Nyquist mode which has wavelength $2h$.

Figure 1a displays the effective wavenumber for various schemes and shows that
the fourth-order Padé scheme has better resolving power than a conventional
sixth order scheme. In general, $k_{\mathrm{eff}}(k)$ is complex and its
imaginary part represents numerical dissipation of the scheme. The fact that
$k_{\mathrm{eff}}(k)$ is real for central schemes means that they have
dispersion error but no dissipation error (for periodic boundary conditions).
As as example of how to use the effective wavenumber diagram, we estimate by
eye that the highest wavenumber for which we can trust the 4th order Padé
scheme is $kh=1.5$ which implies that the smallest number of grid points per
wavelength for which the scheme is accurate is about four.

###  3.3 Padé filtering

It is known that in the presence of non-linearity or non-uniform grids, Padé
schemes produce small $2\Delta$ waves at every time step which can grow if not
controlled; this is true for conventional central schemes as well. Our remedy
is to apply a minimum amount of Padé filtering (Lele, 1992, §C.2) which has a
sharp cut-off and targets the very highest wave numbers. We also use Padé
filtering as an implicit sub-grid treatment; this is discussed in §4 where we
use it for this purpose in an axisymmetric simulation of vertical shear
instability (VSI).

We use the fourth-order filter of Lele (1992) which for periodic boundary
conditions, or in the interior of non-periodic domains, has the form

$aU_{j-1}+U_{j}+aU_{j+1}=P(u_{j-2}+u_{j+2})+Q(u_{j-1}+u_{j+1})+Ru_{j},$ (30)

where $U_{j}$ are the filtered $u_{j}$. The conditions for fourth-order
accuracy are:

$a=-\frac{1}{2}+2Q,\hskip 14.22636ptR=\frac{1}{2}\left(2+3a-3Q\right),\hskip
14.22636ptP=\frac{1}{4}\left(a-Q\right).$ (31)

The reader can verify that there is no filtering when $Q=1/2$. To specify the
strength of the filter, the code uses the parameter
$\epsilon_{\mathrm{filter}}$ such that

$Q=\frac{1}{2}-\frac{1}{4}\epsilon_{\mathrm{filter}}.$ (32)

The transfer function $T(k)$ of the filter versus wavenumber can be obtained
by substituting $u_{j}=e^{\mathrm{i}kx_{j}}$ and
$U_{j}=T(k)e^{\mathrm{i}kx_{j}}$ (with $x_{j}=jh$) into (30). Figure 1b shows
$T(k)$ for various values of $\epsilon_{\mathrm{filter}}$. Axisymmetric
simulations for $\epsilon_{\mathrm{filter}}$ values ranging from $0.01$ to
$0.125$ will be presented in §4.5.

For non-periodic boundaries, we use the boundary treatment developed by Alan
Wray (private communication) that is conservative, i.e., it preserves

$\sum_{j=1}^{N}u_{j}h_{j}.$ (33)

The formulae applied at $j=2$ and $N-1$ are

$\displaystyle BU_{1}+CU_{2}+DU_{3}$
$\displaystyle=Eu_{1}+Fu_{2}+Gu_{3}+Hu_{4},$ (34) $\displaystyle
BU_{N}+CU_{N-1}+DU_{N-2}$ $\displaystyle=Eu_{N}+Fu_{N-1}+Gu_{N-2}+Hu_{N-3}.$
(35)

with

$B=2a,\hskip 8.5359ptC=1+a,\hskip 8.5359ptD=a,\hskip
8.5359ptE=\frac{1}{4}(7a+Q),\hskip 8.5359ptF=\frac{1}{4}(4+7a-3Q),\hskip
8.5359ptG=\frac{1}{4}(a+3Q),\hskip 8.5359ptH=\frac{1}{4}(a-Q)$ (36)

The boundary values $u_{1}$ and $u_{N}$ are unchanged.

The first item of Table IX in Lele (1992) lists the the leading order
truncation error for the filter (30):

$U(x)=\left[1+\frac{1}{16}\left(1-2a\right)h^{4}\partial^{4}/\partial
x^{4}+\cdots\right]u(x).$ (37)

Next consider the expression

$u(x,t+\Delta t)=(1+\nu_{4}\Delta t\partial^{4}/\partial x^{4})u(x),$ (38)

which represents application of an Euler step to

$\frac{\partial u}{\partial t}=\nu_{4}\frac{\partial^{4}u}{\partial x^{4}}$
(39)

Comparing (37) and (38) gives

$\nu_{4}=\frac{1}{16}(1-2a)\frac{h^{4}}{\Delta t}.$ (40)

Hence, to leading order, Padé filtering corresponds to fourth-order
hyperviscosity. Substituting the first member of (31) and (32) into (40) gives

$\nu_{4}=\frac{1}{16}\frac{\epsilon_{\mathrm{filter}}}{\Delta t}h^{4}.$ (41)

Equation (40) implies that as the time step is reduced
$\epsilon_{\mathrm{filter}}$ should be reduced. It should be noted that the
matrix associated with the Padé filter becomes ill-conditioned for very small
values of $\epsilon_{\mathrm{filter}}$ and a sufficiently large number of grid
points. To alleviate this, the filter is applied every $N_{\mathrm{filter}}$
time steps (rather than after every step) and with
$\epsilon_{\mathrm{filter}}$ increased by $N_{\mathrm{filter}}$. In that case
one should replace $\Delta t$ with $\Delta tN_{\mathrm{filter}}$ in (41).

###  3.4 Boundary schemes and global conservation

At the end points of a non-periodic direction ($j=1$ and $N$), the scheme (27)
involves values and derivatives outside the domain. Therefore, at $j=1$ and
$N$ we use the third-order one-sided scheme from equation (4.1.1) in Lele
(1992):

$\displaystyle f^{\prime}_{1}+\alpha_{1}f^{\prime}_{2}$
$\displaystyle=h^{-1}\left(a_{1}f_{1}+b_{1}f_{2}+c_{1}f_{3}\right),$ (42)
$\displaystyle f^{\prime}_{N}+\alpha_{1}f^{\prime}_{N-1}$
$\displaystyle=-h^{-1}\left(a_{1}f_{N}+b_{1}f_{N-1}+c_{1}f_{N-2}\right),$ (43)

with

$\alpha_{1}=2,\hskip 7.11317pta_{1}=-15/16,\hskip 7.11317ptb_{1}=2,\hskip
7.11317pt\mathrm{and}\hskip 7.11317ptc_{1}=1/2.$ (44)

According to a theorem for hyperbolic initial boundary value problems, one can
reduce the order of accuracy at the boundary by one without affecting the
global order of accuracy (Gustafsson, 1981). The scheme used at the interior
points $j\in[2,N-1]$ is (27).

We wish the finite-difference scheme to possess a discrete conservation
property. Specifically, we require that a discrete version of the Leibniz
integral rule be satisfied: a numerical integral of the numerical derivative
should reduce to a difference of boundary values. The is achieved for a Padé
scheme as described in Lele (1992, §4.2) and Brady & Livescu (2019). Its
application to the present scheme is described in Appendix B. Briefly, the
condition to be satisfied is that the row entries in columns 2 to $N-1$ of
matrix $B$ must have a weighted sum of zero. Here $B$ is the matrix
representation of the right-hand-side of (27):

$A\mathbf{f}^{\prime}=\frac{1}{h}B\mathbf{f}.$ (45)

Appendix B shows that this condition is satisfied for the present boundary
scheme. An alternate way to obtain conservation is the “reconstruction by
primitive function” trick referred to above.

###  3.5 Data partitioning for parallelization

Since an entire line of data along $x$ is needed to compute a Padé derivative
along $x$, for any direction $x$, we employ the pencil data structure. Each
processor is assigned a pencil of data that can be thought of as a long brick
with the long side along the direction of differentiation. Most of the work is
done with $z$-pencils (i.e., with $z$ as the long direction). This work
includes initialization, output, and, time integration sub-steps. To compute
$r$ derivatives, we perform a so-called “transpose” such that each processor
also has $r$ pencils, and similarly for $\phi$ derivatives. The flowfield
arrays with $z$, $r$, and $\phi$ pencilling, respectively, are dimensioned as
follows:

    q          (sr:er, sphi:ephi, nz, ndof)         ! pencil along z
    q_r_space  (sphi:ephi, sz_r:ez_r, ndof, nr)     ! pencil along r
    q_phi_space(sr:er, sz_phi:ez_phi, ndof, nphi)   ! pencil along phi

Here ndof refers to the number of degrees of freedom (number of flow
variables) at each grid-point and the other dimensions can be read, for
example, as follows: sz_r:ez_r—starting $z$ index to ending $z$ index for an
$r$ pencil. Partitioning and transpose routines were taken from Alan Wray’s
Stellarbox code.

##  4 Tests

###  4.1 Convergence for two-dimensional advection

Here we solve the equation for a scalar, $f(z,\phi,t)$, uniformly advecting in
the $z$ and $\phi$ directions:

$\frac{\partial f}{\partial t}+\frac{\partial f}{\partial z}+\frac{\partial
f}{\partial\phi}=0,\hskip 8.5359ptz\in[0,2\pi),\hskip
2.84544pt\phi\in[0,2\pi),$ (46)

with periodic boundary conditions and the initial condition

$f(z,\phi,0)=1+\frac{1}{2}\sin 2z\sin 2\phi.$ (47)

The CFL is fixed at unity and the number of grid points $n_{z}$ and $n_{\phi}$
is varied.

Figure 2: Convergence for two-dimensional advection.

Figure 2 plots the rms error at $t=6\pi$ compared with the exact solution at
$t=6\pi$. The rate of convergence is seen to be fourth-order.

###  4.2 One-dimensional Euler equations with shocks

Figure 3: 1D Euler equation tests. (a)–(c) : Shock tube. (a) Density. (b)
Pressure. (c) Velocity. (d) Mach 3 shock propagating into density waves (Shu &
Osher, 1989). The result of the present code is compared with the WENO5-Roe
method. For both tests the coefficient of artificial pressure
$C_{\mathrm{ap}}=2$ and Pade/compact filter strength
$\epsilon_{\mathrm{filter}}=0.05$ (applied at the end of every step).
$\mathrm{CFL}=1$.

Figure 3 presents two test cases for the one-dimensional Euler equations which
was run in the code’s $z$ direction by suppressing the other two. Both tests
employed $n_{z}=512$ grid points and artificial pressure (with
$C_{\mathrm{ap}}=2$) to capture shocks. Figure 3a shows the solution (solid
line) to the 1D shock-tube problem with initial conditions to the left and
right of the diaphragm (located at $x=0.8$) as follows:

$(\rho_{\mathrm{L}},u_{\mathrm{L}},p_{\mathrm{L}})=(8,0,10/\gamma),\hskip
14.22636pt(\rho_{\mathrm{R}},u_{\mathrm{R}},p_{\mathrm{R}})=(1,0,1/\gamma),$
(48)

with zero velocity everywhere and $\gamma=1.4$. The computed solution is
accurate; however, small oscillations are present in the post-shock region

Shu & Osher (1989) introduced the problem of a Mach 3 shock propagating
through density waves as way of testing a method’s ability to both capture
shocks and resolve non-shock wavy features without excessive dissipation. The
initial condition is:

$(\rho,u,p)=\begin{cases}(3.857143,2.629369,10.33333),&x<4;\\\ (1+0.2\sin
5z,0,1),&x\geq 4.\end{cases}$ (49)

The result of this test is shown in Figure 3 where the baseline comparison
(dashed line) was obtained using a fifth-order WENO scheme with $n_{z}=1600$
points, the Roe flux, and reconstruction along characteristics. The present
code performs very well.

###  4.3 Kelvin-Helmholtz instability

Here we consider the benchmark for viscous Kelvin-Helmholtz instability with a
density gradient starting with smooth initial conditions constructed by
Lecoanet et al. (2016, hereafter L2016). The full domain size is $L_{x}\times
L_{z}=1\times 2$, however, only the lower half of the $z$ domain will be
shown, since the rest is shift-symmetric. The resolution is $n_{x}\times
n_{z}$ with $n_{z}=2n_{x}$. It should be noted that L2016 write the heat
conductivity as $k=\rho\chi$, where $\chi$ is the heat diffusivity. The actual
definition (which Padé uses) is $k=\rho\chi c_{p}$. Therefore, to match L2016
we needed to divide our $k$ by $c_{p}$ which equals $5/2$ for the set-up of
$L1016$, which assumes that the gas constant $R=c_{p}-c_{v}=1$ and
$\gamma\equiv c_{p}/c_{v}=5/3$. We chose time step $\Delta t$ so that the
Courant-Friedrichs-Lewy number was 1.5. We initially applied the Padé filter
after every step with $\epsilon_{\mathrm{filter}}=0.0025$, however, the
associated matrix becomes ill-conditioned for very small values of
$\epsilon_{\mathrm{filter}}$ when the number of grid points is sufficiently
large. To alleviate this, the strength of the Padé filter was applied every 40
time steps with $\epsilon_{\mathrm{filter}}=0.10$. The simulation was run on a
laptop with an Apple M2 Pro chip using 8 cpus. The cpu time per step was 0.3
s.

Figure 4: Density field at $t=6$ for the Kelvin-Helmholtz instability
benchmark of Lecoanet et al. (2016). (a) Padé code with $n_{x}=1024$. (b)
Pseudo-spectral Dedalus code with $n_{x}=4096$. Data provided by D. Lecoanet.

Figure 4a shows the density field at $t=6$ obtained from Padé with
$n_{x}=1024$. It is compared with the result of L2016 obtained using their
pseudo spectral code Dedalus with $n_{x}=4096$; we thank Prof. D. Lecoanet
(Northwestern Univ.) for sending us the data. The agreement is very good.

###  4.4 Taylor-Couette flow

A natural test case in cylindrical coordinates is Taylor-Couette flow, i.e.,
viscous flow driven by rotating inner and outer cylinders with radii
$r_{\mathrm{i}}$ and $r_{\mathrm{o}}$, respectively. The corresponding
rotation rates are $\Omega_{\mathrm{i}}$ and $\Omega_{\mathrm{o}}$ giving
corresponding rotational speeds
$U_{\mathrm{i}}\equiv\Omega_{\mathrm{i}}r_{\mathrm{i}}$ and
$U_{\mathrm{o}}\equiv\Omega_{\mathrm{o}}r_{\mathrm{o}}$. Two of the four non-
dimensional parameters are the inner and outer Reynolds numbers
$\mathrm{Re}_{\mathrm{i}}\equiv U_{\mathrm{i}}d/\nu$ and
$\mathrm{Re}_{\mathrm{o}}\equiv U_{\mathrm{o}}r_{\mathrm{i}}d/\nu$, where
$d=r_{\mathrm{o}}-r_{\mathrm{i}}$ is the gap width. The third and fourth non-
dimensional parameters are the ratio $\eta\equiv
r_{\mathrm{i}}/r_{\mathrm{o}}$, and the non-dimensional vertical domain size
$\lambda\equiv L_{z}/d$.

The present code solves the compressible equations while most simulations
reported in the literature are for incompressible flow. To approximate
incompressible simulations the rotational Mach number of the inner cylinder is
set to $0.1$, the equation of state is isothermal, and isothermal boundary
conditions are applied at the two walls. The initial density is set to a
uniform value $\rho_{0}$. Code units are such that
$U_{\mathrm{i}}=d=\rho_{0}=1$.

Torques, $G_{\mathrm{i}}$ and $G_{\mathrm{o}}$, per unit axial length exerted
on the fluid by the inner and outer cylinders, respectively, are computed as
diagnostics.

$\mathrm{torque}=\mathrm{area\ averaged\ shear\ stress}\times
r\times\mathrm{area}.$ (50)

The area averaged shear stress is given by

$\tau_{r\phi}=\left<\mu r\frac{\partial}{\partial
r}\left(\frac{u_{\phi}}{r}\right)\right>,$ (51)

where $\mu$ is the dynamic viscosity and $\left<.\right>$ denotes an average
over the surface. Therefore,

$\displaystyle G_{\mathrm{i}}$ $\displaystyle=-2\pi
r_{\mathrm{i}}{{}^{3}}\left<\frac{\partial}{\partial
r}\left(\mu\frac{u_{\theta}}{r}\right)\right>_{r=r_{\mathrm{i}}},$ (52)
$\displaystyle G_{\mathrm{o}}$ $\displaystyle=+2\pi
r_{\mathrm{o}}{{}^{3}}\left<\frac{\partial}{\partial
r}\left(\mu\frac{u_{\theta}}{r}\right)\right>_{r=r_{\mathrm{o}}}.$ (53)

The negative sign in (52) arises from the fact that at $r=r_{\mathrm{i}}$, the
normal to the fluid surface is in the $-r$ direction. For steady flow,
conservation of angular momentum implies that $G_{\mathrm{i}}$ and
$G_{\mathrm{o}}$ must be equal and opposite in sign.

$\eta\equiv r_{\mathrm{i}}/r_{\mathrm{o}}$ | $L_{z}/d$ | $\mathrm{Re}_{\mathrm{i}}$ | $G_{\mathrm{i}}$ | $G_{\mathrm{o}}$ | Published $G_{\mathrm{i}}$ | Ref.
---|---|---|---|---|---|---
0.875 | 2.5 | 139.22 | 3.3485 | -3.3482 | 3.3539 | Marcus (1984)
0.50 | 1.988 | 78.8 | 1485 | -1485 | 1487 | Moser et al. (1983)
Table 1: Computed torques $G_{\mathrm{i}}$ and $G_{\mathrm{o}}$ per unit axial
length exerted on the fluid by the inner and outer cylinders, respectively,
for steady axisymmetric Taylor-Couette flow. These values are compared with
values obtained in the cited references. The values are normalized differently
in the two references as described in the text.

We first consider two axisymmetric cases that produce a steady flow with
counter-rotating vortices. The inner cylinder rotates at angular velocity
$\Omega_{\mathrm{i}}$ while the outer cylinder is fixed. Each case was run
using a $32\times 32$ grid ($n_{r}\times n_{z}$) with
$\epsilon_{\mathrm{filter}}=0.005$. The first two entries of Table 1 shows
that the computed torque per unit length agrees with previous published
results. Marcus (1984) uses units in which $\rho_{0}=d=\Omega
r_{\mathrm{i}}=1$. Moser et al. (1983) use the same units but normalize $G$ by
$\rho_{0}\nu^{2}$ where $\nu=\mu/\rho_{0}$. The values in Table 1 use the same
conventions.

Finally, a case of 3D unsteady counter-rotating Taylor-Couette flow is
considered following Dong (2008). The inner and outer cylinders rotate
counter-clockwise and clockwise, respectively with
$\mathrm{Re}_{\mathrm{i}}=-\mathrm{Re}_{\mathrm{o}}=500$, $\eta=0.5$ and
$L_{z}/d=2\pi$. Dong (2008) defines the non-dimensional torque coefficient for
the inner cylinder as

$\left(C_{\mathrm{T}}\right)_{\mathrm{i}}=\frac{G_{\mathrm{i}}}{\frac{1}{2}\pi\rho_{0}U_{\mathrm{i}}^{2}r_{\mathrm{i}}^{2}L_{z}},$
(54)

and similarly for the outer cylinder. The number of grid points is $48\times
96^{2}$ ($n_{r},n_{z},n_{\phi}$) and the strength of the Padé filter was set
to $\epsilon_{\mathrm{filter}}=0.005$.

Figure 5: Counter-rotating Taylor-Couette flow at
$\mathrm{Re}_{\mathrm{i}}=-\mathrm{Re}_{\mathrm{o}}=500$ following Dong
(2008). The mesh is $48\times 96^{2}(n_{r}\times n_{z}\times n_{\phi})$. (a) &
(b) Radial velocity contours at $U_{\mathrm{i}}t/d=798.8$. (a) The meridional
plane $\phi=0$. (b) Horizontal plane at mid-height ($z=\pi$) to show that the
flow is three-dimensional. (c) Torque coefficients (black) compared with
curves (red) digitized from Dong (2008, Figure 3) and shifted backward in time
by $U_{\mathrm{i}}\Delta t/d=802$. The negative of the torque coefficient is
plotted for the outer cylinder.

Figures 5a and b shows the radial velocity in a meridional and horizontal
plane and reveals the three-dimensionality of the flow. Figure 5c shows the
torque coefficients for the inner and outer cylinders after the flow has
reached statistical stationarity. The values agree with those shown in Figure
3 of Dong (2008).

###  4.5 Vertical shear instability and the effect of varying the strength of
the Padé filter

Here we present results for axisymmetric vertical shear instability (VSI) at
an early stage of evolution following the set up of Nelson et al. (2013).
Detailed results will be presented in a forthcoming publication which will
also include 3D results.

Parameter | Value
---|---
Orbital period, $T_{0}$, at $r_{0}$ | $1$
Density, $\rho_{0}$, at $r_{0}$ | $1$
Scale height, $H_{0}$, at $r_{0}$ | $1$
Density exponent, $p$ | $-3/2$
Temperature exponent, $q$ | $-1$
Disk aspect ratio, $H_{0}/r_{0}$ | $0.10$
Radial domain (including sponge), $[r_{\mathrm{min}}/H_{0},r_{\mathrm{max}}/H_{0}]$ | $[6.5,13.5]$
Vertical domain (including sponge), $[z_{\mathrm{min}}/H_{0},z_{\mathrm{max}}/H_{0}]$ | $[-3.5,3.5]$
Width of sponge at domain border, $\delta_{\mathrm{sponge}}/H_{0}$ | 0.5
Decay period, $t_{\mathrm{sponge}}$, for sponge, | 20 time steps
Number of grid points, $n_{r}\times n_{z}$ | $512^{2}$
Strength of Padé filter, $\epsilon_{\mathrm{filter}}$ | 0.01-0.125
Table 2: Parameters for the axisymmetric vertical shear instability run.
$H_{0}$ is the disk scale height at mid-radius and was set to unity. Subscript
‘0’ refers to a quantity evaluated at mid-radius.

The parameters of the set-up are given in Table 2. A locally isothermal
equation of state

$p=\rho c_{\mathrm{i}}^{2}(r),$ (55)

is used which represents the case of infinitely rapid relaxation of
temperature to the basic state. The square of the sound speed, which is
proportional to the temperature, is specified as a power-law:

$c_{\mathrm{i}}^{2}(r)=c_{0}^{2}\left(r/r_{0}\right)^{q},$ (56)

where the temperature exponent $q=-1$, $r_{0}$ is the mid-radius of the
computational domain, and $c_{0}$ is chosen to make the scale height
$H_{0}\equiv H(r_{0})=1$.

Zero normal velocity boundary conditions are applied at all four domain edges.
Each edge abuts a sponge region of width $\delta_{\mathrm{sponge}}=0.5H_{0}$
where the flow relaxes back to the basic state with a characteristic period
$t_{\mathrm{sponge}}$ given in the table.

Figure 6: Effect of varying the strength of the Padé filter,
$\epsilon_{\mathrm{filter}}$ when the Padé filter is used as an implicit sub-
grid scale treatment for an axisymmetric VSI simulation. Completely white
pixels are where $\omega_{\phi}$ exceeds the range of the legend. The mesh
size is $512^{2}$. The plots in the right-hand column are intended to show the
level of numerical $2\Delta$ oscillations which appear as sawtooth features.

Figure 6 shows the azimuthal vorticity, $\omega_{\phi}$ at $t/T_{0}=37$ which
is just after saturation of the linear phase of the instability. We caution
the reader that the statistically stationary state is quite different and
reached much later at $t/T_{0}\approx 350$; this will be reported in a later
publication. The result of using three different values for the strength,
$\epsilon_{\mathrm{filter}}$, of the Padé filter is shown.

Figure 7: (a) Vertical velocity for the axisymmetric vertical shear
instability simulation using a $512^{2}$ mesh. The units of velocity are scale
heights ($H_{0}$) per orbit time ($T_{0}$). (b) Similarly, the perturbation
azimuthal velocity.

The azimuthal vorticity consists of pairs of shear layers of opposite sign
which induce up and down jets of vertical velocity (Figure 7a) that are
characteristic of VSI (Nelson et al., 2013). The shear layers roll up into
discrete eddies via the Kelvin-Helmholtz instability. It is important to note
that VSI also produces shear layers with $\partial\delta u_{\phi}/\partial r$
shear, i.e., radial gradients of perturbation angular velocity. These are
shown in Figure 7b which plots the azimuthal velocity perturbation, $\delta
u_{\phi}$, and shows that is it about $2/3$ the value of the vertical velocity
perturbation. This shear will also be subject to Kelvin-Helmholtz instability,
however, it will be modified by the presence of mean Keplerian shear.

We now discuss how a user should choose the filter strength,
$\epsilon_{\mathrm{filter}}$. Direct numerical simulations (DNS) of turbulent
flow are performed with molecular viscosity and all the scales of the
turbulence down to the dissipation scale are well resolved. In this case, the
very minimum value of $\epsilon_{\mathrm{filter}}$ should be chosen. For
instance in the Taylor-Couette simulations we chose
$\epsilon_{\mathrm{filter}}=0.005$.

However, the Reynolds number in protoplanetary disks is too large for
numerical simulations to be able to resolve all the turbulent scales.
Therefore, some treatment of the unresolved scales is required. For
engineering and geophysical flows, the most common practice is to use an
explicit model for the sub-grid stresses, the simplest being the Smagorinsky
model. Another approach, followed for all protoplanetary disk simulations to
date, is to simply let the dissipation inherent in the numerical method damp
scales near the grid cut-off. This procedure is referred to as implicit large-
eddy simulation (ILES) and was first articulated by Boris et al. (1992).
Comparison with direct numerical simulations (DNS, in which all scales are
resolved) have since shown that it is accurate (Ritos et al., 2018). In our
approach, we use the dissipation provided by the Padé filter as an implicit
sub-grid treatment. To leading order, the Padé filter corresponds to a fourth-
order hyperviscosity; see equation (30) in Shariff & Wray (2018) and the
discussion following it.

When the Padé filter is used an ILES treatment, $\epsilon_{\mathrm{filter}}$
should be chosen to balance the desire to capture as wide a range of small
scales possible (with a fixed grid size) by lowering
$\epsilon_{\mathrm{filter}}$, while at the same not allowing to much energy to
to pile-up in $2\Delta$ waves. For illustration, Figure 6 shows the result of
varying the strength, $\epsilon_{\mathrm{filter}}$, of the filter on the
vorticity field. The reader may refer back to Figure 1b which shows the filter
transfer function for different $\epsilon_{\mathrm{filter}}$ choices. The left
hand column of plots in Figure 6 shows that with reduced
$\epsilon_{\mathrm{filter}}$ more finer scales of the flow are resolved. The
right-hand column of plots zoom in to a square region with sides equal to $H$.
For the lowest filter strength ($\epsilon_{\mathrm{filter}}=0.01$), $2\Delta$
(sawtooth) oscillations can be observed in thin vortex layers oriented at
$45^{\circ}$ to the mesh. For $\epsilon_{\mathrm{filter}}=0.04$ and 0.125,
smaller amplitude oscillations are present in one vorticity layer whose width
is about one grid diagonal, which is very thin indeed. These oscillations are
not visually detectable in plots of the velocity field even for
$\epsilon_{\mathrm{filter}}=0.01$. We conclude that $0.125$ would be a
conservative choice for $\epsilon_{\mathrm{filter}}$.

Finally, we would like to discuss some physical effects that manifest as
diffusivity represented by the filter is reduced and the effective resolution
is increased. (a) The rolled-up vortices are smaller and there are more of
them. This is explained as follows. The shear-layer thickness, $\delta$, is
reduced with smaller diffusivity. The most amplified Kelvin-Helmholtz (KH)
wavelength $\lambda\approx 2\pi\delta$ is therefore also reduced, and with it
the size of the vortices. The number of vortices increases because there are
more waves per unit length. (b) The shear-layer vorticity increases. We have
that $\omega_{\phi}\sim\Delta U/\delta$ where $\Delta U$, the jump in velocity
across each shear layer, is independent of diffusivity. Therefore a reduction
in $\delta$ with diffusivity leads to increased $\omega_{\phi}$. (c) The
rolled-up vortices appear earlier. The KH growth-rate $\propto\Delta
Uk_{\mathrm{max}}$, where $k_{\mathrm{max}}$ is the most amplified wavenumber.
Since $k_{\mathrm{max}}$ increases with reduced thickness, the KH vortices
develop earlier with reduced diffusivity. (d) The vorticity in the the rolled-
up vortex cores increases. For each vortex we have
$\omega_{\phi}\sim\Gamma/\mathrm{area}$, where $\Gamma\sim\Delta U\lambda$ is
its circulation. Since the vortex area $\sim\lambda^{2}$, we get that
$\omega_{\phi}\sim\Delta U/\lambda$ which increases since $\lambda$ decreases.

###  4.6 3D vertical shear instability: comparison of Fargo and non-Fargo runs

In Shariff & Wray (2018) a comparison (with plots of error) was made between
runs with and without the Fargo treatment for the case of two co-rotating
vortices in a razor thin disk. Here, we perform a similar comparison for 3D
vertical shear instability (VSI). The simulation was first run till
$t/T_{0}=300.27$ with Fargo activated. Next runs were made with and without
Fargo for one orbital period ($T_{0}$) at mid-radius. The run parameters are
the same as for the axisymmetric run presented in Table 2. The only
differences are a resolution of $512\times 512\times 1024$ ($n_{r}\times
n_{z}\times n_{\phi}$) with an azimuthal domain of $\phi\in[0,2\pi)$, and a
Padé filter strength of $\epsilon_{\mathrm{filter}}=0.125$. For the non-Fargo
run, the Padé filter was applied every other time step due to the fact that
this run required about twice as many time steps as the run with Fargo.

Figure 8: A portion of the horizontal plane $(r,\phi)$ cutting through the
disc at $z/H_{0}=2$ comparing perturbation vertical vorticity,
$\delta\omega_{z}$, for a 3D VSI run with and without Fargo. Each run was
started at $t/T_{0}=300.27$ and run for $0.5T_{0}$ up to $t/T_{0}=300.77$. (a)
With Fargo. (b) Without Fargo.

Figure 8 compares the vertical vorticity ($\delta\omega_{z}$) perturbation
(relative to the basic state) in the $z/H_{0}=2$ plane after a time of
$0.5T_{0}$, i.e., at $t/T_{0}=300.77$. Here $T_{0}$ is the orbital period at
the mid-radius of the computational domain. The difference in the two
solutions at this time is very small. Due to the chaotic nature of the flow,
the error due to the different time steps, $\Delta t$, chosen for the two
simulations grows with time and differences become more apparent. The CFL
number was chosen to be 1.5 and based on this, the time step selected by the
code for the non-Fargo run was $\Delta t/T_{0}=2.25\times 10^{-4}$. This
choice was constrained by Keplerian advection. For the Fargo run, the code
selected $\Delta t/T_{0}=5.49\times 10^{-4}$ for the same CFL number which
represents a better than factor of two improvement. In the run with Fargo, the
choice of time step was constrained by the characteristic-wave speed and grid
size in the radial direction. The cpu time for the non-Fargo and Fargo runs
was $1.50$ and $1.76$ seconds per step. This represents a 17.3% overhead for a
more than factor of two gain in time step. The Intel Broadwell processor was
used for these runs.

We close with a brief description of the physics observed in the 3D VSI runs.
The Keplerian mean flow is counter-clockwise in Figure 8 and the vertical
vorticity perturbation consists of layers of cyclonic $\delta\omega_{z}$ (red
bands) which induce across them, a positive jump in specific angular momentum,
$j_{\phi}=u_{\phi}r$ as $r$ increases. These layers are formed by the vertical
transport of basic state angular momentum by the vertical jets that are the
main feature of VSI. This will be discussed in more detail in a forthcoming
paper. In between the cyclonic layers, one observes weaker and more diffuse
anti-cyclonic $\delta\omega_{z}$ which reduces $j_{\phi}$ compared to the
basic state. In the inner portion of the disk, anti-cyclonic
$\delta\omega_{z}$ structures take the form of smaller aspect ratio
structures.

###  4.7 Parallel scaling tests for 3D vertical shear instability runs

Following the suggestion of the referee, the efficiency $\eta$ of the code is
defined as the ratio of the useful work to the resources consumed. The useful
work is the number $N_{\mathrm{grid}}$ of grid points evolved while the
resources consumed is the cpu time $t_{\mathrm{cpu}}$ (per time step in secs.)
times the number of cores, $N_{\mathrm{cores}}$. We also include the ratio
$\Delta t_{\mathrm{actual}}/\Delta t_{\mathrm{Fargo}}$ to account for the
reduced time step in non-Fargo runs:

$\eta\equiv\frac{N_{\mathrm{grid}}/2^{17}}{t_{\mathrm{cpu}}N_{\mathrm{cores}}}\frac{\Delta
t_{\mathrm{actual}}}{\Delta t_{\mathrm{Fargo}}},$ (57)

where we have normalized the number of grid points to $2^{17}$. When there is
perfect scaling, $\eta$ should be a constant as $N_{\mathrm{cores}}$ is
increased.

All the tests in this sub-section were performed using Intel Haswell nodes
which use the E5-2680V3 (Xeon) processor. The variable MPI_IB_RAILS used by
the InfiniBand network was set to 2; this causes two InfiniBand (IB) fabrics
to be used for communication and results in reduced cpu time. Figures 9a and b
plot $\eta$ for a 3D vertical shear instability (VSI) setup with Fargo
activated; therefore the ratio $\Delta t_{\mathrm{actual}}/\Delta
t_{\mathrm{Fargo}}$ in (57) equals unity.

Figure 9: Scaling tests for the 3D vertical shear instability (VSI) set-up
performed using the Intel Haswell nodes. The environment variable MPI_IB_RAILS
used by the InfiniBand network was set to 2 for all runs. (a) Efficiency,
$\eta$, for the weak scaling test where the problem size per processor is
fixed. (b) Efficiency for the strong scaling test where the total problem size
is fixed at $n_{r}\times n_{z}\times n_{\phi}=512\times 512\times 1024$. (c)
Efficiency, $\eta$, with Fargo deactivated. The measure, $\eta$, accounts for
the decrease in time step when Fargo is deactivated.

Figure 9a is for a weak scaling test in which the grid size per core is kept
fixed while $N_{\mathrm{cores}}$ is increased. In other words, both
$N_{\mathrm{cores}}$ and the total grid size increase simultaneously. In the
present test, the number of grid points for the one core run is $64\times
32\times 64$ ($n_{r}\times n_{z}\times n_{\phi}$) and each direction is
successively doubled in resolution as $N_{\mathrm{cores}}$ doubles. Figures 9a
plots two curves, one which used the actual cpu time (solid) and another
(dashed) for which the time taken to perform transposes was subtracted out.
Both curves show an initial rapid decrease in $\eta$ up to
$N_{\mathrm{cores}}=8$. The main cause of this is likely the fact that each
Xeon processor in a Haswell node has 12 cores which share a single memory and
level 3 cache. Contention for both resources increases as the number of cores
increases from 1 to 12. Therefore, we focus on the region
$N_{\mathrm{cores}}\geq 12$. At $N_{\mathrm{cores}}=2048$, the efficiency has
decreased to 38% of its value at $N_{\mathrm{cores}}=16$. When the cpu time
for transposes is subtracted out, the efficiency remains relatively flat. This
indicates that the loss in efficiency is due to communication intensive
transposes. Indeed, the fraction of time (not shown) taken for transposes
increases from $0.38$ to $0.73$ as $N_{\mathrm{cores}}$ increases from 16 to
2048.

Figure 9b is for a strong scaling test in which the total problem size is
fixed at $512\times 512\times 1024$ ($n_{r}\times n_{z}\times n_{\phi}$) and
the number of cores is varied. At $N_{\mathrm{cores}}=2048$, the efficiency
has decreased to 67% of its value at $N_{\mathrm{cores}}=64$. The cpu fraction
taken for transposes increases from 0.57 to 0.73 (not shown) in this range.
When the cpu time for transposes is subtracted out, the efficiency increases
slowly. This is explained as follows. In the strong test, the total number of
memory fetches is constant with increasing $N_{\mathrm{cores}}$, however, the
number of cache hits (when needed data is found to be already in cache) is
likely to be statistically higher since the number of cache slots per fetch is
higher with more cores.

Figure 9c compares the efficiency of the Fargo versus the non-Fargo scheme. It
shows that there is at least a factor of 3.7 advantage to using Fargo, i.e.,
the overhead of the Fargo method is more than compensated by an increase in
time step.

In conclusion, the all-to-all communication of transposes leads to a
significant loss in efficiency as the number of cores is increased. To reduce
this communication overhead, an effort is underway to use a parallel Padé
algorithm (Kim et al., 2021).

##  5 Closing Remarks

A code has been developed that uses a fourth-order Padé scheme to simulate
hydrodynamic turbulence in protoplanetary disks. Padé schemes are non-
dissipative and have high resolving power. Thus, with the same grid
resolution, they are better able to capture fine scale vortical features
compared to the dissipative shock capturing schemes employed in most
astrophysics codes. They also have better resolving power then central finite
difference schemes of the same order.

Suggested improvements are as follows. (i) To eliminate communication
intensive transposes, consider using parallel tridiagonal matrix algorithms
(Kim et al., 2021, and the references therein) which require much less
communication. Kim et al. (2021) demonstrate good scaling as the number of
cores is increased. (ii). For simulations that require long radial domains
(more than $\approx 6$ scale heights) it would be better to use spherical
rather than cylindrical coordinates. An option for this could be provided.
(iii) The sixth-order tridiagonal Padé scheme

$\alpha f^{\prime}_{j-1}+f^{\prime}_{j}+\alpha
f^{\prime}_{j+1}=\frac{a}{2h}\left(f_{j+1}-f_{j-1}\right)+\frac{b}{4h}\left(f_{j+2}-f_{j-2}\right),$
(58)

with $\alpha=1/3$, $a=14/9$, and $b=1/9$ could be implemented. (iv) The
capability to track Lagrangian solid particles could be provided. (v) An
intercomparison effort with other codes could be performed.

The code is available at https://github.com/NASA-Planetary-Science/Pade-disk-
code and https://zenodo.org/records/11114378. The citation for Zenodo is
Shariff (2024).

I am grateful to Alan Wray (NASA Ames) for providing several routines that
have been used in Padé, including directional transposes, tridiagonal solvers,
and the Padé filter. He also helped to explain features observed in the
scaling tests. I am grateful to Jeff Cuzzi (NASA Ames) for performing the
internal review and for his encouragement, Debanjan Sengupta (New Mexico State
Univ., Las Cruces) for performing the internal review and suggesting the weak
scaling test; Orkan Umurhan (SETI Inst.) for his encouragement and advice on
the vertical shear instability runs, Prof. Ali Mani (Stanford Univ.) for
advice on the implementation of artificial bulk viscosity, and Prof. Daniel
Lecoanet (Northwestern Univ.) for providing the result of the Kelvin-Helmholtz
test case run with Dedalus. Finally, I thank the referee for many useful
suggestions that improved the paper. This work was partly supported by NASA’s
Internal Scientist Funding Model (ISFM) through the “Origins” group in NASA
Ames’ Space Science Division.

##  Appendix A Molecular viscous force, viscous heating, and heat conduction

The code provides the option to add terms that represent molecular viscosity
and heat conduction. Similar terms arise in models of sub-grid turbulence. For
these we have coded the models due to Smagorinsky (1963) and Vreman (2004).
However, since we have not tested them, this section describes the
implementation for the molecular/laminar case only. If the Fargo option has
been activated, Fargo chain-rule is applied wherever needed.

###  A.1 Viscous force

The viscous stress tensor is

$T_{ij}=2\mu
S_{ij}+\left(\mu_{\mathrm{b}}-\frac{2}{3}\mu\right)S_{kk}\delta_{ij},$ (A1)

where $\mu$ and $\mu_{\mathrm{b}}$ are the shear and bulk viscosities,
respectively, and

$S_{ij}=\frac{1}{2}\left(\partial_{j}u_{i}+\partial_{i}u_{j}\right)$ (A2)

is the strain rate tensor.

Equation (A1) is implemented in subroutine laminar_stress_and_heat_flux. The
components of the strain tensor (which is symmetric) are computed in
subroutine strain_tensor as follows

$\displaystyle S_{zz}=\partial_{z}u_{z},\hskip 7.11317ptS_{\phi
z}=\frac{1}{2}\left(\partial_{z}u_{\phi}+\frac{1}{r}\partial_{\phi}u_{z}\right),\hskip
7.11317ptS_{zr}=\frac{1}{2}\left(\partial_{r}u_{z}+\partial_{z}u_{r}\right),$
(A3) $\displaystyle S_{rr}=\partial_{r}u_{r},\hskip
7.11317ptS_{r\phi}=\frac{1}{2}\left(\frac{1}{r}\partial_{\phi}u_{r}+\partial_{r}u_{\phi}-\frac{u_{\phi}}{r}\right),\hskip
7.11317ptS_{\phi\phi}=\frac{1}{r}\partial_{\phi}u_{\phi}+\frac{u_{r}}{r}.$
(A4)

The above expressions are from Aris (1989, p.181) and agree with
Batchelor(1967, p. 602).

The viscous force is the divergence of the viscous stress tensor and given by
Aris (1989, p. 179) for orthogonal coordinates as follows (in his notation):

$F^{\mathrm{v}}_{i}=T(ij,j)=\frac{h_{i}}{h_{1}h_{2}h_{3}}\frac{\partial}{\partial
x_{j}}\left[\frac{h_{1}h_{2}h_{3}}{h_{i}h_{j}}T(ij)\right]+\frac{h_{i}}{h_{j}h_{k}}\left\\{\begin{subarray}{c}i\\\
jk\end{subarray}\right\\}T(jk),\hskip 7.11317pt(\mathrm{no\ sum\ on\ }i).$
(A5)

Here $(x_{1},x_{2},x_{3})=(r,\phi,z)$ and corresponding scale factors are
$h_{1}=h_{r}=1$, $h_{2}=h_{\phi}=r$, and $h_{3}=h_{z}=1$. The quantities in
braces are Christoffel symbols and the only non-zero ones are

$\left\\{\begin{subarray}{c}2\\\
12\end{subarray}\right\\}=\left\\{\begin{subarray}{c}2\\\
21\end{subarray}\right\\}=r\mathrm{\ and\ }\left\\{\begin{subarray}{c}1\\\
22\end{subarray}\right\\}=-r.$ (A6)

The symbolic algebra package maxima was used to verify the correctness of
Aris’ expression (A5) by writing Cartesian velocities in terms of cylindrical
quantities (velocities and coordinates) and computing Cartesian viscous forces
in terms of cylindrical quantities using the chain rule throughout. These can
be be rotated to obtain the forces in cylindrical coordinates in terms of
cylindrical quantities and compared with (A5). The relevant maxima script can
be found in check_Aris.mac in the Symbolic_algebra sub-directory.

The final expressions for the viscous force area:

$\displaystyle F^{\mathrm{v}}_{z}$ $\displaystyle=$
$\displaystyle\frac{1}{r}\left[\partial_{r}(rT_{zr})+\partial_{\phi}(T_{z\phi})+\partial_{z}(rT_{zz})\right],$
(A7) $\displaystyle F^{\mathrm{v}}_{r}$ $\displaystyle=$
$\displaystyle\frac{1}{r}\left[\partial_{r}(rT_{rr})+\partial_{\phi}(T_{r\phi})+\partial_{z}(rT_{rz})\right]-\frac{1}{r}T_{\phi\phi},$
(A8) $\displaystyle F^{\mathrm{v}}_{\phi}$ $\displaystyle=$
$\displaystyle\partial_{r}(T_{\phi
r})+\frac{1}{r}\partial_{\phi}(T_{\phi\phi})+\partial_{z}(T_{\phi
z})+\frac{2}{r}T_{\phi r}.$ (A9)

These agree with the expressions given on p. 739B of Bird et al. (1960).

###  A.2 Viscous heating

The equation for total energy $e=e_{\mathrm{int}}+\rho u_{i}u_{i}$ (per unit
volume) has a term for the work done by shear stresses (Liepmann & Roshko,
2001, p. 335),

$W^{\mathrm{shear}}=\frac{\partial}{\partial x_{j}}(T_{ij}u_{i}).$ (A10)

For the internal energy equation which we solve, we must subtract the kinetic
energy dissipation:

$D=u_{i}\frac{\partial}{\partial x_{j}}T_{ij}.$ (A11)

This gives the viscous heating term for the internal energy

$Q^{\mathrm{v}}=W^{\mathrm{shear}}-D=T_{ij}\frac{\partial u_{i}}{\partial
x_{j}}.$ (A12)

Now $T_{ij}$ is symmetric so only the symmetric part of $\partial_{j}u_{i}$
survives. Then substituting the constitutive equation (A1) for $T_{ij}$ into
(A12) one obtains

$Q^{\mathrm{v}}=2\mu
S_{ij}S_{ij}+\left(\mu_{\mathrm{b}}-\frac{2}{3}\mu\right)S_{kk}^{2}.$ (A13)

###  A.3 Heat conduction

The flux of internal energy due to molecular conductivity is given by
Fourier’s law:

$\vec{q}_{\mathrm{cond}}=-\frac{k}{c_{v}}\nabla(c_{v}T),$ (A14)

where we have multiplied and divided by $c_{v}$ inside and outside the
gradient, respectively. This assumes that $c_{v}$ is constant, i.e., that we
have a calorically perfect gas. Now $c_{v}T$ is simply
$e_{\mathrm{int}}/\rho$, and using the definition of the Prandtl number
$\mathrm{Pr}\equiv\mu c_{p}/k$ we get

$\frac{k}{c_{v}}=\frac{\mu\gamma}{\mathrm{Pr}}.$ (A15)

##  Appendix B Discrete conservation

###  B.1 Theory

Here we describe how to choose boundary schemes to ensure that the overall
scheme possesses a discrete conservation property. We follow Lele (1992) and
Brady & Livescu (2019) and offer two clarifications: (1) There is a
distinction between provisional and final weights; the weights given in §4.2
of Lele (1992) are provisional. (2) The weights cannot be specified a priori
and must be determined as part of the solution. Specifically, one needs to
verify that the final weights provide a reasonable discrete conservation law.
In general, the final weights will not correspond exactly to a quadrature
rule.

For a system of partial differential equations in more than one dimension, we
compute derivatives of fluxes along each direction separately. Hence, it is
sufficient to consider the one-dimensional partial differential equation

$\frac{\partial u}{\partial t}+\frac{\partial f}{\partial x},\hskip
14.22636ptx\in[0,L].$ (B1)

Upon integration over the domain, (B1) gives the conservation law:

$\frac{d}{dt}\int_{0}^{L}u(x,t)\,dx=f(0)-f(L).$ (B2)

Padé differencing applied to (B1) should possess a discrete analog of (B2); in
its absence, a long time solution can drift and fail to achieve statistical
stationarity. Padé difference schemes have the form

$A\mathbf{f}^{\prime}=B\mathbf{f}.$ (B3)

where $A$ and $B$ are banded matrices and henceforth, lower case bold letters
will be used to denote column vectors.

We now state a result that was stated by Lele (1992) without a proof, which
was later provided by Brady & Livescu (2019).

######  Proposition 1.

To obtain a discrete analog of (B2), columns 2 through $N-1$ of matrix $B$
must have a weighted sum of zero, i.e.,

$\mathbf{w}^{T}\mathbf{b}_{i}=0\hskip 7.11317pt\mathrm{for\ }i\in[2,N-1],$
(B4)

where $\mathbf{w}$ is a column vector of weights and $\mathbf{b}_{i}$ is the
$i$th column of matrix $B$. The weights $\mathbf{w}$ are provisional; final
weights will be given below.

######  Remark.

Not all the weights, $\mathbf{w}$, can be specified a priori but must be
obtained as part of the process of satisfying (B4).

###### Proof.

We assume that grid points $x_{i},i=1,\ldots,N$ have been laid out according
to a smooth analytic mapping $x(\xi)$ such that $\xi_{i}=i$, i.e., $\xi$ is a
smooth grid index variable. In implementations, it is convenient to take
derivatives with respect to $\xi$ and then use the chain rule:

$\left(\frac{\partial f}{\partial x}\right)_{i}=\left(\frac{\partial
f}{\partial\xi}\right)_{i}\left(\frac{d\xi}{dx}\right)_{i}.$ (B5)

For simplicity we use the notation

$f_{i}^{\prime}(\xi)\equiv\left(\frac{\partial
f}{\partial\xi}\right)_{i},\hskip
14.22636pth_{i}^{-1}\equiv\left(\frac{d\xi}{dx}\right)_{i}.$ (B6)

Then the spatially discretized version of (B2) is

$\frac{du_{i}}{dt}+\frac{1}{h_{i}}\sum_{j,k}\left[A^{-1}\right]_{ij}B_{jk}f_{k}=0.$
(B7)

Defining the vector
$\mathbf{U}=\left[h_{1}u_{1},h_{2}u_{2},\ldots,h_{N}u_{N}\right]^{T}$, we can
write (B7) as

$A\frac{d\mathbf{U}}{dt}+B\mathbf{f}=0.$ (B8)

A weighted sum is applied to (B8) to mimic the integration in (B2):

$\frac{d}{dt}\mathbf{w}^{T}A\mathbf{U}=-\mathbf{w}^{T}B\mathbf{f}.$ (B9)

Following Brady & Livescu (2019), let
$B=[\mathbf{b}_{1},\mathbf{b}_{2},\ldots,\mathbf{b}_{N}]$ where
$\mathbf{b}_{i}$ denotes the $i$th column vector of the matrix $B$. Now

$\displaystyle\mathbf{w}^{T}B$ $\displaystyle=$
$\displaystyle\mathbf{w}^{T}[\mathbf{b}_{1},\mathbf{b}_{2},\ldots,\mathbf{b}_{N}],$
(B10) $\displaystyle=$
$\displaystyle[\mathbf{w}^{T}\mathbf{b}_{1},\mathbf{w}^{T}\mathbf{b}_{2},\ldots,\mathbf{w}^{T}\mathbf{b}_{N}],$
(B11)

so that

$\mathbf{w}^{T}B\mathbf{f}=\sum_{i=1}^{N}\mathbf{w}^{T}\mathbf{b}_{i}f_{i}.$
(B13)

Hence (B9) can be written

$\frac{d}{dt}\mathbf{w}^{T}A\mathbf{U}=-\mathbf{w}^{T}\mathbf{b}_{1}f_{1}-\mathbf{w}^{T}\mathbf{b}_{N}f_{N}-\sum_{i=2}^{N-1}\mathbf{w}^{T}\mathbf{b}_{i}f_{i}.$
(B14)

In order to arrive at a discrete conservation law analogous to (B2), let us
try imposing

$\sum_{i=2}^{N-1}\mathbf{w}^{T}\mathbf{b}_{i}f_{i}=0.$ (B15)

For (B15) to be true for arbitrary $f_{i}$ we must have that

$\mathbf{w}^{T}\mathbf{b}_{i}=0\hskip 7.11317pt\mathrm{for\ }i\in[2,N-1].$
(B16)

In other words we want columns 2 through $N-1$ of the matrix $B$ to have a
weighted sum of zero.

We now show, finally, that the trial condition (B16) does indeed lead to a
discrete conservation law. Equation (B16) gives the set of conditions one uses
to solve for the provisional weights. Then (B14) becomes

$\frac{d}{dt}\widetilde{\mathbf{w}}^{T}\mathbf{U}=-\mathbf{w}^{T}\mathbf{b}_{1}f_{1}-\mathbf{w}^{T}\mathbf{b}_{N}f_{N},$
(B17)

where we have defined a new set of weights

$\widetilde{\mathbf{w}}^{T}\equiv\mathbf{w}^{T}A.$ (B18)

Now $\mathbf{w}^{T}\mathbf{b}_{1}$ is a number and if we use the same boundary
schemes and quadrature weights at the left and right boundaries then
$-\mathbf{w}^{T}\mathbf{b}_{1}=+\mathbf{w}^{T}\mathbf{b}_{N}$ and we can
divide (B17) by it. Finally, putting back $U_{i}=u_{i}h_{i}$ (B17) we get the
desired discrete conservation law

$\frac{d}{dt}\sum_{i=1}^{N}\widehat{w}_{j}h_{j}u_{j}=f_{1}-f_{N},,$ (B19)

where the final weights are

$\widehat{\mathbf{w}}=-\frac{\mathbf{w}^{T}A}{\mathbf{w}^{T}\mathbf{b}_{1}}.$
(B20)

Note that this expression is homogeneous in the provisional weights
$\mathbf{w}$. Hence we may scale the provisional weights as we wish when we
solve equations (B16) to determine them. In particular, it it is convenient to
choose them to be unity in the interior of the domain. ∎

###  B.2 Application to the present differentiation scheme

We consider matrix entries near the left end of the domain; the weights near
the right end are obtained by symmetry. Referring back to §3.4, the first five
columns of the right-hand-side matrix of the present Padé scheme are

$B=\left(\begin{matrix}a_{1}&b_{1}&c_{1}&&\\\ -\widehat{a}&0&\widehat{a}&&\\\
&-\widehat{a}&0&\widehat{a}&\\\ &&-\widehat{a}&0&\widehat{a}\\\
&&&-\widehat{a}&0\\\ &&&&-\widehat{a}\\\ \end{matrix}\right),$ (B21)

where we have defined $\widehat{a}=a/2$. The condition that the weighted sum
of the fifth column equals zero implies that $w_{4}=w_{6}$ which we set equal
to unity. The same holds true for the rest of the interior weights. The
condition on the fourth column implies that $w_{3}=w_{5}$ which we also set
equal to unity. Using the fact that $w_{3}=1$, the second column gives the
condition

$w_{1}b_{1}-\widehat{a}=0,$ (B22)

or

$w_{1}=\widehat{a}/b_{1}=(3/4)/2=3/8.$ (B23)

The third column gives

$w_{1}c_{1}+w_{2}\widehat{a}-w_{4}\widehat{a}=0.$ (B24)

Using known information, we get

$w_{2}=1-c_{1}/b_{1}=3/4.$ (B25)

We now have all the provisional weights.

The final weights are obtained from (B20) where the first four columns of the
LHS matrix for the present scheme is

$A=\left(\begin{matrix}1&\alpha_{1}&0&\\\ \alpha&1&\alpha&\\\
&\alpha&1&\alpha\\\ &&\alpha&1&\\\ \ &&&\alpha\end{matrix}\right)$ (B26)

Hence (B20) gives the final weights as

$\displaystyle\widehat{w}_{1}$
$\displaystyle=-\frac{w_{1}+w_{2}\alpha}{w_{1}a_{1}-w_{2}\widehat{a}}=\frac{3}{8}$
(B27) $\displaystyle\widehat{w}_{2}$
$\displaystyle=-\frac{w_{1}\alpha_{1}+w_{2}+w_{3}\alpha}{w_{1}a_{1}-w_{2}\widehat{a}}=\frac{7}{6}$
(B28) $\displaystyle\widehat{w}_{3}$
$\displaystyle=-\frac{w_{2}\alpha+w_{3}+w_{4}\alpha}{w_{1}a_{1}-w_{2}\widehat{a}}=\frac{23}{24}$
(B29) $\displaystyle\widehat{w}_{j}$
$\displaystyle=-\frac{w_{j-1}\alpha+w_{j}+w_{j+1}\alpha}{w_{1}a_{1}-w_{2}\widehat{a}}=1,\hskip
14.22636ptj\in[4,N-3].$ (B30)

It is interesting that $\widehat{w}_{1}+\widehat{w}_{2}+\widehat{w}_{3}=5/2$
which would be the case for trapezoidal rule and indeed the three weights
approximate the weights $1/2$, $1$, and $1$ for trapezoidal rule.

## References

  * Adams & Shariff (1996) Adams, N., & Shariff, K. 1996, J. Chem. Phys., 127, 27, doi: https://doi.org/10.1006/jcph.1996.0156
  * Aris (1989) Aris, R. 1989, Vectors, Tensors, and the Basic Equations of Fluid Mechanics (New York: Dover)
  * Batchelor (1967) Batchelor, G. 1967, An Introduction to Fluid Dynamics (Cambridge, UK: Cambridge Univ. Press)
  * Benítez-Llambay & Masset (2016) Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11, doi: 10.3847/0067-0049/223/1/11
  * Bird et al. (1960) Bird, R., Stewart, W., & Lightfoot, E. 1960, Transport phenomena (New York: John Wiley & Sons, Inc.)
  * Boris et al. (1992) Boris, J. P., Grinstein, F. F., Oran, E. S., & Kolbe, R. L. 1992, Fluid Dyn. Res., 10, 199, doi: 10.1016/0169-5983(92)90023-P
  * Brady & Livescu (2019) Brady, P., & Livescu, D. 2019, Computers & Fluids, 183, 84 , doi: https://doi.org/10.1016/j.compfluid.2018.12.010
  * Brandenburg et al. (2021) Brandenburg, A., Johansen, A., Bourdin, P., et al. 2021, J. Open Source Soft., 6, 2807, doi: 10.21105/joss.02807
  * Cook & Cabot (2005) Cook, A. W., & Cabot, W. H. 2005, J. Comp. Phys., 203, 379, doi: 10.1016/j.jcp.2004.09.011
  * Dong (2008) Dong, S. 2008, J. Fluid Mech., 615, 371–399, doi: 10.1017/S0022112008003716
  * Frigo & Johnson (2005) Frigo, M., & Johnson, S. 2005, Proceedings of the IEEE, 93, 216
  * Gustafsson (1981) Gustafsson, B. 1981, SIAM J. Num. Anal., 18, 179, doi: 10.1137/0718014
  * Harten et al. (1987) Harten, A., Engquist, B., Osher, S., & Chakravarthy, S. 1987, J. Comp. Phys., 71, 231, doi: https://doi.org/10.1016/0021-9991(87)90031-3
  * Hu et al. (2013) Hu, X., Adams, N., & Shu, C.-W. 2013, J. Comp. Phys., 242, 169, doi: https://doi.org/10.1016/j.jcp.2013.01.024
  * Kim et al. (2021) Kim, K.-H., Kang, J.-H., Pan, X., & Choi, J.-I. 2021, Comp. Phys. Comm., 260, 107722, doi: https://doi.org/10.1016/j.cpc.2020.107722
  * Kopal (1961) Kopal, Z. 1961, Numerical Analysis (London: Chapman & Hall)
  * Lecoanet et al. (2016) Lecoanet, D., McCourt, M., Quataert, E., et al. 2016, MNRAS, 455, 4274, doi: 10.1093/mnras/stv2564
  * Lele (1992) Lele, S. 1992, J. Comp. Phys., 103, 16
  * Lerat et al. (2007) Lerat, A., Falissard, F., & Sidès, J. 2007, J. Chem. Phys., 225, 635, doi: https://doi.org/10.1016/j.jcp.2006.12.025
  * Lesur et al. (2022) Lesur, G., Ercolano, B., Flock, M., et al. 2022, in Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura (Astro. Soc. of the Pacific), 465–500.  https://arxiv.org/abs/2203.09821
  * Liepmann & Roshko (2001) Liepmann, H., & Roshko, A. 2001, Elements of Gasdynamics (Dover Publications)
  * Mani et al. (2009) Mani, A., Larsson, J., & Moin, P. 2009, J. Comp. Phys., 228, 7368, doi: doi.org/10.1051/aas:2000116
  * Marcus (1984) Marcus, P. S. 1984, J. Fluid Mech., 146, 65–113, doi: 10.1017/S0022112084001774
  * Masset (2000) Masset, F. 2000, ApJS, 141, 165, doi: 10.1051/aas:2000116
  * Merriman (2003) Merriman, B. 2003, J. Sci. Comp., 19, 309, doi: 10.1023/A:1025312210724
  * Mignone (2021) Mignone, A. 2021, Pluto Users Guide, Available at http://plutocode.ph.unito.it/userguide.pdf
  * Mignone et al. (2007) Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228, doi: 10.1086/513316
  * Moser et al. (1983) Moser, R., Moin, P., & Leonard, A. 1983, J. Chem. Phys., 52, 524, doi: https://doi.org/10.1016/0021-9991(83)90006-2
  * Nelson et al. (2013) Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610, doi: 10.1093/mnras/stt1475
  * Neufeld & Hollenbach (1994) Neufeld, D. A., & Hollenbach, D. J. 1994, ApJ, 428, 170, doi: 10.1086/174230
  * Pirozzoli (2002) Pirozzoli, S. 2002, J. Chem. Phys., 178, 81, doi: https://doi.org/10.1006/jcph.2002.7021
  * Ritos et al. (2018) Ritos, K., Kokkinakis, I., & Drikakis, D. 2018, Computers & Fluids, 173, 307, doi: https://doi.org/10.1016/j.compfluid.2018.01.030
  * Seligman & Laughlin (2017) Seligman, D., & Laughlin, G. 2017, ApJ, 848, 54, doi: 10.3847/1538-4357/aa8e45
  * Seligman & Shariff (2019) Seligman, D., & Shariff, K. 2019, ApJ, 877, 113, doi: 10.3847/1538-4357/ab13a7
  * Shariff (2024) Shariff, K. 2024, Pade: A code for protoplanetary disk turbulence based on Pade differencing., 1.0, Zenodo, doi: 10.5281/zenodo.11114378
  * Shariff & Wray (2018) Shariff, K., & Wray, A. 2018, ApJS, 238, 12, doi: 10.3847/1538-4365/aad907
  * Shu & Osher (1989) Shu, C.-W., & Osher, S. 1989, J. of Comp. Phys., 83, 32, doi: 10.1016/0021-9991(89)90222-2
  * Smagorinsky (1963) Smagorinsky, J. 1963, Mon. Weather Rev., 91
  * Stone et al. (2008) Stone, J., Gardiner, T., Teuben, P., Hawley, J., & Simon, J. 2008, ApJS, 178, 137, doi: 10.1086/588755
  * Stone et al. (2020) Stone, J., Tomida, K., White, C., & Felker, K. 2020, ApJS, 249, 4, doi: 10.3847/1538-4365/ab929b
  * Vreman (2004) Vreman, A. W. 2004, Phys. Fluids, 16, 3670, doi: 10.1063/1.1785131
  * Yee & Sjögreen (2018) Yee, H., & Sjögreen, B. 2018, Computers & Fluids, 169, 331, doi: https://doi.org/10.1016/j.compfluid.2017.08.028