{'fraction_non_alphanumeric': 0.035465955306773954, 'fraction_numerical': 0.015437092753749946, 'mean_word_length': 4.631695603989656, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Poloidal-toroidal decomposition in a finite cylinder: II. Discretization, regularization and validation

The Navier-Stokes equations in a finite cylinder are written in terms of poloidal and toroidal potentials in order to impose incompressibility. Regularity of the solutions is ensured in several ways: First, the potentials are represented using a spectral basis which is analytic at the cylindrical axis. Second, the non-physical discontinuous boundary conditions at the cylindrical corners are smoothed using a polynomial approximation to a steep exponential profile. Third, the nonlinear term is evaluated in such a way as to eliminate singularities. The resulting pseudo-spectral code is tested using exact polynomial solutions and the spectral convergence of the coefficients is demonstrated. Our solutions are shown to agree with exact polynomial solutions and with previous axisymmetric calculations of vortex breakdown and of nonaxisymmetric calculations of onset of helical spirals. Parallelization by azimuthal wavenumber is shown to be highly effective.


Introduction
The von Kármán flow owes its name to T. von Kármán [1] who in 1921 first studied the flow in the semi-infinite domain bounded by a single rotating disk using a similarity transformation. In 1951, Batchelor [2] extended the problem to the flow confined between two infinite rotating disks. For rotating disks of finite radius, the configuration can be described by three control parameters: the ratio s of the angular velocity of the two disks, the height-to-radius ratio h and a Reynolds number Re based on the radius and the azimuthal velocity of one of the disks. The variation of these three parameters (Re, s, h) has proved to yield a rich variety of qualitatively different accessible flows, even before the onset of turbulence. The symmetries influence the transitions that the flow can undergo. This configuration is extensively studied in the context of transition to complex and turbulent flows. All of these properties explain why the von Kármán flow is increasingly considered as one of the classical hydrodynamic configurations and why the scientific community is interested in further exploring its complex behavior.
The first numerical studies were necessarily devoted to axisymmetric flows and their stability. In the rotor-stator configuration (s = 0), vortex breakdown forming characteristic recirculation bubbles was observed by number of authors (Lugt and Abboud [3], Daube and Sorensen [4], Lopez [5], Daube [6]). This now well-documented configuration became a benchmark for testing axisymmetric codes. Following Lopez and Shen [7], Speetjens [8] and a number of other authors, we will validate our method in the axisymmetric configuration by reproducing the stationary state at Re = 1800 and, for Re = 2800, the non-stationary, oscillating flow, for which we will compare the bifurcation threshold and the oscillation frequency against previous findings.
Increasing computational power has made it possible to study three-dimensional instabilities. Breaking of axisymmetry has been the subject of several studies, of which we mention these of Gauthier et al. [9], Gelfgat et al. [10], Blackburn and Lopez [11], Lopez et al. [12] and Nore et al. [13]. Three-dimensional instability precedes axisymmetric instability for h < 1.6 and h > 2.8. As the test problem for validating our code in three dimensions, we have selected a configuration with (s = −1, h = 3.5, Re = 2150), where the perturbation of the axisymmetric flow takes the form of a helical spiral. In section 3.5 we compare our results with those of Lopez et al. [14] and Gelfgat et al. [10].
Interesting phenomena can also be observed in the turbulent regime. Turbulence and large-scale structures may coexist. In experiments in a turbulent counter-rotating configuration, Marié [15] and Ravelet et al. [16] discovered that a two-cell mean flow with a shear layer at the cylinder mid-plane undergoes switching to a one-cell mean flow whose shear layer is adjacent to the less rapidly rotating disk. This transition can be observed at Reynolds numbers which are numerically accessible.
A comprehensive classification of the solutions for different values of the parameters (s, h, Re) is beyond the scope of this work. Our main purpose here is to develop a mathematical and algorithmic tool which can be applied to von Kármán flow and to rotating turbulence, and which can be extended to the magnetohydrodynamic configuration of the VKS experiment [17]. The major component of our algorithm is the poloidal-toroidal decomposition [18,19], which insures incompressibility by construction, at the price of increasing the order of the governing equations. When applied to the Navier-Stokes equation in a finite cylinder, the resulting system has boundary conditions which are coupled and of high order. In a companion article [20], we showed that this system could be reduced to the solution of a set of nested Helmholtz and Poisson problems with uncoupled Dirichlet boundary conditions, whose solutions could be superposed via the influence matrix technique. The purpose of the present article is to describe our method for solving these elliptic problems using a spectral representation which exploits azimuthal symmetry of the system and which is regular at the cylindrical axis, and to demonstrate the validity of the resulting hydrodynamic code.
More specifically, it was shown by Marques and co-workers [18,19] that the Navier-Stokes equations in a finite cylinder with boundary conditions and with toroidal and poloidal potentials defined by u = ∇ × (ψê z ) + ∇ × ∇ × (φê z ) (1.3) are equivalent to the two scalar equations Our article [20] was devoted to showing how the problem (1.4)-(1.5) can in turn be reduced to a sequence of five nested parabolic and elliptic problems, each with Dirichlet boundary conditions: for the two potentials ψ, φ and three intermediate fields g, f ψ and f φ . The influence matrix technique [21], a generalization of the usual separation into particular and homogeneous solutions, is used to determine boundary values σ f , σ g , σ ± g such that the boundary conditions present in (1.5) but not in (1.6) are satisfied, i.e. such that In this article, we describe the numerical implementation of this algorithm. We first present the spatial discretization of the fields, using a set of basis functions [22] which is regular at the cylindrical axis r = 0, and regularizing the discontinuous boundary conditions at the corners r = 1, z = ±h/2. We then explain the methods we have used for solving equations (1.6), in particular for stably and economically solving the Helmholtz problems resulting from time discretization and for evaluating the nonlinear terms. Finally, we describe the validation of the implementation, comparing results from our code to an analytic polynomial solution and to previously published two-and threedimensional test cases.

Spatial and temporal discretization
The spectral discretization that we use is In (2.1), we do not introduce new notation for Fourier coefficients, or for coefficients in the 3D tensor-product basis, using instead the number and type of superscripts and subscripts to distinguish between functions in physical space and spectral space coefficients.
The basis functions in the azimuthal and axial directions are standard [23,24]: Fourier modes e imθ and Chebyshev polynomials T k (2z/h), respectively. In cylindrical geometries, it is the radial direction which is most problematic and on which we will focus. To represent this direction, we use the basis functions Q m n (r) developed by Matsushima and Marcus [22]. In sections 2.1 and 2.2, we will discuss the means by which we impose regularity at the origin r = 0 and at the corners r = 1.

Regular basis of radial polynomials
A function f on the disk is analytic at the origin if the radial dependence of the Fourier coefficient f m (r) multiplying the Fourier mode e imθ is of the form where p is a polynomial. Examples of functions which violate (2.2) are given in figures 1, 2 and 3 and compared with regular functions obeying (2.2) on the right of each figure.
Various approaches used in spectral methods to represent functions in polar coordinates are surveyed by Boyd [25,26] and by Canuto et al. [24]. A common practice [14,21,27,28] has been to impose some degree of continuity, such as C 3 , but not complete analyticity C ∞ . Although basis functions which are not analytic at the origin generally do not pollute the fields, retaining such functions wastes CPU time and memory which could be better used to increase resolution.
The condition (2.2) is stated in terms of monomials, the use of which is excluded because of their poor numerical properties. The polynomial basis developed by Matsushima and Marcus [22] respects these conditions and yet is numerically well-conditioned. These polynomials Q m n (α, β; r) are solutions to the singular Sturm-Liouville equation: 3) defined over r ∈ [0, 1]. In (2.3), 0 ≤ |m| ≤ n, α ∈ [0, 1] and β is a positive integer. With the special choice α = β = 1, Q m n (1, 1; r) are related to Legendre and shifted Jacobi polynomials used by Leonard and Wray [29]; similar functions were also derived by Verkley [30]. The functions Q m n (α, β; r) are complete and orthogonal over [0, 1] with respect to   . Left: f m (r, θ) = r 2 + 32r 2 cos(mθ), which is not smooth at 0 for m > 2. Right: f m (r, θ) = r 2 + 32r m cos(mθ), which is smooth at 0. the inner product: The n th order polynomials Q m n (α, β; r) have the following explicit expression: (2.5) but they, as well as the normalizing coefficients I m n (α, β), can be calculated in O(n − |m|) operations using recursion relations given by Matsushima and Marcus [22]. Recursion relations also exist for the operators r∂ r , r 2 , (r∂ r ) 2 − m 2 , (r∂ r ) 2 + λr 2 (2.6) expressed in the Q m n polynomial basis, meaning that for any of the operators H in (2.6), there exist banded matrices L and R such that H = R −1 L. Thus, reducing the time for multiplication by H or H −1 from quadratic to linear in the number of radial modes or gridpoints. The existence of recursion relations is a general property of differential operators represented in polynomial bases, for reasons explained by Tuckerman [31]. Recursion relations will be further discussed in section 2.4.
The radial function f m (r) associated with Fourier mode m and its coefficients f m n in the polynomial basis are related by the transform pair: if the boundary points are to be excluded; otherwise, they are the roots of a slightly more complicated expression [22]. Once the {r i } are determined, the weights {w i } are computed by recursion relations [22].
Equations (2.8) specify N m ≡ N − [m/2] coefficients from values at N quadrature points via a rectangular matrix. Since the basis is orthonormal, the inverse transformation is obtained from the transpose of this rectangular matrix. The smaller size of the spectral representation compared to the grid representation is a consequence of the fact that the functions in (2.8) are not arbitrary functions of r, but belong to the restricted subspace of functions obeying the regularity conditions (2.2).

