Openpipeflow 1.02b
Openpipeflow 1.02b
Openpipeflow 1.02b
, u
z
Field vectors.
If K ,=K or M ,=M, harmonics are truncated or
zeros appended.
1.4 Output
1.4.1 Snapshot data
Data saved every i_save_rate1 timesteps:
state????.cdf.dat
vel_spectrum.dat
All output is sent to the current directory, and state les are numbered 0000, 0001, 0002,. . . .
Each le can be renamed state.cdf.in should a restart be necessary. To list times t for
Chapter 1. Getting started 5
each saved state le,
grep state OUT
The spectrum les are overwritten each save as they are retrievable from the state data.
To verify sucient truncation, a quick prole of the energy spectrum can be plotted with
gnuplot> set log
gnuplot> plot vel_spectrum.dat w lp
1.4.2 Time-series data
Data saved every i_save_rate2 timesteps:
tim_step.dat t, t, t
f
, t
CFL
current and limiting step sizes
vel_energy.dat t, E energies and viscous dissipation
vel_friction.dat t, U
b
or , u
z
(r = 0)
z
, u
(x
0
) +
(x
k
x
0
)
2
2!
f
(x
0
) +. . .
.
.
.
f(x
k
) = f(x
0
) + (x
k
x
0
) f
(x
0
) +
(x
k
x
0
)
2
2!
f
(x
0
) +. . .
(2.1)
The expansions (2.1) can be written
f = Adf, df = A
1
f, f =
_
_
f(x
k
)
.
.
.
f(x
k
)
_
_, df =
_
_
f(x
0
)
f
(x
0
)
f
(x
0
)
.
.
.
_
_
(2.2)
Derivatives are calculated using weights from the appropriate row of A
1
.
2.1.2 Integration
Integrating (2.1) the indenite integral may be approximated about x
0
_
f(x) dx =
(x x
0
) f(x
0
) +
(xx
0
)
2
2!
f
(x
0
) +
(xx
0
)
3
3!
f
(x
0
) +. . .
(2.3)
= [
(x x
0
)
(xx
0
)
2
2!
. . .
] A
1
f (2.4)
The total integral is approximated by
1
2
n
_
x
n+1
x
n1
f(x) dx (2.5)
Chapter 2. Spatial representation 9
2.1.3 Matrix representation of linear operators
Linear equations can be written in the matrix form
L p = Mq +s. (2.6)
L and M are N N matrices representing linear operators. If q and s are known then the
right-hand side is simple to evaluate. Consider the equation
L p = q, (2.7)
where p is unknown. If the matrix operators are formed by linear combinations of the
weights in 2.1.1, using only k neighbouring points to each side, they are banded. Boundary
conditions are formed in the same way and are placed in the end rows. The equation is
solved by forward-backward substitution, using the banded LU-factorisation of L.
BC
p q L
=
LU
2.2 Fourier representaion , z
2.2.1 Expansion of variables
A(, z) =
km
A
km
e
i(kz+m
0
m)
(2.8)
A real implies A
km
= A
k,m
(coecients are conjugate-symmetric). It is sucient to keep
only m 0, and for m = 0 keep k 0.
2.2.2 Orthogonality
Exponentials
_
2
0
e
i(m+n)
d = 2
m,n
,
_ 2
0
e
i(k+j)
dz =
2
k,j
(2.9)
2.2.3 Energy integral
E =
1
2
_
A AdV =
2
2
km
_
[A
km
[
2
r dr (2.10)
E =
m=0
E
m
, E
m
=
_
_
2
2
_
A
2
00
r dr +
4
2
k>0
_
[A
2
k0
[ r dr, m = 0,
4
2
k
_
[A
km
[
2
r dr, m > 0.
(2.11)
Chapter 2. Spatial representation 10
2.2.4 Volume integral
E =
1
2
_
AdV =
4
2
_
A
00
r dr (2.12)
2.3 Enforced conditions at the axis
The geometry enforces conditions on each variable at the axis when expanded over Fourier
modes in . Given a vector A, each mode for A
z
is even in r if m is even, and is odd if m
is odd. Each mode for A
r
and A
= u
R/ = (2 Re
m
(1 +))
1
2
= (2 Re)
1
2
3.1.3 Evolution equations
Fixed ux,
(
t
+u )u = p +
4
Re
m
(1 +) z +
1
Re
m
2
u (3.1)
Fixed pressure
(
t
+u )u = p +
4
Re
z +
1
Re
2
u (3.2)
Let u = W(r) z +u
2
) u
= u
( u
)
dW
dr
u
z
z W
z
u
+
4
Re
m
z p
. (3.3)
3.2 Decoupling the equations
The equations for u
r
and u
= u
r
i u
, (3.4)
Chapter 3. Model 12
for which the are considered respectively. Original variables are easily recovered
u
r
=
1
2
(u
+
+u
), u
=
i
2
(u
+
u
). (3.5)
Governing equations are then
(
t
2
) u
= N
(p)
, (3.6)
(
t
2
) u
z
= N
z
(p)
z
, (3.7)
where
=
2
1
r
2
i
r
2
(3.8)
3.3 Simple projection method
The basic projection method oers simplicity and exibility. The error measured by the
predictor-corrector is second order (on the nonlinear terms), but the method is rst order
overall. See following section for a better method.
The following steps are performed, where q refers to time t
q
:
u
u
q
t
= N
q+
1
2
, (3.9)
u
q+1
= u
, (3.10)
where
2
= u
. (3.11)
3.4 PPE-formulation with correct boundary conditions
Write the time-discretised NavierStokes equations in the form
_
Xu
q+1
= Yu
q
+N
q+
1
2
p ,
2
p = (Yu
q
+N
q+
1
2
),
(3.12)
where q denotes time t
q
, which is sixth order in r for u
q+1
and second order for p, where the
solenoidal condition is not explicitly imposed. Symmetry conditions provide the conditions
at the axis. The diculty is in imposing the remaining four this system should be
inverted, in principle, simultaneously for p and u
q+1
with boundary conditions u
q+1
= 0
and u
q+1
= 0 on r = R (Rempfer 2006). In practice it would be preferable to invert for
p rst then for u
q+1
, but the boundary conditions to not involve p directly.
Note that the Yu
q
term has been included in the right-hand side of the pressure-Poisson
equation, the divergence of which should be small. Assume that pressure boundary condi-
tion is known: the right-hand side of the NavierStokes equation is then projected onto the
space of solenoidal functions though p and hence after inversion, u
q+1
will be solenoidal.
Consider the bulk solution, u, p, obtained from solution of the following:
_
X u = Yu
q
+N
q+
1
2
p ,
2
p = (Yu
q
+N
q+
1
2
),
(3.13)
Chapter 3. Model 13
with boundary conditions u = 0 and
r
p = 0. Introduce the following systems:
_
Xu
= p
2
p
= 0,
(3.14)
with boundary conditions u = 0 and
r
p = 1 on r = R, and
_
X u = 0, (3.15)
with boundary conditions u
+
= 1, u
= 1, u
z
= i on r = R. The system (3.14) provides
a linearly independent function u
j
that may be added to u without aecting the right-
hand side in (3.13), but altering (to correct) the boundary condition applied. Similarly the
system (3.15) provides a further three functions. The superposition
u
q+1
= u +
4
j=1
a
j
u
j
(3.16)
may be formed in order to satisfy the four boundary conditions, u
q+1
= 0 and u
q+1
= 0
on r = R. Substituting (3.16) into the boundary conditions, they may be written
Aa = g( u), (3.17)
where A = A(g(u
j
(r), the matrix A and its inverse may all be precomputed. The boundary
conditions for u
z
is pure imaginary, and
A is real. For each timestep, this application of the inuence matrix technique requires
only evaluation of the deviation from the boundary condition, multiplication by a 44 real
matrix, and the addition of only two functions to each component of u, each either pure
real or pure imaginary. Compared to the evaluation of nonlinear terms, the computational
overhead is negligible.
3.5 Predictor-corrector
The model equation for each Fourier mode is
(
t
2
)f = N, (3.18)
where nonlinear terms have been evaluated on each radial point by the transform method,
and is solved as in 2.1.3. The predictor at time t
q
, with Euler nonlinear terms and implic-
itness c
f
q+1
1
f
q
t
_
c
2
f
q+1
1
+ (1 c)
2
f
q
_
= N
q
, (3.19)
_
1
t
c
2
_
f
q+1
1
=
_
1
t
+ (1 c)
2
_
f
q
+N
q
. (3.20)
Corrector iterations are
_
1
t
c
2
_
f
q+1
j+1
=
_
1
t
+ (1 c)
2
_
f
q
+c N
q+1
j
+ (1 c) N
q
, (3.21)
Chapter 3. Model 14
or equivalently, for the correction f
corr
= f
q+1
j+1
f
q+1
j
_
1
t
c
2
_
f
corr
= c N
q+1
j
c N
q+1
j1
, (3.22)
where j = 1, 2, . . . and N
q+1
0
= N
q
. The size of the correction |f
corr
| must reduce at each
iteration. For c =
1
2
the scheme is second order such that |f
corr
| t
2
.
3.6 Timestep control
t = C min(/ [v[), 0 < C < 1, (3.23)
where C is the Courant number.
The timestep should also be small enough such that the corrector norm |f
corr
| is satis-
factorily small. This may be important for integrating initial transients, but ought not to
be the limiting factor in general.
4
THE CODE
4.1 Naming conventions
4.1.1 Parameters
Parameters dened within parameters.f90 are given a prex indicating the data type.
For example:
N i_N integer
t d_timestep double
xed ux? b_fixed_flux boolean
Note that the naming convention frees the variable name for dummy variables, e.g.:
integer :: n
do n = 1, i_N
...
4.1.2 Public variables
By default variables declared within a module are accessible to any other module or
subroutine that uses the module. To ensure that it is clear where a public variable is
declared, variables are given a prex. Examples:
t tim_dt timestep.f90
u
r
vel_r velocity.f90
4.2 Ordering the Fourier modes
A =
km
A
km
e
i(kz+m
0
m)
, (4.1)
Coecients are stored as in the following example for m
0
=Mp= 2, M= 4 and K= 4:
Chapter 4. The code 16
m m_ nh
3 6 18 19 20 21 22 23 24
2 4 11 12 13 14 15 16 17
1 2 4 5 6 7 8 9 10
0 0 0 1 2 3
-3 -2 -1 0 1 2 3 k
Rather than two indices k and m, the single index nh labels the modes. In the code this
ordering corresponds to the typical loop
nh = -1
do m = 0, M-1
m_ = m*Mp
do k = -(K-1), (K-1)
if(m==0 .and. k<0) cycle
nh = nh+1
...
NOTE: this loop, and its parallel equivalent, are predended macros in parallel.h.
4.3 Data types
The data types below consist of logical data groups to facilitate the passing of data
between functions. For clarity, in the following text they are dened as though we were
working on a single processor. See 4.4 regarding subtle dierences in the denitions due
to parallelisation.
4.3.1 coll, spec storage of coecients
The collocated data type (coll) is the principle type used for mode-independent oper-
ations.
type coll
double precision :: Re(i_N, 0:i_H1)
double precision :: Im(i_N, 0:i_H1)
end type coll
The value Re(n,nh) is the real part of the coecient for the nh
th
harmonic at r
n
. H1= 24
for the example in 4.2 and N is the number of radial points.
The spectral type (spec) is the data type used for operations independent at each radial
point. It is principally a transitory data format between the coll and phys types.
type spec
double precision :: Re(0:i_H1, i_N)
double precision :: Im(0:i_H1, i_N)
end type spec
Conversions between the coll and spec types involve only a transpose, var_coll2spec()
and var_spec2coll().
Chapter 4. The code 17
4.3.2 phys real space data
The type (phys) is used for data in real space.
type phys
double precision :: Re(0:i_Z-1, 0:i_Th-1, i_N)
end type phys
The element Re(k,m,n) refers to the value at z
k
,
m
and r
n
.
4.3.3 mesh, lumesh storage of banded matrices
The nite-dierence stencil has nite width (2.1) and so matrix operations involve ma-
trices that are banded. The type (mesh), dened in meshs.f90, is used for matrices on
the radial mesh involving kl points to either side of each node.
type mesh
double precision :: M(2*i_KL+1, i_N)
end type mesh
For a matrix A, the element M(KL+1+n-j, j) = A(n,j). See the man page for the lapack
routine dgbtrf. The LU-factorisation of a banded matrix is also banded:
type lumesh
integer :: ipiv(i_N)
double precision :: M(3*i_KL+1, i_N)
end type lumesh
The element M(2*KL+1+n-j, j) = A(n,j). Pivoting information is stored along with the
matrix.
4.3.4 rdom the radial domain
The type (rdom) contains information the radial domain. The public variable mes_D
is declared in meshs.f90. Data includes the number of radial points N; powers of r,
r
p
n
=r(n,p); weights for integration rdr; matrices for taking the p
th
derivative dr(p);. . .
type rdom
integer :: N
double precision :: r(i_N,i_rpowmin:i_rpowmax)
double precision :: intrdr(i_N)
type (mesh) :: dr(i_KL)
end type rdom
4.4 Modications for the parallel implementation
The code has been parallelised using the Message Passing Interface (MPI). Sections of the
code with special MPI function calls are separated from the serial sections by preprocessor
directives. The number of processors is given by the parameter _Np set in parallel.h. The
rank of the local processor is given by mpi_rnk declared in mpi.f90, and mpi_sze=_Np.
Data type changes for the model equation (see equations 2.7, 3.18) are as follows:
Chapter 4. The code 18
q
*
s = p L
mesh
= g
lumesh
on radial subdomain
all real space values
for radial subdomain
all spectral coefficients
selection of modes
all radial points for coll
spec
transpose
transform
g
coll
spec
phys
multiply
invert
p
coll
coll
spec
phys
factorise
transform
for radial subdomain
all spectral coefficients
selection of modes
all radial points for
transpose
phys
Data is split over the spherical harmonics for the linear parts of the code evaluating
curls, gradients and matrix inversions for the timestepping. These linear operations do
not couple modes. Here all radial points for a particular mode are located on the same
processor, separate modes may be located on separate processors. In the denition of
type (coll) the parameter i_H1 is replaced by i_pH1 = (_Np+i_H1)/_Np-1. The variable
var_H has elements to determine which harmonics are on which processor. The elements
pH0 and pH1 are the zeroth index and one-minus the number of harmonics on the local
processor. The corresponding values for all processors are stored in the arrays pH0_(rank)
and pH1_(rank). For a variable of type (coll), valid values for the harmonic index in
Re(n,nh) are nh=0..pH1, which refer to the indices (4.2) nh=pH0..pH0+pH1.
Data is split radially when calculating Fourier transforms and when evaluating products
in real space. In the denitions of type (spec) and type (phys) the parameter i_N is
replaced by i_pN = (_Np+i_N-1)/_Np. The location of the radial subdomain on the local
processor is given by additional elements of the variable mes_D. The elements pNi and pN
are the location of the rst (inner) radial point and the number of points. The values
for all processors are stored in the arrays pNi_(rank) and pN_(rank). For a variable of
type (spec), valid values for the radial index in Re(nh,n) are n=1..pN which refer to the
indices of r
n
, where n=pNi..pNi+pN-1.
The bulk of communication between processors occurs during the data transposes in
the functions var_coll2spec() and var_spec2coll(). For _Np processors the data is
transfered in _Np1 steps. At each step each processor sends and receives a block of data,
the source and destination are selected systematically as set in the following loop:
do step = 1, mpi_sze-1
dest = modulo(mpi_rnk+step, mpi_sze)
src = modulo(mpi_rnk-step+mpi_sze, mpi_sze)
...
Chapter 4. The code 19
0 1 2 3
2
3
step=1
src / dest
data
sections
nr / nh
The sketch is for the case _Np= 4. The size of each block is approximately i_pN(i_pH1+1),
and the data within each block must also be transposed.
A
DIFFERENTIAL OPERATORS IN CYLINDRICALPOLAR
COORDINATES
Gradient
(f)
r
=
r
f, (f)
=
1
r
f, (f)
z
=
z
f. (A.1)
Laplacian
2
f = (
1
r
r
+
rr
)f +
1
r
2
f +
zz
f. (A.2)
Laplacian of vector
(
2
A)
r
=
2
A
r
2
r
2
A
r
r
2
,
(
2
A)
=
2
A
+
2
r
2
A
r
A
r
2
, (A.3)
(
2
A)
z
=
2
A
z
.
Divergence
A = (
1
r
+
r
)A
r
+
1
r
+
z
A
z
. (A.4)
Curl
( A)
r
=
1
r
A
z
z
A
,
( A)
=
z
A
r
r
A
z
, (A.5)
( A)
z
= (
1
r
+
r
)A
1
r
A
r
.
Directional derivative
(A B)
r
= A
r
r
B
r
+
A
B
r
+A
z
z
B
r
A
r
,
(A B)
= A
r
r
B
+
A
+A
z
z
B
+
A
B
r
r
, (A.6)
(A B)
z
= A
r
r
B
z
+
A
B
z
+A
z
z
B
z
.
Appendix A. Dierential operators in cylindricalpolar coordinates 21
Toroidalpoloidal decompositions
Axial
A = ( z) + ( z), (A.7)
= (r, , z), = (r, , z).
2
h
f = (
1
r
r
+
rr
)f +
1
r
2
f. (A.8)
A
r
=
1
r
+
rz
,
A
=
r
+
1
r
z
, (A.9)
A
z
=
2
h
.
Radial
A =
0
+
0
z + (r) + (r), (A.10)
0
=
0
(r),
0
=
0
(r), = (r, , z), = (r, , z).
2
c
f =
1
r
2
f +
zz
f. (A.11)
A
r
= r
2
c
,
A
=
0
+r
z
+
r
, (A.12)
A
z
=
0
+ (2 +r
r
)
z
.