Regularization of the corners
The boundary conditions on u θ stated in (1.2) are The equations have been nondimensionalized using the radius as the unit of length and the inverse angular velocity 1/ω − (with ω − > 0 and |ω + | ≤ ω − ) as the unit of time. Therefore, ω − = 1 and −1 ≤ s = ω + ≤ 1, in particular ω + is 0, +1 or −1, for the rotor-stator, exactly corotating, or exactly counter-rotating configurations, respectively. (It is also possible to set the velocity on the cylinder to some other constant value instead of 0; for example to simulate the flow in a rotating cylinder, we would set Here, we focus instead on the difficulties in implementing the discontinuous boundary conditions (2.9).) Boundary conditions (2.9) are discontinuous at the corner points r = 1, z = −h/2 and possibly r = 1, z = +h/2: one or both disks rotate while the lateral boundary remains fixed. Mathematically, a PDE with a finite number of singular points can have a solution which is smooth except at these points. However, spectral methods then do not converge exponentially because series of smooth functions cannot converge uniformly to a discontinuous solution. If nothing is done to prevent it, the Gibbs phenomenon will lead to spurious oscillations which propagate into the whole domain from the neighborhood of the singularity. For finite difference methods, the discontinuity will affect only the neighborhood of the singular point, on the order of the grid interval, and therefore does not pose a severe problem. Finite volume methods have a local integral formulation and so the discontinuity presents an even less serious problem. The filtering intrinsic to local methods is, however, intrinsically related to the high numerical diffusion which in turn makes local methods less precise.
In some cases, even local methods do not sufficiently filter singularities. Georgiou et al. [32] discuss the issue of spurious oscillations in the context of finite element methods.
In the solutocapillary problem studied by Martin-Witkowski and Walker [33], the authors were required to explicitly filter the solution to achieve acceptable convergence even in a finite difference calculation.
Spectral methods must always explicitly filter strong singularities like (2.9a)-(2.9b). We have chosen to do this by approximating the discontinuous function at the boundary by a steep but smooth profile. This procedure can be justified by arguing that we are not interested in finding the solution to the singular problem. In any real experiment, the boundary conditions are not discontinuous: a small gap must necessarily exist between the rotating disks and the stationary lateral boundary. In our algorithm, we replace the constant angular velocities in (2.9a) by continuous functions, as illustrated by figure 4a. Two possibilities are: where µ is an arbitrary but large even integer (e.g. µ = 10), and δ, ω ± are constants.
In order to be represented in a radial polynomial basis, the exponential regularization (2.10a) must be approximated by a polynomial, so that both expressions above are effectively polynomials. The steepness of the profiles are adjusted by varying δ or µ, whose possible values are limited by the radial polynomial order N.
The advantage of exponential regularization is that, for a given N, a steeper profile can be achieved by the polynomial approximation to (2.10a) than by the polynomial (2.10b), as can be seen in figure 5a. In this way, the deviation from the idealized profile is minimized, without unduly increasing the radial resolution N over that required to resolve the field in the interior. The polynomial approximation to (2.10a) differs from 1 by more than 10% only over the small range 1 − 2δ r < 1.
Regularization of the boundary condition imposes a lower bound on the spectral resolution -the spectral approximation must be able to represent the regularization profiles smoothly. In fact, Lopez and Shen [7] observed that the actual resolution should be approximately twice the minimal resolution sufficient for representing the regularization profiles, due to generation of higher wavenumbers by the nonlinear terms. As was shown by Lopez and Shen [7], for comparing with results obtained using different methods and for benchmark purposes it is sufficient to use δ ≈ 0.005. In practice, we use 0.005 < δ < 0.05, illustrated in figure 5b.
Care must be taken in approximating (2.10a) by a polynomial expression, in order to satisfy all of the conditions that we require of u θ and of ψ. The procedure we use is as follows. We evaluate (2.10a) on the collocation points. Since u θ is an odd function of r (see (2.30)), we apply the transform (2.8) using the (odd) polynomials associated with the m = 1 Fourier mode. The basis of odd radial functions insures that this ap-proximation to u θ on the bounding disks is zero at r = 0, while the use of Gauss-Radau collocation points which include the cylinder boundary insures that that it is also zero at r = 1. This approximation to u θ is integrated over r to obtain an approximation to ψ, with the integration constant chosen in order to satisfy ψ = 0 at r = 0.
Another choice is to apply a filter to the lateral boundary, replacing (2.9b) by while keeping the uniform angular velocity profiles (2.9a) on the disks unchanged; see figure 4b. This kind of regularization is similar to that of Lopez and Shen [7].
Other, quite different, approaches to the treatment of singularities exist. One is singularity subtraction. The form of the singular part of the solution is determined analytically, and the solution is written as a sum of the singular solution and an unknown regular part. Only the regular part is treated numerically. The effect of the singular solution on the numerical one can be filtered down to the scales representable by the spatial resolution. The main advantage of this method is that it recovers the convergence of the scheme and at the same time approaches the exact solution. Recent applications of this method are to the driven cavity problem [34] and to injection of fluid into a cylindrical channel [35]. The results obtained are generally of high precision and often provide a benchmark for a particular problem. The main drawback is that it requires knowledge of the solution near the singular point. For the 2D driven cavity problem, the nature of the singularity was given by Dean and Montagnon [36] and Moffatt [37] for a Stokes flow. For most inertial (Navier-Stokes) flows, and for 3D flows Hills and Moffatt [38], the analytic form of the singular solution is unknown. Note that even when the velocity boundary conditions are continuous, lower-order singularities of purely geometric origin are present at the corners. We will return to this in section 3.2.
Another approach is to derive a physically justified model which is no longer singular. Methods from molecular dynamics reflect the microscopic nature of the fluid at the smallest scales but are very hard to adapt to problems containing both large and small scales. Several continuous (macroscopic) approaches have been proposed as a compromise between a continuous and a molecular description. These all introduce a spatially limited physical effect which effectively removes the singularity. In this category are methods based on variable slippage, as well as the surface viscosity or dynamic surface tension applicable to free-surface problems. A comprehensive review of physically justified models, as well as other regularization techniques, is provided by Nguyen and co-workers [39,40].

Time integration
We use an implicit scheme for the linear diffusive terms while treating the nonlinear terms explicitly. Spectral methods require that the coefficients representing the solution decay with their index or wavenumber. The nonlinear term in the Navier-Stokes equation can be seen as a generator and amplifier of high wavenumbers while the viscous term damps these high wavenumbers. The intensity of this damping depends on the particular time-integration scheme and on the way the Laplacian is evaluated and must be strong enough to oppose the effect of the nonlinear term. In our case, high-wavenumber modes are needed to represent both the thin boundary layers created near the rotating disks and the steep regularized boundary profile. Fortunately, these effects are most pronounced in the proximity of the boundaries, where the ax-ial Chebyshev and radial polynomial grid is finest. However, in the counter-rotating case, the central shear layer can also require high wavenumber modes in order to be well-represented.
We use the first-order backward Euler scheme for linear terms because it attenuates high frequencies faster than all other methods. Tests performed with the Crank-Nicolson method confirmed that for this scheme, nonlinear simulation was unstable even for quite small Reynolds numbers Re ≈ 300. This behavior was also observed for the von Kármán flow by Speetjens [8] and Lopez et al. [14] and by Marcus [41] in Taylor-Couette flow. The nonlinear term is treated by the second-order explicit Adams-Bashforth scheme, so that the backward Euler/Adams-Bashforth time-integration scheme for the potentials ψ and φ takes the following form: where S ψ and S p si are defined in (1.4). Equations (2.12) can be written as the nested system of equations: As explained in section 1, the boundary conditions imposed on (2.13) are Dirichlet conditions with boundary values calculated via the influence matrix in such a way as to satisfy the more complicated coupled boundary conditions given in (1.5).
The maximal time step δt depends on the Reynolds number. Typically starting from state u = 0 requires a δt which is 4-10 times smaller than that which can be used for evolving a fully developed state at the same Reynolds number. This is because the state u = 0 is incompatible with the boundary condition (2.9a). In the first few iterations a boundary layer is created near the rotating cylinder lids, requiring higher spatial resolution. This can be avoided by performing about 100 initial steps of the linear Stokes solver, i.e. with S ψ = S φ = 0. In

Viscous terms
We now describe the way in which the Helmholtz and Poisson problems with Dirichlet boundary conditions in (1.6) or (2.13) are solved. The azimuthal Fourier representation in (2.1) decomposes each 3D elliptic problem into a set of 2D problems, each of which is associated to a single azimuthal Fourier mode m. The reflection symmetry in z leads to further decoupling between the set of modes that are symmetric or antisymmetric in z. Each of the 2M resulting elliptic problems corresponds to a single azimuthal Fourier mode m and axial parity p ∈ {s, a}, within which each equation corresponds to a value of (k, n), the indices of the axial and radial basis functions. The number of axial modes of each parity is K/2, and the number the number of radial basis functions corresponding to Fourier mode m is N m ≡ N − [ m 2 ]; in the remainder of this section, we will take m = 0, so that N m = N.
In the full cylinder, regularity at the axis serves as one of the boundary conditions and is imposed by the use of the Q n radial basis in (2.1), leaving only the boundary condition at r = 1 to be imposed. In the axial direction, the boundary conditions at the two disks can be recombined to yield one condition for each parity. Thus, one radial and one axial boundary condition remains to be imposed on each 2D problem. These are imposed via the τ method [23,24], so that the equations corresponding to the highestwavenumber modes in each direction are replaced by the boundary conditions. Gaussian elimination is especially economical for systems resulting from the spectral discretization of differential equations, whose structure is barely altered when boundary conditions are imposed via the τ method, since recursion relations reduce the solution time from quadratic to linear in the number of modes if the resulting systems are diagonally dominant. This can only be done in one direction however. In geometries with more than one non-periodic direction, the remaining directions must be treated by diagonalization. Here, we treat the axial direction by incorporating the boundary condition via Schur decomposition into the matrices representing ∂ 2 z for each parity and diagonalizing [24,42], leading to decoupled problems for each axial eigenvalue λ z . The operation count at each timestep, dominated by multiplication by the eigenvector matrix, is quadratic in K/2 and linear in N.
Thus, the 2D and 3D problems in (2.13) are all decomposed into a set of one-dimensional radial problems: which we will write in practice as where m is the Fourier mode. The scalar λ is 0 or λ z for the Poisson problems (2.13b), (2.13d) or (2.13e), or λ z + Re/δt for the Helmholtz problems (2.13a) and (2.13c) (with the multiplicative factor −Re/δt incorporated into g). For nonzero λ, L is pentadiagonal and R is tridiagonal. Two obstacles must be surmounted before (2.16) can be solved in a time which is linear in the number of radial modes. First, L is not diagonally dominant, so that stable Gaussian elimination would require pivoting, destroying the banded structure. Second the radial boundary condition must be included. A method which overcomes these difficulties was presented by Matsushima and Marcus [22]. Here, we will recast this method in terms of the Sherman-Morrison-Woodbury formula, which can be shown [21] to underly a large class of transformations between coupled and uncoupled systems.
System (2.16) must be altered in two ways. First, since the largest element of L is located on the first super-digonal, permuting its rows leads to a matrix PL which is diagonally dominant, but is no longer banded. Second, the boundary conditions must be imposed. The tau method replaces (2.15) by which effectively discards g N by adjusting it with the extra unknown τ, introduced along with the boundary condition (2.17b). B T is the row vector which represents the discretized boundary condition andê N is the unit vector which selects the N th component. Multiplying by r 2 , using (2.16) and permuting rows leads to: This is rewritten in matrix form as where the subscript lo refers to all indices lower than N. System (2.19) is solved by using the Sherman-Morrison-Woodbury formula which relates the inverses of two matrices differing by a low-rank transformation vw T , in particular differing by a few rows or columns. We define A, v and w T for (2.19) as follows. The matrix in (2.19) is replaced by another matrix A, which is more easily inverted since it is block upper triangular: where a is an arbitrary value whose order of magnitude is that of the dominant values of L andê T 1 is a unit row vector corresponding to the lowest radial wavenumber present for this m. The upper-left matrix in (2.21) is banded and diagonally dominant and so can be stably inverted without any pivoting. The matrices in (2.21) and (2.19) differ only in their first and last rows, so their difference vw T is of rank two: We define: In addition to the ability to solve (2.21), the Sherman-Morrison-Woodbury formula (2.20) requires only the inversion of the following 2 × 2 matrix: where we have used the fact that aê T 1 (L ′ ) −1ê 1 = 1/a.

Nonlinear terms
To compute the nonlinear terms we use the pseudo-spectral method [23], in which fields are transformed into physical space, the nonlinear terms are carried out via pointwise multiplication, and the results transformed back into spectral space. Computing the nonlinear term S in the rotational form −u × ω requires only 9 spectral↔physical transforms as compared to the 15 transforms required by the convective form (u · ∇)u. The difference between them is annihilated by the curls taken in (2.25), so we will write S = −u × ω.
The calculation of the nonlinear terms presents two difficulties. The first involves radial parity and appears when creating S. We have sought to use only scalar fields which can be represented by expansions of type (2.1) and constructed using radial operators such as (2.6) which preserve radial parity and can be implemented via recursion relations. The components of velocity and vorticity, defined in cylindrical coordinates using the toroidal and poroidal potentials as do not have this property. We therefore construct modified fields: whose components have the same parity as ψ and φ and so can be created and acted upon using the differential operators in (2.6). The modified fields u * and ω * are transformed into physical space, where their cross product is taken to form The second difficulty appears when differentiating S in (2.25) and involves regularity. A vector function which is analytic at the origin must obey conditions analogous to (2.2), namely where p r , p θ and p z are polynomials. We require not only regularity at the origin of S, but also regularity of its curl and double curl. When S = −u × ω is calculated analytically, this is in fact the case. However, the numerical transforms to and from physical space introduce aliasing errors which destroy this property. Full dealiasing would multiply the time necessary for evaluating of the nonlinear term by a factor ≈ 4.5. Matsushima and Marcus [22] suggest instead that all terms that could potentially suffer in spectral space from singular operations (like dividing by r) be evaluated in physical space (at collocation points excluding the coordinate origin) and transformed back to the spectral space using the radial transform, ensuring the correct polynomial order for a given Fourier mode. We have generalized this approach to the evaluation of S ψ and S φ . For each Fourier mode m, we write The purpose of the decomposition in (2.31) can be explained as follows. Consider the first terms on the right-most-sides of (2.31), those which do not contain (S * r + iS * θ ). Because they do not generate terms of lower polynomial order, these terms preserve regularity and can be carried out in spectral space. In contrast, although (S * r + iS * θ ) should be divisible by r 2 , aliasing error in the multiplication (2.26) can generate terms of lower polynomial order. (More specifically, S * and iS * θ are not individually divisible by r 2 and aliasing perturbs the cancellation in the sum S * + iS * θ .) This division by r 2 is therefore carried out in physical space. This increases the number of spectral↔physical transformations only from 9 to 10. Further details can be found in [43].

Parallelization
The separability of almost the entire algorithm (except for the nonlinear term) between the Fourier modes makes parallelization of the code quite straightforward. Our code was parallelized using the MPI protocol which made it possible to run even very time-consuming three-dimensional simulations with resolutions such as M × K × N = 128 × 160 × 90. Spectral methods are often considered to be poorly suited for parallelization as they require the exchange of all the data at each timestep of the simulation. In our code, all the necessary data exchange is done within two calls to the MPI_Alltoall MPI subroutine treating, in total, 10 three-dimensional fields at each time step. Even though this may seem to be a large operation, on the IBM Power4 architecture with 64 processors we found that the time overhead per timestep due to the data exchange is counterbalanced by more efficient usage of the processor cache memory: each processor of the parallelized code treats smaller data portions which can more easily fit into processor's fast internal memory (cache). We observed that the total CPU time used by the parallel code is often smaller than that used by the serial code treating the same problem, as shown in table 2. The efficiency of the parallel code depends, however, on the speed and the latency of the inter-processor network: the IBM Power4 architecture has particularly fast connection between the nodes which use the mixed model fast network/shared memory communication between processors. We conclude that for modern massively parallel computers, parallelization of a pseudo-spectral code does not necessarily degrade its efficiency but can actually enhance it. Additional technical information about the MPI parallelization of the code is given in [43].  Table 2 CPU timings on IBM Power4, as the azimuthal resolution M and number of processors is varied.

Tests and validation
We now describe the ways in which we have validated the hydrodynamic code described here and in our companion paper [20]. We have obtained an exact polynomial solution to the nested Helmholtz-Poisson solver, which is the by far most complicated portion of the code; we present its form in the hopes it may prove useful to other researchers. For non-polynomial solutions, we have analyzed the effect of the corner singularity on the error and verified the spectral convergence. We have tested the full nonlinear time-dependent program by simulating 2D and 3D stationary and time-periodic rotor-stator flows which are well documented in the literature.

Polynomial solutions
There are no polynomial solutions which exactly satisfy the linear differential equations (1.4) and boundary conditions (1.5) with which to compare a numerical solution. We can, however, formulate a polynomial solution to a related problem: the timediscretized equations relating the four fields ψ, φ, rhs ψ , rhs φ and subject to boundary conditions (1.5). Under these conditions, no error is introduced in imposing the boundary conditions via the τ method (see section 2.4). Comparison between this analytic solution and a numerical solution provides a test of the linear Helmholtz/Poisson solver, including the enforcement of the boundary conditions via the influence matrix method, which is independent of the evaluation of the nonlinear terms, the temporal integration and the spatial resolution.
Using a symbolic algorithm implemented in Maple, we have calculated polynomial solutions which contain several radial and axial basis functions (wavenumbers) and which correspond to a realistic profile ω ± , such as ω poly ± = ±(1 − r 6 ). Polynomial solutions were calculated for Fourier modes ranging from m = 0 to m = 5. For m = 0, we sought a solution containing two recirculation rolls separated by the mid-plane z = 0, leading to the potentials, right-hand-sides, and velocities: For all polynomial solutions tested, we obtained numerical solutions for which the maximum relative errors in ψ, φ and in satisfaction of the boundary conditions, as well as the equations, are typically O 10 −14 , i.e. machine precision, and never exceed O 10 −12 .

Non-polynomial solutions: error analysis
If the right-hand-sides rhs ψ and rhs φ are not polynomials that can be exactly represented within the spatial discretization, then the equations cannot be exactly satisfied. Note that the tau error introduced in favor of satisfying the boundary conditions at each level of the nested system of equations is necessarily propagated to the next level. That is, ψ n+1 is not an exact solution to the Poisson problem (2.13b) and in addition, the right-hand-side f ψ is not an exact solution to the Helmholtz problem (2.13a). This implies that the error is not isolated in the equations corresponding to the highest wavenumbers, but is distributed among all the equations. However, we will see below that non-satisfaction of the equations for low wavenumbers does not have severe consequences.
In figures 6 and 7, we display the relative error in physical space (r, z) for the m = 0 and m = 2 modes. In this subsection, as in the previous one, we restrict ourselves to the Stokes problem, whose errors behave like those of the Navier-Stokes equations. As is usually the case, the error is concentrated in the neighborhood of the boundaries, decaying rapidly away from them. These errors represent quite severe criteria, since these measure the satisfaction of the curl and double curl of the original equations; the error in satisfying the Stokes equations themselves is considerably lower, on the order of ǫ ψ divided by O(K, N) and ǫ φ divided by O(K 2 , N 2 ) near the boundary. This estimate relies on the fact that the spec-trum of the error in satisfying the Stokes equations is approximately uniform, so that its derivatives are dominated by its high wavenumber terms.
The error depends significantly on the right-hand-sides rhs ψ and rhs φ . For arbitrary right-hand-sides, the boundary conditions can be very constraining and can lead to a nearly singular solution suffering from spurious oscillations, with an error near the boundary that is O(1). However, ǫ ψ , ǫ φ are considerably smaller when the right-handsides are calculated from the solution at the previous timestep, especially for the Stokes problem for which the nonlinear term is zero.
For the axisymmetric modes shown in figure 6, ǫ ψ is O(10 −6 ) on the boundaries and O(10 −10 ) for the internal points. This is much less than ǫ φ , which reaches O(1) at the cylinder corners. To understand this, we recall that in the axisymmetric case, the toroidal flow described by ψ is azimuthal and the poloidal flow described by φ is in the (r, z) plane. The azimuthal flow described by ψ follows smooth paths, while the flow described by φ must abruptly change direction near the corners. This poloidal flow in fact resembles the analytic asymptotic solution derived by Moffatt [37] for 2D flow in a rectangular container, which is weakly singular in that its vorticity behaves like ρ 1.74 (for a small distance ρ from the corner). This in turn implies that the axial component of the Laplacian of the Stokes equation measured by ǫ φ diverges for exact solutions to the continuous 2D Stokes problem, in contrast with the numerically computed solution which is forced to be finite and regular. Satisfaction of the Stokes equation itself follows from the satisfaction of the boundary condition which our code imposes to precision O (10 −14 ).
For the non-axisymmetric modes shown in figure 7, typical errors at the corners are O(0.01) for ψ and O(0.1) for φ and the same corner singularity is observed for both. This is because the 3D Stokes solutions are also weakly singular at the corners Hills and Moffatt [38] and ψ and φ are coupled for m = 0. The relative error decreases rapidly away from the boundaries: it is O(10 −4 ) only two gridpoints away and reaches reaches O(10 −6 ) for the interior points.
The accuracy could be further enhanced -or confined to certain modes -by completing the method with the τ-correction [21], which takes into account the high-wavenumber residuals resulting from the imposition of the boundary conditions for each Helmholtz or Poisson problem.

Spectral convergence
The important indicator of spatial convergence for spectral methods is the decay rate of high-wavenumber coefficients in the solution fields. For a well-behaved solver, in the absence of volume and boundary singularities, the magnitude of spectral coefficients should decay rapidly with wavenumber. For a laminar flow this decay rate should be exponential, but the presence of thin boundary layers can significantly influence the convergence. Figure 8 shows the spectral convergence for a full Navier-Stokes simulation at Re = 750 at T = 100δt = 1 from initial conditions of u θ = ω + (r)2z/h, u z = u r = 0. We show the (r, z) spectral coefficients for Fourier mode m = 0; those for m = 0 are similar. The convergence of the spectra can be qualified as quasi-exponential, meaning that the high-wavenumber spectral coefficients seem not to decrease below a level O(10 −12 ). This can almost certainly be attributed to the singular character of the solution to the Stokes equation near the cylinder corners. In support of this explanation, we note that the spectrum obtained after 100 timesteps of the linear (Stokes) solver behaves similarly, except for ψ m=0 , which displays true explonential decay, down to levels of 10 −22 for the same resolution.

Axisymmetric rotor-stator configuration
The code was first tested on the well-documented axisymmetric rotor-stator configuration with aspect ratio h = 2. The first test is the reproduction of the characteristic steady state for Re ≈ 1850 where the flow exhibits two recirculation bubbles (one large and the other much smaller) situated approximately at (r = 0, z = 1/2) and (r = 0, z = 0). The contour plot of the Stokes streamfunctions defined as presented on figure 10 matches that presented by Daube [6] (figure 9) and is similar to those of Lopez and Shen [7] obtained for the slightly larger aspect ratio h = 2.5. Quantitative agreement between our results and the previous calculations by Daube Fig. 9. Contours of Stokes poloidal streamfunction η from [6] for Re = 1850, h = 2. [6] and Lugt and Abboud [3] is established by comparing the profiles of axial velocity u z on the cylinder axis (see figure 12). This test shows excellent agreement between our results obtained using the poloidal-toroidal formulation (fig 12b) and the velocityvorticity formulation ( fig. 12a).
It was observed experimentally by Escudier [44] and numerically by Daube and Sorensen [4], Lopez [5], Daube [6], Gelfgat et al. [10] and Speetjens [8] that this flow undergoes a Hopf bifurcation toward a flow oscillating with a frequency approximately 0.25 times the rotation frequency. This transition occurs at Re near 2600 with a period of T = 26.55. These values are in the ranges previously found for this configuration (see table 3); deviations can probably be attributed to the differences in the treatment of the boundary conditions (regularization). The simulation was performed with δt = 0.01 and with high spatial resolution K × N = 140 × 70 in order to well represent the sharp regularization profile corresponding to δ = 0.06 imposed on the lateral boundary r = 1 (see section 2.2) as proposed by Lopez and Shen [7] and also used by Speetjens [8]. The time evolution u θ (r = 0.5, z = 0, t), along with the normalized power spectrum, are shown on figures 13 and 14.

First 3D instability
We have tested the non-axisymmetric aspects of the code using another rotor-stator configuration. For the aspect ratio of h = 3.5 we found that the first bifurcating mode has wavenumber m = 3 and critical Reynolds number Re c = 2116. This result is in good agreement with Gelfgat et al. [10], where the critical Reynolds number was estimated at Re = 2131. A characteristic spiral analogous to that visualized by Lopez et al. [14] 0

Conclusion
Motivated by a need for a numerical tool adequate for investigating cylindrical von Kármán flow, we have written a spectral code which solves the Navier-Stokes equation for an incompressible fluid in a finite cylinder. This task encompasses a number of algorithmic challenges, which we list in order of decreasing difficulty. The first is to impose incompressibility, a goal which we achieved by using the poloidal-toroidal decomposition. In a finite cylinder, this formulation leads to differential equations of higher order and coupled boundary conditions [18,19]. In a companion paper [20], we described the way in which we used the influence matrix technique to transform these equations and boundary conditions into an equivalent set of decoupled Helmholtz or Poisson problems, each with Dirichlet boundary conditions. The second challenge, which has been the main focus of this article, is the treatment of the cylindrical axis. The singularity engendered by the use of cylindrical coordinates is only apparent and should not be transmitted to the fields. We have dealt with the axis singularity by using a basis of radial polynomials developed by Matsushima and Marcus [22] which are analytic at the axis and which have properties similar to those of the Legendre polynomials. We have extended [22] in several ways. (i) The radial basis was developed for two dimensions, i.e.polar coordinates. Here, we have used it to represent fields in a three-dimensional cylinder of finite (non-periodic) axial length. This fairly straightforward extension was implemented by diagonalizing the differential operators in the axial direction. (ii) Numerical inaccuracy in the pseudo-spectral evaluation of the nonlinear term can generate terms which are not analytic. We have generalized to the high-order equations resulting from the poloidal-toroidal decomposition the procedure developed by Matsushima and Marcus [22] to avoid this problem. Care must be taken to preserve order and parity at each stage of the calculation. (iii) Differential operators expressed in this basis can, as in most such cases [31], be reformulated as recursion relations, which can be used to reduce the time for action or inversion of differential operators. Using the Sherman-Morrison-Woodbury formula, we have formalized the stable algorithm given in [22] for solving Poisson and Helmholtz problems in a time proportional to the number or radial gridpoints or modes.
The third challenge is the genuine singularity at the corners of the domain, where the disks and the cylinder which bound the domain meet. Finite difference and finite element codes intrinsically smooth the singularity; in contrast, spectral expansions attempt to converge to the discontinuity, leading to spurious oscillations. In our code, we replaced the discontinuous boundary conditions by a profile on the disks which is steep but continuous as the corner is approached. The geometrical singularity remains, but it is weak and does not prevent spectral convergence.
Tests performed for analytic polynomial solutions to the Helmholtz problem with an appropriate right-hand side showed that the solver reproduces exact solutions to nearly machine precision. For a non-polynomial solution, the solver displays exponential convergence of spectral coefficients of the solution. The potential equations (corresponding to the curl and double curl of the Navier-Stokes equations) are satisfied, with an error of only O(10 −10 ) for the interior points. The numerical code was parallelized using the MPI protocol. This made it possible to simulate nearly turbulent flow for Re = 5000 with a spatial resolution of 128 × 160 × 90.
Finally, we have validated the hydrodynamic code by testing it against well-documented problems in the literature, demonstrating the feasibility of calculating solutions to the time-dependent Navier-Stokes equations in a finite cylindrical geometry which are both analytic and divergence-free.