== Summary ==
First distrib version of Heated pipe code.

MediaWiki:Sidebar

2025-11-11T14:00:46Z

Apwillis:

MediaWiki:Sidebar

2025-08-28T13:31:38Z

Apwillis:

2025-02-07T10:35:58Z

Apwillis:

Download - Heated pipe

2025-02-07T10:32:04Z

Apwillis:

Related Codes

Getting started

2024-12-03T17:23:08Z

Apwillis: /* Typical setup commands */

{{latexPreamble}}

== Installation ==

If you just want a quick try, or if you are not running a linux environment, then you could try [http://www.virtualbox.org virtualbox] and the [[Xubuntu-Openpipeflow_virtual_disk_image]].

The following few commands are for installation in a linux environment. (On a Mac you could try 'brew install ...' [https://brew.sh/].)

If you've not already, [[download]] the code and unpack the tarball:
:~> tar -xvvzf Openpipeflow-x.xx.tgz

If you have administrator privilages, with linux try
:~> sudo apt-get install gfortran
:~> sudo apt-get install liblapack-dev
:~> sudo apt-get install libfftw3-dev
:~> sudo apt-get install libnetcdf-dev
:~> sudo apt-get install libnetcdff-dev
# optionally, for parallel MPI use:
:~> sudo apt-get install libopenmpi-dev

Then
:~> cd openpipeflow-x.xx/
:~> make
If this does not end with an error message, then you can skip the section on libraries! Next
:~> make clean

'''NOTE: If you wish to run in parallel, then you will need an mpi-fortran compiler. BUT, the code is sufficiently parallelized that it does NOT require parallel versions of LAPACK, netCDF or FFTW3. Each thread calls the serial implementation of each library.'''

'''Next''' try the [[Tutorial]], or take a look below if there were problems with libraries.

== Overview of files ==

[[File:Modules4.png]]

* '''<tt>Makefile</tt>''' is likely to require modification for your compiler and libraries (see [[#Libraries]]). It has been set up for <tt>gfortran</tt>, but compile flags for several compilers (<tt>g95, gfortran, ifort, pathf90, pgf90</tt>) can be found at the top of the file.
* '''<tt>parallel.h</tt>''' Ensure <tt>_Nr</tt> and <tt>_Ns</tt> are both set to <tt>1</tt> if you do not have MPI. This file contains macros for parallelisation that are only invoked if the number of processes <tt>_Np=_Nr*_Ns</tt> is greater than 1. Unless you need many cores, vary <tt>_Nr</tt> only, keep <tt>_Ns=1</tt> .
* '''<tt>program/parameters.f90</tt>'''. Reynolds number, resolution, timestep, etc. See [[#Parameters]]
* '''<tt>utils/</tt>''' contains a number of utilities. Almost anything can be done in a util, both post-processing and analysing data at runtime. There should be no need to alter the core code. See [[#Making_utils]]
* '''<tt>matlab/</tt>''', a few scripts. See <tt>matlab/Readme.txt</tt> .

== Makefile ==
A quick description of the Makefile below. There are suggested flags for several compilers at the top of this file, but the last group (compiler and paths) is the most important section if compiling for the first time. Note the include/link paths specified with <tt>-I<path></tt> and <tt>-L<path></tt>

INSTDIR = ./install/ [no need to change; where to put main.out and main.info files]
PROGDIR = ./program/ [no need to change; where code is located]
UTILDIR = ./utils/ [no need to change; where utility programs are located]
UTIL = prim2matlab [See [[#Making_utils]]; name of utility to build]

TRANSFORM = fftw3 [no need to change; name of FFT library]
MODSOBJ = io.o meshs.o mpi.o parameters.o \
timestep.o transform.o variables.o velocity.o

#COMPILER = g95 #-C [sample flags for several compilers]
#COMPFLAGS = -cpp -c -O3
#COMPILER = ifort -i-dynamic #-C #-static
#COMPFLAGS = -cpp -c -O3 -heap-arrays 1024 -mcmodel=medium
#COMPILER = pgf90 #-C
#COMPFLAGS = -Mpreprocess -c -fast #-mcmodel=medium
#COMPILER = pathf90 #-pg -C
#COMPFLAGS = -cpp -c -O3 -OPT:Ofast -march=opteron -fno-second-underscore

COMPILER = gfortran [***CHANGE TO mpif90 FOR PARALLEL USE***]
COMPFLAGS = -ffree-line-length-none -x f95-cpp-input -c -O3 \
-I/usr/include \ [***MIGHT NEED TO UPDATE INCLUDE PATH***]
#-C #-pg [ See info in following section ]
LIBS = \
-L/home/ash/lib \ [***MIGHT NEED TO UPDATE LIBRARY PATH***]
cheby.o -lfftw3 -llapack -lnetcdff -lnetcdf \
# -lblas -lcurl

== Libraries ==

'''The code uses freely availably libraries LAPACK, netCDF and FFTW3. If compiling for the first time, consider trying the <tt>gfortran</tt> + ''precompiled binaries'' combination.'''

'''NOTE: The code is sufficiently parallelized that it does not require parallel versions of LAPACK, netCDF or FFTW3. Each thread can calls the serial implementation.'''

In the openpipeflow directory, try <tt>make</tt>, looking for errors indicating missing libraries
> make
...
Fatal Error: Can't open module file 'netcdf.mod'
/usr/bin/ld: cannot find -llapack
/usr/bin/ld: cannot find -lfftw3

If you get something like
Error: Can't open included file 'mpif.h'
it is possible that <tt>_Np (=_Nr*_Ns)</tt> is not <tt>1</tt> in <tt>parallel.h</tt>, but the <tt>COMPILER</tt> has not been changed to <tt>mpif90</tt> in <tt>Makefile</tt>.

=== Precompiled binaries ===

If using a linux environment, install the LAPACK, FFTW3 and netCDF packages with your package manager or 'software center'. If you have sudo privileges, it may be sufficient to try
:~> sudo apt-get install ...
commands listed at the top of this page. Then
:~> cd openpipeflow-x.xx/
:~> make

If libraries are missing for netCDF, find the location of <tt>netcdf.mod</tt>
:~> cd /
:/> ls */netcdf.mod
:/> ls */*/netcdf.mod
/usr/include/netcdf.mod
Ensure that the file is in the include path in <tt>Makefile</tt>, here <tt>-I/usr/include</tt>

Often the package manager supplies several versions of a package, distinguished by extra extensions. If a library has been installed but is not found, then search for the library, e.g.
:/> ls */*/liblapack* */*/*/libfftw3* */*/*/libnetcdf*
etc/alternatives/liblapack.so.3
etc/alternatives/liblapack.so.3gf
usr/lib/i386-linux-gnu/libfftw3_omp.so.3
usr/lib/i386-linux-gnu/libfftw3_omp.so.3.3.2
usr/lib/i386-linux-gnu/libfftw3.so.3
usr/lib/i386-linux-gnu/libfftw3.so.3.3.2
usr/lib/i386-linux-gnu/libnetcdff.so.6
usr/lib/i386-linux-gnu/libnetcdf.so.11
Create a symbolic link to one of the versions for each package
:~> mkdir /home/ash/lib
:~> ln -s /etc/alternatives/liblapack.so.3 /home/ash/lib/liblapack.so
:~> ln -s /usr/lib/i386-linux-gnu/libfftw3.so.3 /home/ash/lib/libfftw3.so
:~> ln -s /usr/lib/i386-linux-gnu/libnetcdff.so.6 /home/ash/lib/libnetcdff.so
:~> ln -s /usr/lib/i386-linux-gnu/libnetcdf.so.11 /home/ash/lib/libnetcdf.so
Ensure that they are in the link path in <tt>Makefile</tt>, <tt>-L/home/ash/lib</tt>

'''Note:''' Some extra libraries may or may not be needed, depending on compiler/distribution. Try simply omitting the link flags
<tt>-lblas</tt> and/or <tt>-lcurl</tt> at first, for example.

=== Compiling libraries ===

If it is not possible to link with precompiled libraries, or if special flags are necessary, e.g. for a very large memory model,
then it may be necessary to build LAPACK and NetCDF with the compiler and compile flags that will be used for the main simulation code.

The default procedure for building a package (applicable to FFTW3 and netCDF) is

<pre>tar -xvvzf package.tar.gz
cd package/
[set environment variables if necessary]
./configure --prefix=<path to put libraries>
make
make install</pre>
– '''FFTW3'''. This usually requires no special treatment. Install with your package manager or build with the default settings.

– '''LAPACK'''. If using gfortran you might be able to use the binary version supplied for your linux distribution. Otherwise, edit the file make.inc that comes with LAPACK, setting the fortran compiler and flags to those you plan to use. Type ‘make’. Once finished, copy the following binaries into your library path (see Makefile LIBS -L$<$path$>$/lib/)<br />
<tt>cp lapack.a <path>/lib/liblapack.a</tt><br />
<tt>cp blas.a <path>/lib/libblas.a</tt>

– '''netCDF'''. If using gfortran you might be able to use the supplied binary version for your linux distribution. Several versions can be found at <tt>http://www.unidata.ucar.edu/downloads/netcdf/current/index.jsp</tt>.

Version '''4.1.3''' is relatively straight forward to install, the following typical environment variables required to build netCDF should be sufficient:
<br />
<tt>CXX=""</tt><br />
<tt>FC=/opt/intel/fc/10.1.018/bin/ifort</tt><br />
<tt>FFLAGS="-O3 -mcmodel=medium"</tt><br />
<tt>export CXX FC FFLAGS</tt><br />
After building, ensure that files <tt>netcdf.mod</tt> and <tt>typesizes.mod</tt> appear in your include path (see Makefile COMPFLAGS -I$<$path$>$/include/).

For more recent versions, netCDF installation is slightly trickier [currently netcdf-4.3.0.tar.gz 2015-07-20]. First<br />
<tt>./configure --disable-netcdf-4 --prefix=<path></tt><br />
which disables HDF5 support (not currently required, see comment on [[Parallel_i/o]]). Also, Fortran is no longer bundled, so get netcdf-fortran-4.2.2.tar.gz or a more recent version from here<br />
<tt>http://www.unidata.ucar.edu/downloads/netcdf/index.jsp</tt><br />
Build with the environment variables above, and in addition<br />
<tt>CPPFLAGS=-I<path>/include</tt><br />
<tt>export CPPFLAGS</tt><br />
Finally, add the link flag <tt>-lnetcdff</tt> immediately ''before'' <tt>-lnetcdf</tt> in your <tt>Makefile</tt>.

== Parameters ==

'''<tt>./parallel.h</tt>''':<br />
<tt> _Nr</tt> split in radius<br />
<tt> _Ns</tt> split axially<br />
Number of cores is <tt>_Np = _Nr * _Ns</tt>.

'''Set both <tt>_Nr</tt> and <tt>_Ns</tt> to <tt>1</tt> for serial use. MPI not required in that case.'''

Unless a large number of cores is required, it is recommended to vary only <tt>_Nr</tt> and to keep <tt>_Ns=1</tt>, i.e. split radially only. Prior to openpipeflow-1.10, only <tt>_Np</tt> was available, equivalent to the current <tt>_Nr</tt>.

<tt>_Nr</tt> is 'optimal' if it is a divisor of <tt>i_N</tt> or only slightly larger than a divisor.
To keep code simple, <tt>_Ns</tt> must divide both <tt>i_Z</tt> and <tt>i_M</tt> parameters below (easier to leave <tt>_Ns=1</tt> if not needed).

<br />
'''<tt>./program/parameters.f90</tt>''':<br />
<table>
<tr class="odd">
<td align="left"><tt>i_N</tt></td>
<td align="left">Number of radial points $n\in[1,N]$</td>
</tr>
<tr class="even">
<td align="left"><tt>i_K</tt></td>
<td align="left">Maximum k (axial), $k\in(-K,\,K)$</td>
</tr>
<tr class="odd">
<td align="left"><tt>i_M</tt></td>
<td align="left">Maximum m (azimuthal), $m\in[0,\,M)$</td>
</tr>
<tr class="even">
<td align="left"><tt>i_Mp</tt></td>
<td align="left">Azimuthal periodicity, i.e. $m=0,M_p,2M_p,\dots,(M-1)M_p$</td>
</tr>
<tr class="odd">
<td align="left"></td>
<td align="left">(set =1 for no symmetry assumption)</td>
</tr>
<tr class="even">
<td align="left"><tt>d_Re</tt></td>
<td align="left">Reynolds number $Re$ or $Rm_m$</td>
</tr>
<tr class="odd">
<td align="left"><tt>d_alpha</tt></td>
<td align="left">Axial wavenumber $\alpha=2\pi/L_z$</td>
</tr>
<tr class="even">
<td align="left"><tt>b_const_flux</tt></td>
<td align="left">Enforce constant flux $U_b=\frac1{2}$.</td>
</tr>
<tr class="odd">
<td align="left"><tt>i_save_rate1</tt></td>
<td align="left">Save frequency for snapshot data files (timesteps between saves)</td>
</tr>
<tr class="even">
<td align="left"><tt>i_save_rate2</tt></td>
<td align="left">Save frequency for time-series data</td>
</tr>
<tr class="odd">
<td align="left"><tt>i_maxtstep</tt></td>
<td align="left">Maximum number of timesteps (no limit if set <tt>=-1</tt>)</td>
</tr>
<tr class="even">
<td align="left"><tt>d_cpuhours</tt></td>
<td align="left">Maximum number of cpu hours</td>
</tr>
<tr class="odd">
<td align="left"><tt>d_time</tt></td>
<td align="left">Start time (taken from <tt>state.cdf.in</tt> if set <tt>=-1d0</tt>)</td>
</tr>
<tr class="even">
<td align="left"><tt>d_timestep</tt></td>
<td align="left">Fixed timestep (typically <tt>=0.01d0</tt> or dynamically controlled if set <tt>=-1d0</tt>)</td>
</tr>
<tr class="odd">
<td align="left"><tt>d_dterr</tt></td>
<td align="left">Maximum corrector norm, $\|f_{corr}\|$ (typically <tt>=1d-5</tt> or set <tt>=1d1</tt> to avoid extra corrector iterations) </td>
</tr>
<tr class="even">
<td align="left"><tt>d_courant</tt></td>
<td align="left">Courant number $\mathrm{C}$ (unlikely to need changing)</td>
</tr>
<tr class="odd">
<td align="left"><tt>d_implicit</tt></td>
<td align="left">Implicitness $c$ (unlikely to need changing)</td>
</tr>
</table>
'''Note the default cases''', usually if the parameter is set to -1.

== Input files ==

State files are stored in the NetCDF data format are binary yet can be transferred across different architectures safely. The program <tt>main.out</tt> runs with the compiled parameters (see <tt>main.info</tt>) but will load states of other truncations. For example, an output state file <tt>state0018.cdf.dat</tt> can be copied to an input <tt>state.cdf.in</tt>, and when loaded it will be interpolated if necessary.

Sample initial conditions are available in the [[Database]].

<tt>state.cdf.in</tt>:<br />
$t$ – '''Start time'''. Overridden by <tt>d_time</tt> if <tt>d_time</tt>$\ge$0<br />
$\Delta t$ – '''Timestep'''. Ignored, see parameter <tt>d_timestep</tt>.<br />
$N, M_p, r_n$ – '''Number of radial points, azimuthal periodicity, radial values''' of the input state.<br />
– If the input radial points differ from the runtime points, then the fields are interpolated onto the new points automatically.<br />
$u_r,\,u_\theta,\, u_z$ – '''Field vectors'''.<br />
– If $K\ne\,$<tt>i_K</tt> or $M\ne\,$<tt>i_M</tt>, then Fourier modes are truncated or zeros appended.

== Output ==

=== Snapshot data ===

Data saved every <tt>i_save_rate1</tt> timesteps:

<pre> state????.cdf.dat
vel_spec????.dat</pre>
All output is sent to the current directory, and <tt>????</tt> indicates numbers <tt>0000</tt>, <tt>0001</tt>, <tt>0002</tt>,…. Each state file can be copied to a <tt>state.cdf.in</tt> should a restart be necessary. To list times $t$ for each saved state file,<br />
<tt> > grep state OUT</tt><br />
The spectrum files are overwritten each save as they are retrievable from the state data. To verify sufficient truncation, a quick profile of the energy spectrum can be plotted with<br />
<tt> gnuplot> set log</tt><br />
<tt> gnuplot> plot 'vel_spec0002.dat' w lp</tt><br />

=== Time-series data ===

Data saved every <tt>i_save_rate2</tt> timesteps:

<table>
<tr class="odd">
<td align="left"><tt>tim_step.dat</tt></td>
<td align="left">$t$, $\Delta t$, $\Delta t_{\|f\|}$, $\Delta t_{CFL}$</td>
<td align="left">current and limiting step sizes</td>
</tr>
<tr class="even">
<td align="left"><tt>vel_energy.dat</tt></td>
<td align="left">$t$, $E$, $E_{k=0}$, $E_{m=0}$</td>
<td align="left">energies. Streamwise-dependent component of energy = $E-E_{k=0}$.</td>
</tr>
<tr class="odd">
<td align="left"><tt>vel_friction.dat</tt></td>
<td align="left">$t$, $U_b$ or $\beta$, $\langle u_z(r=0)\rangle_z$, $u_\tau$</td>
<td align="left">bulk speed or pressure measure $1+\beta=Re/Re_m$, mean centreline speed, friction vel.</td>
</tr>
<tr class="even">
<td align="left"><tt>vel_totEID.dat</tt></td>
<td align="left">$t$, $E/E_{lam}$, $I/I_{lam}$, $D/D_{lam}$</td>
<td align="left">total energy, input, dissipation; each normalised by laminar value.</td>
</tr>
</table>

== Typical usage ==

'''See the [[Tutorial]] for the fundamentals of setting up, launching, monitoring and ending a job.'''

==== Serial and parallel use ====

* To launch a serial job
> nohup ./main.out > OUT 2> OUT.err &
* The number of cores is set in <tt>parallel.h</tt> (all other parameters are in <tt>program/parameters.f90</tt>).
* For serial use, both <tt>_Nr</tt> and <tt>_Ns</tt> should be <tt>1</tt> in <tt>parallel.h</tt>.

* For parallel use, in <tt>Makefile</tt> the <tt>COMPILER</tt> needs to be changed to e.g. <tt>mpif90</tt>.
* MPI is required if, and only if, <tt>_Np=_Nr*_Ns</tt> is greater than <tt>1</tt> in <tt>parallel.h</tt>.
* For modest parallelisation, <tt>_Np</tt><=32, it is recommended to change <tt>_Nr</tt> only. Leave <tt>_Ns</tt> set to <tt>1</tt> (i.e. <tt>_Np==_Nr</tt>).
* For high parallelisation, <tt>_Nr</tt> and <tt>_Ns</tt> should be of the same order.
* To launch a parallel job on 8 cores
> nohup mpirun -np 8 ./main.out > OUT 2> OUT.err &
* Note that, rather than keep changing <tt>COMPILER</tt> when switching between parallel and serial use, it may be more convenient to simply launch with <tt>-np 1</tt> in the parallel launch command.

==== Initial conditions ====

The best initial condition is the state with most similar parameters. Some are supplied in the [[Database]].
Any output state, e.g. state0012.cdf.dat can be copied to <tt>state.cdf.in</tt> to be used as an initial condition. If resolutions do not match, they are '''automatically interpolated or truncated'''. If none is appropriate, then, it may be necessary to create a utility. See [[Utilities]].

==== Typical setup commands ====

> nano program/parameters.f90 [set up parameters and compile]
> make
> make install
> [prepare a directory for the job]
> mkdir ~/runs
> mv install ~/runs/job0009
> cd ~/runs/job0009
> [get input state.cdf.in from somewhere... e.g. another job]
> diff ../job0008/main.info main.info [compare parameters with a similar job - a good check]
> cp ../job0008/state0052.cdf.dat state.cdf.in
> [or get a state.cdf.in from
> http://www.openpipeflow.org/index.php?title=Database]
>
> [Run!...]
> nohup ./main.out > OUT 2> OUT.err & [run serial]
> nohup mpirun -np 8 ./main.out > OUT 2> OUT.err & [run parallel]
>
> head OUT [check for warnings/errors]
> tail OUT [see how far job has got]
>
> rm RUNNING [end the job cleanly]

If a job does not launch correctly, the problem can usually be found in the output
> less OUT ['q' to quit]
> less OUT.err

== Data processing ==

Almost anything can be done either at runtime or as post-processing by creating a utility.
It very rare that the core code in <tt>program/</tt> should need to be changed.

There are many examples in <tt>utils/</tt>. Further information can be found on the [[Utilities]] page.

Equations and parameters

2024-11-05T14:52:21Z

Apwillis: /* Non-dimensionalisation / scales */

{{latexPreamble}}
(As implemented in openpipeflow.org. For a reminder of the code parameter names see [[Getting_started#parameters]].)

== Governing equations ==

=== Non-dimensionalisation / scales ===

The scales used are
* $R$, the radius of the pipe. ($R=D/2$, where $D$ is the diameter.)
* $U_{cl}$, the centre-line velocity for laminar flow.
* $R/U_{cl}$ for time.

In computational units,
* the non-dimensional radius is 1,
* the non-dimensional laminar flow is $W(r)=1-r^2$,
* the non-dimensional bulk speed, when fixed, is $\frac{1}{2}$.

Note that when the mean axial flow speed $U_b$ is constant, we have $U_{cl}=2U_b$.

The streamwise wavenumber is $\alpha=2\pi/L$, where $L$ is the periodic length in units $R$. In papers where $D$ has been used as the scale, $\alpha$ may have been defined to give consistent values, i.e. $\alpha=\pi/L$.

For 'lab-units', based on $D$ and $U_b$, see [[Table of unit conversions]].
1 advection time unit $D/U_b$ is equivalent to 4 code time units $R/U_{cl}$.

=== Dimensionless parameters ===

Reynolds number, fixed flux, $Re_m = 2 U_b R / \nu = DU_b / \nu$, where the kinematic viscosity $\nu = \mu / \rho$.

Reynolds number, fixed pressure, $Re = U_{cl} R / \nu$.

For fixed flux (constant flow rate), $U_{cl}=2\,U_b$ at all times. The Reynolds number $Re_m$ is more commonly defined in terms of the constant mean speed $U_b$.

For fixed pressure gradient, $U_b$ is a time-dependent quantity that depends on the flow pattern. We define the Reynolds number $Re$ in terms of the unique $U_{cl}$ for the given pressure gradient.

$1+\beta = Re / Re_m$ is an observed quantity. For fixed flux, $1+\beta=\langle\partial p/\partial z\rangle \,/\, (dP/dz)$, where $(dP/dz)$ is the laminar pressure gradient and $\langle\partial p/\partial z\rangle$ is the average pressure gradient observed. For fixed pressure, $1+\beta=U_{cl}/(2U_b)$, where $U_b$ is the observed bulk speed.

The 'wall-Reynolds number' $Re_\tau=u_\tau R/\nu$, where $u_\tau$ is the 'wall-velocity' [https://en.wikipedia.org/wiki/Shear_velocity],
is given by $Re_\tau = (2\,Re_m\,(1+\beta))^\frac{1}{2}=(2\,Re)^\frac{1}{2}$.

=== Evolution equations ===

Fixed flux,

$ (\partial_{t} + \vec{u}\cdot\bnabla) \vec{u}
= -\bnabla \hat{p} + \frac{4}{\Rey_m}\,(1+\beta)\vechat{z}
+ \frac{1}{\Rey_m}\bnabla^2 \vec{u} $
and
$\bnabla\cdot\vec{u}=0$.

Fixed pressure

$ (\partial_{t} + \vec{u}\cdot\bnabla) \vec{u}
= -\bnabla \hat{p} + \frac{4}{\Rey}\vechat{z}
+ \frac{1}{\Rey}\bnabla^2 \vec{u} $
and
$\bnabla\cdot\vec{u}=0$.

Let $\vec{u}=W(r)\vechat{z}+\vec{u}'$. Using the scaling above, the laminar flow is $W(r) = 1-r^2$. The
equation, in rotational form, for the evolution of the perturbation $\vec{u}'$ is then

$ (\partial_{t} - \frac{1}{\Rey_m}\bnabla^2)\,\vec{u}'
= \vec{u}' \wedge (\bnabla \wedge\vec{u}')
- \frac{\mathrm{d}W}{\mathrm{d}r}\,u'_r \vechat{z}
- W\,\partial_{z}\vec{u}' + \frac{4\,\beta}{\Rey_m}\vechat{z} - \bnabla\hat{p}' \, . $

== Boundary conditions ==

The no-slip boundary conditions are $\vec{u}=\vec{0}$ at the wall, $r=1$.
There is no boundary condition explicitly on the pressure. Indirectly, the pressure must ensure that
$\bnabla\cdot\vec{u}=0$ is satisfied everywhere, i.e. also on the boundary.

At the axis $r=0$, symmetry implies that functions are odd or even across the axis. For a Fourier mode with azimuthal index $m$, each mode is odd/even if $m$ is odd/even for the variables $u_z$ and $p$ (and other scalars). For $u_r$ and $u_\theta$, each mode is even/odd if $m$ is odd/even.

== Decoupling the equations ==

The equations for $u_r$ and $u_\theta$ are coupled in the [[Differential_operators_in_cylindrical_coordinates#Laplacian|Laplacian]]. They can be separated in a Fourier decompositon by considering

$u_\pm = u_r \pm \mathrm{i} \, u_\theta,$

for which the $\pm$ are considered respectively. Original variables are easily recovered

$u_r = \frac{1}{2} ( u_+ + u_-),
\qquad
u_\theta = -\,\frac{\mathrm{i}}{2}(u_+ - u_- ) .$

Governing equations are then decoupled in the linear part and take the form

$\begin{eqnarray*}
(\partial_{t} - \nabla^2_\pm)\, u_\pm
& = & N_\pm - (\bnabla p)_\pm , \\
(\partial_{t} - \nabla^2 )\, u_z
& = & N_z - (\bnabla p)_z ,\end{eqnarray*}$

where

$\nabla^2_\pm = \nabla^2 - \frac{1}{r^2}
\pm \frac{2\,\mathrm{i}}{r^2}\partial_{\theta}$

Equations and parameters

2024-11-05T14:47:53Z

Apwillis: /* Non-dimensionalisation / scales */

{{latexPreamble}}
(As implemented in openpipeflow.org. For a reminder of the code parameter names see [[Getting_started#parameters]].)

== Governing equations ==

=== Non-dimensionalisation / scales ===

The scales used are
* $R$, the radius of the pipe.
* $U_{cl}$, the centre-line velocity for laminar flow.
* $R/U_{cl}$ for time.

In computational units,
* the non-dimensional radius is 1,
* the non-dimensional laminar flow is $W(r)=1-r^2$,
* the bulk velocity is $\frac{1}{2}$.

Note that $R=D/2$, where $D$ is the diameter. When the bulk flow rate is fixed, so that the mean axial flow speed $U_b$ is constant, we have $U_{cl}=2U_b$.

The streamwise wavenumber is $\alpha=2\pi/L$, where $L$ is the periodic length in units $R$. In papers where $D$ has been used as the scale, $\alpha$ may have been defined to give consistent values, i.e. $\alpha=\pi/L$.

For 'lab-units', based on $D$ and $U_b$, 1 advection time unit $D/U_b$ is equivalent to 4 code time units $R/U_{cl}$. See [[Table of unit conversions]].

=== Dimensionless parameters ===

Reynolds number, fixed flux, $Re_m = 2 U_b R / \nu = DU_b / \nu$, where the kinematic viscosity $\nu = \mu / \rho$.

Reynolds number, fixed pressure, $Re = U_{cl} R / \nu$.

For fixed flux (constant flow rate), $U_{cl}=2\,U_b$ at all times. The Reynolds number $Re_m$ is more commonly defined in terms of the constant mean speed $U_b$.

For fixed pressure gradient, $U_b$ is a time-dependent quantity that depends on the flow pattern. We define the Reynolds number $Re$ in terms of the unique $U_{cl}$ for the given pressure gradient.

$1+\beta = Re / Re_m$ is an observed quantity. For fixed flux, $1+\beta=\langle\partial p/\partial z\rangle \,/\, (dP/dz)$, where $(dP/dz)$ is the laminar pressure gradient and $\langle\partial p/\partial z\rangle$ is the average pressure gradient observed. For fixed pressure, $1+\beta=U_{cl}/(2U_b)$, where $U_b$ is the observed bulk speed.

The 'wall-Reynolds number' $Re_\tau=u_\tau R/\nu$, where $u_\tau$ is the 'wall-velocity' [https://en.wikipedia.org/wiki/Shear_velocity],
is given by $Re_\tau = (2\,Re_m\,(1+\beta))^\frac{1}{2}=(2\,Re)^\frac{1}{2}$.

=== Evolution equations ===

Fixed flux,

$ (\partial_{t} + \vec{u}\cdot\bnabla) \vec{u}
= -\bnabla \hat{p} + \frac{4}{\Rey_m}\,(1+\beta)\vechat{z}
+ \frac{1}{\Rey_m}\bnabla^2 \vec{u} $
and
$\bnabla\cdot\vec{u}=0$.

Fixed pressure

$ (\partial_{t} + \vec{u}\cdot\bnabla) \vec{u}
= -\bnabla \hat{p} + \frac{4}{\Rey}\vechat{z}
+ \frac{1}{\Rey}\bnabla^2 \vec{u} $
and
$\bnabla\cdot\vec{u}=0$.

Let $\vec{u}=W(r)\vechat{z}+\vec{u}'$. Using the scaling above, the laminar flow is $W(r) = 1-r^2$. The
equation, in rotational form, for the evolution of the perturbation $\vec{u}'$ is then

$ (\partial_{t} - \frac{1}{\Rey_m}\bnabla^2)\,\vec{u}'
= \vec{u}' \wedge (\bnabla \wedge\vec{u}')
- \frac{\mathrm{d}W}{\mathrm{d}r}\,u'_r \vechat{z}
- W\,\partial_{z}\vec{u}' + \frac{4\,\beta}{\Rey_m}\vechat{z} - \bnabla\hat{p}' \, . $

== Boundary conditions ==

The no-slip boundary conditions are $\vec{u}=\vec{0}$ at the wall, $r=1$.
There is no boundary condition explicitly on the pressure. Indirectly, the pressure must ensure that
$\bnabla\cdot\vec{u}=0$ is satisfied everywhere, i.e. also on the boundary.

At the axis $r=0$, symmetry implies that functions are odd or even across the axis. For a Fourier mode with azimuthal index $m$, each mode is odd/even if $m$ is odd/even for the variables $u_z$ and $p$ (and other scalars). For $u_r$ and $u_\theta$, each mode is even/odd if $m$ is odd/even.

== Decoupling the equations ==

The equations for $u_r$ and $u_\theta$ are coupled in the [[Differential_operators_in_cylindrical_coordinates#Laplacian|Laplacian]]. They can be separated in a Fourier decompositon by considering

$u_\pm = u_r \pm \mathrm{i} \, u_\theta,$

for which the $\pm$ are considered respectively. Original variables are easily recovered

$u_r = \frac{1}{2} ( u_+ + u_-),
\qquad
u_\theta = -\,\frac{\mathrm{i}}{2}(u_+ - u_- ) .$

Governing equations are then decoupled in the linear part and take the form

$\begin{eqnarray*}
(\partial_{t} - \nabla^2_\pm)\, u_\pm
& = & N_\pm - (\bnabla p)_\pm , \\
(\partial_{t} - \nabla^2 )\, u_z
& = & N_z - (\bnabla p)_z ,\end{eqnarray*}$

where

$\nabla^2_\pm = \nabla^2 - \frac{1}{r^2}
\pm \frac{2\,\mathrm{i}}{r^2}\partial_{\theta}$

Equations and parameters

2024-11-05T14:45:31Z

Apwillis: /* Non-dimensionalisation / scales */

{{latexPreamble}}
(As implemented in openpipeflow.org. For a reminder of the code parameter names see [[Getting_started#parameters]].)

== Governing equations ==

=== Non-dimensionalisation / scales ===

The scales used are
* $R$, the radius of the pipe.
* $U_{cl}$, the centre-line velocity for laminar flow.
* $R/U_{cl}$ for time.

In computational units,
* the non-dimensional radius is 1 and
* the non-dimensional laminar flow is $W(r)=1-r^2$.

Note that $R=D/2$, where $D$ is the diameter. When the bulk flow rate is fixed, so that the mean axial flow speed $U_b$ is constant, we have $U_{cl}=2U_b$.

The streamwise wavenumber is $\alpha=2\pi/L$, where $L$ is the periodic length in units $R$. In papers where $D$ has been used as the scale, $\alpha$ may have been defined to give consistent values, i.e. $\alpha=\pi/L$.

For 'lab-units', based on $D$ and $U_b$,
* 1 advection time unit $D/U_b$ is equivalent to 4 code time units $R/U_{cl}$,
* the bulk velocity corresponds to $\frac{1}{2}$ in code units.
[[Table of unit conversions]].

=== Dimensionless parameters ===

Reynolds number, fixed flux, $Re_m = 2 U_b R / \nu = DU_b / \nu$, where the kinematic viscosity $\nu = \mu / \rho$.

Reynolds number, fixed pressure, $Re = U_{cl} R / \nu$.

For fixed flux (constant flow rate), $U_{cl}=2\,U_b$ at all times. The Reynolds number $Re_m$ is more commonly defined in terms of the constant mean speed $U_b$.

For fixed pressure gradient, $U_b$ is a time-dependent quantity that depends on the flow pattern. We define the Reynolds number $Re$ in terms of the unique $U_{cl}$ for the given pressure gradient.

$1+\beta = Re / Re_m$ is an observed quantity. For fixed flux, $1+\beta=\langle\partial p/\partial z\rangle \,/\, (dP/dz)$, where $(dP/dz)$ is the laminar pressure gradient and $\langle\partial p/\partial z\rangle$ is the average pressure gradient observed. For fixed pressure, $1+\beta=U_{cl}/(2U_b)$, where $U_b$ is the observed bulk speed.

The 'wall-Reynolds number' $Re_\tau=u_\tau R/\nu$, where $u_\tau$ is the 'wall-velocity' [https://en.wikipedia.org/wiki/Shear_velocity],
is given by $Re_\tau = (2\,Re_m\,(1+\beta))^\frac{1}{2}=(2\,Re)^\frac{1}{2}$.

=== Evolution equations ===

Fixed flux,

$ (\partial_{t} + \vec{u}\cdot\bnabla) \vec{u}
= -\bnabla \hat{p} + \frac{4}{\Rey_m}\,(1+\beta)\vechat{z}
+ \frac{1}{\Rey_m}\bnabla^2 \vec{u} $
and
$\bnabla\cdot\vec{u}=0$.

Fixed pressure

$ (\partial_{t} + \vec{u}\cdot\bnabla) \vec{u}
= -\bnabla \hat{p} + \frac{4}{\Rey}\vechat{z}
+ \frac{1}{\Rey}\bnabla^2 \vec{u} $
and
$\bnabla\cdot\vec{u}=0$.

Let $\vec{u}=W(r)\vechat{z}+\vec{u}'$. Using the scaling above, the laminar flow is $W(r) = 1-r^2$. The
equation, in rotational form, for the evolution of the perturbation $\vec{u}'$ is then

$ (\partial_{t} - \frac{1}{\Rey_m}\bnabla^2)\,\vec{u}'
= \vec{u}' \wedge (\bnabla \wedge\vec{u}')
- \frac{\mathrm{d}W}{\mathrm{d}r}\,u'_r \vechat{z}
- W\,\partial_{z}\vec{u}' + \frac{4\,\beta}{\Rey_m}\vechat{z} - \bnabla\hat{p}' \, . $

== Boundary conditions ==

The no-slip boundary conditions are $\vec{u}=\vec{0}$ at the wall, $r=1$.
There is no boundary condition explicitly on the pressure. Indirectly, the pressure must ensure that
$\bnabla\cdot\vec{u}=0$ is satisfied everywhere, i.e. also on the boundary.

At the axis $r=0$, symmetry implies that functions are odd or even across the axis. For a Fourier mode with azimuthal index $m$, each mode is odd/even if $m$ is odd/even for the variables $u_z$ and $p$ (and other scalars). For $u_r$ and $u_\theta$, each mode is even/odd if $m$ is odd/even.

== Decoupling the equations ==

The equations for $u_r$ and $u_\theta$ are coupled in the [[Differential_operators_in_cylindrical_coordinates#Laplacian|Laplacian]]. They can be separated in a Fourier decompositon by considering

$u_\pm = u_r \pm \mathrm{i} \, u_\theta,$

for which the $\pm$ are considered respectively. Original variables are easily recovered

$u_r = \frac{1}{2} ( u_+ + u_-),
\qquad
u_\theta = -\,\frac{\mathrm{i}}{2}(u_+ - u_- ) .$

Governing equations are then decoupled in the linear part and take the form

$\begin{eqnarray*}
(\partial_{t} - \nabla^2_\pm)\, u_\pm
& = & N_\pm - (\bnabla p)_\pm , \\
(\partial_{t} - \nabla^2 )\, u_z
& = & N_z - (\bnabla p)_z ,\end{eqnarray*}$

where

$\nabla^2_\pm = \nabla^2 - \frac{1}{r^2}
\pm \frac{2\,\mathrm{i}}{r^2}\partial_{\theta}$

Table of unit conversions

2024-11-05T14:39:37Z

Apwillis:

{{latexPreamble}}

To convert the value of a dimensionless variable in 'code' units to 'lab' units, multiply by $C$ from the following table
('code' units are based on $R$ and $U_{cl}$, 'lab' units are based on $D$ and $U_b$).

E.g. for the variable $z$: $~~z_\mathrm{lab} \,(D) = z_\mathrm{code} \,(R) = z_\mathrm{code} \,(\frac{1}{2}D)
~~\Rightarrow~~z_\mathrm{lab} = z_\mathrm{code} \times \frac{1}{2},~~$
i.e. $~z_\mathrm{lab} = z_\mathrm{code} \times C$ with $C=\frac{1}{2}$.

Note: For consistent streamwise wavenumber $\alpha$, papers using the length scale $R$ write $L=2\pi/\alpha$, while papers using $D$ usually write $L=\pi/\alpha$.

$\begin{array}{ccccl}
variable & \mbox{'code' units} & \mbox{'lab' units} & \mbox{conversion factor}~C & \mbox{comment}\\
\hline
r,z & R & D & \frac{1}{2} & \mbox{length}\\
\vec{u} & U_{cl} & U_b & 2 & \mbox{speed} \\
t & R/U_{cl} & D/U_b & \frac{1}{4} & \mbox{time} \\
\sigma & U_{cl}/R & U_b/D & 4 & \mbox{growth rate} \\
E & \rho\,U_{cl}^2\, R^3 & \rho\,U_b^2\,D^3 & \frac{1}{2} & \mbox{kinetic energy} \\
D & \rho\,U_{cl}^3\, R^2 & \rho\,U_b^3\,D^2 & 2 & \mbox{dissipation rate} \\
E' & \rho\,U_{cl}^2\, R^2 & \rho\,U_b^2\,D^2 & 1 & \mbox{energy per unit length} \\
\tilde{E} & \rho\,U_{cl}^2 & \rho\,U_b^2 & 4 & \mbox{energy density}
\end{array}$

Core implementation

2024-03-27T22:59:45Z

Apwillis: /* How the parallelisation works */

{{latexPreamble}}

= Spatial representation =

== Finite difference stencil — $r$ ==

A function $p(r)$ is represented by the values $p_n=p(r_n)$, where $p(r)$ is evaluated on the $N$ radial points $r_n$; $n=1..N$.
There is no point on the axis $r=0$ to avoid singularities in $1/r$ terms.
The points are located at the roots of a Tchebyshev polynomial, such that they are clustered near the wall, and slightly so near the axis to compensate for loss of order in the calculation of derivatives. (They are calculated using banded matrices, implying use of fewer points at the axis and boundary).

=== Differentiation ===

Consider Taylor expansions about central point $x_0$, with $k$ neighbouring points each side

[[File:fd_deriv.png|500px]]

These expansions can be written

[[File:fd_deriv_vec.png||500px]]

Derivatives are calculated using weights from the appropriate row of $A^{-1}$.

=== Integration ===

Integrating the expansions above, from $x_0$ to $x$, the indefinite integral may be approximated about $x_0$:

[[File:fd_int.png|500px]]

The total integral is approximated by

$ \frac1{2} \sum_n \int_{x_{n-1}}^{x_{n+1}} f(x) \, {\mathrm d}x$

=== Matrix representation of linear operators ===

Linear equations can be written in the matrix form

$\mat{L} \,\vec{p} = \mat{M}\, \vec{q} + \vec{s}.$

$\mat{L}$ and $\mat{M}$ are $N\times N$ matrices representing linear operators. If $\vec{q}$ and $\vec{s}$ are known then the right-hand side is simple to evaluate. Consider the equation

$\label{eq:matmod}
\mat{L}\, \vec{p} = \vec{q},$

where $\vec{p}$ is unknown. If the matrix operators are formed by linear combinations of the weights in [[#Differentiation]], 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 $\mat{L}$.

figure=LU.eps, scale=0.65

== Fourier representaion — $\theta$, $z$ ==

=== Expansion of variables ===

$A(\theta,z) =
\sum_{|k|<K}\sum_{|m|<M}\, A_{km} {\mathrm e}^{{\mathrm i} (\alpha kz + m_0 m\theta)}$

$A$ real implies $A_{km} = A_{-k,-m}^*$ (coefficients are conjugate-symmetric). It is sufficient to keep only $m\ge 0$, and for $m=0$ keep $k\ge 0$.

=== Implied conditions on the axis ===

The geometry enforces conditions on each variable at the axis when expanded over Fourier modes in $\theta$. For scalars and the $z$-component of a vector $A_z$, each mode is even in $r$ if $m$ is even, and is odd if $m$ is odd. Each mode for $A_r$ and $A_\theta$ is even in $r$ if $m$ is odd, and is odd if $m$ is even. This implies that either the function or its derivative is zero on the axis.

=== Orthogonality ===

Exponentials

$\int_0^{2\pi} {\mathrm e}^{{\mathrm i}(m+n)} \, {\mathrm d}\theta
= 2\pi \, \delta_{m,-n}\, ,
\qquad
\int_0^{\frac{2\pi}{\alpha}} {\mathrm e}^{{\mathrm i}\alpha(k+j)} \, {\mathrm d}z
= \frac{2\pi}{\alpha} \, \delta_{k,-j}$

=== Volume integral ===

$E \, = \,
\int A \, {\mathrm d}V
\, = \,
\frac{4\pi^2}{\alpha} \, \int A_{00} \, r\,{\mathrm d}r$

=== Energy integral ===

$E \, = \,
{\textstyle \frac1{2}} \int \vec{A} \cdot \vec{A} \, {\mathrm d}V
\, = \,
\frac{2\pi^2}{\alpha} \, \sum_{km}
\int |\vec{A}_{km}|^2 \, r\,{\mathrm d}r$

$E = \sum_{m\ge 0} E_m,
\qquad
E_m =
\left\{
\begin{array}{ll}
\displaystyle
\frac{2\pi^2}{\alpha}\int\vec{A}_{00}^2 \, r \, {\mathrm d}r
+ \frac{4\pi^2}{\alpha} \sum_{k>0}
\int|\vec{A}_{k0}^2| \, r \, {\mathrm d}r,
& m=0, \\[15pt]
\displaystyle
\frac{4\pi^2}{\alpha} \sum_k
\int |\vec{A}_{km}|^2 \, r \, {\mathrm d}r,
& m>0 .
\end{array}
\right.$

= Temporal discretisation =

Temporal discritisation is via a second-order Predictor-Corrector scheme, with Euler predictor for the nonlinear terms and Crank-Nicolson corrector. The linear viscous term is treated implicitly with the Crank-Nicolson term.

== PPE-formulation with correct boundary conditions ==

An influence-matrix technique is used to avoid issues in satisfying the boundary conditions. See [[The_PPE_formulation]].

== Predictor-corrector ==

The model equation for each Fourier mode is

$\label{eq:harmmod}
(\partial_{t} - \nabla^2) f = N,$

where nonlinear terms have been evaluated on each radial point by the transform method, and is solved as in [[#Matrix_representation_of_linear_operators]]. The predictor at time $t_q$, with Euler nonlinear terms and implicitness $c$

$\frac{f_1^{q+1}-f^q}{\Delta t}
- \left( c\nabla^2 f_1^{q+1} + (1-c)\nabla^2 f^q \right)
= N^q ,$

$\left( \frac1{\Delta t} - c\nabla^2 \right) f_1^{q+1}
= \left( \frac1{\Delta t} + (1-c)\nabla^2 \right) f^q + N^q .$

Corrector iterations are

$\left( \frac1{\Delta t} - c\nabla^2 \right) f_{j+1}^{q+1}
= \left( \frac1{\Delta t} + (1-c)\nabla^2 \right) f^q
+ c\,N_j^{q+1} + (1-c)\, N^q ,$

or equivalently, for the correction $f_{corr} = f_{j+1}^{q+1} - f_j^{q+1}$

$\left( \frac1{\Delta t} - c\nabla^2 \right) f_{corr}
= c\,N_j^{q+1} - c\, N_{j-1}^{q+1} ,$

where $j=1,2,\dots$ and $N_0^{q+1}=N^q$. The size of the correction $\|f_{corr}\|$ must reduce at each iteration. For $c=\frac1{2}$ the scheme is second order such that $\|f_{corr}\| \sim \Delta t^2$. Stability is improved without observable degradation in accuracy using $c=0.501$ .

== Timestep control ==

$\Delta t = \mathrm{C} \, \min(\Delta \, / \, |\vec{v}| ) ,
\qquad
0<\mathrm{C}<1,$

where $\mathrm{C}$ is the Courant number.

The timestep should also be small enough such that the corrector norm $\|f_{corr}\|$ is satisfactorily small. This may be important for integrating initial transients, but ought not to be the limiting factor in general.

= The code structure =

[[File:modules4.png]]

== Overview of program files ==

* '''<tt>parameters.f90</tt>''' declaration and values of parameters.
* '''<tt>parallel.h</tt>''' macros for parallelisation; definition of loops for parallel+serial cases <tt>_loop_km_begin</tt> etc.
* '''<tt>mpi.f90</tt>''' mpi initialisation.
* '''<tt>main.f90</tt>''' core initialisation, main timestepping loop, clean up at end.
* '''<tt>meshs.f90</tt>''' definition of radial discretisation.
* '''<tt>variables.f90</tt>''' definition of derived types for spectral- and physical-space data.
* '''<tt>transform.fftw3.f90</tt>''' transform from collocated spectral space to physical space. Uses fftw3 library.
* '''<tt>timestep.f90</tt>''' subroutines to set up time stepping matrices.
* '''<tt>velocity.f90</tt>''' functions related to the velocity field; implementation of the timestepping and influence matrix method.
* '''<tt>io.f90</tt>''' Input-Output: loading and saving of data.

Each of the files contains a '<tt>module</tt>', declaring variables that may be used in other code via '<tt>use modulename</tt>'. With the exception of the fundamental parameters, the first three letters of the module name are used as a prefix in the name of public variables, to make the location of their declaration clear (see following section).

== Naming conventions ==

=== Parameters ===

Parameters defined within <tt>parameters.f90</tt> are given a prefix indicating the data type. For example:

<table>
<tr class="odd">
<td align="left">$N$</td>
<td align="left"><tt>i_N</tt></td>
<td align="left">integer</td>
</tr>
<tr class="even">
<td align="left">$\Delta t$</td>
<td align="left"><tt>d_timestep</tt></td>
<td align="left">double</td>
</tr>
<tr class="odd">
<td align="left">fixed flux?</td>
<td align="left"><tt>b_fixed_flux</tt></td>
<td align="left">boolean</td>
</tr>
</table>

Note that the naming convention frees the variable name for dummy variables, e.g.:

<pre> integer :: n
do n = 1, i_N
...</pre>

=== 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 prefix. Examples:

<table>
<tr class="odd">
<td align="left">$\Delta t$</td>
<td align="left"><tt>tim_dt</tt></td>
<td align="left"><tt>timestep.f90</tt></td>
</tr>
<tr class="even">
<td align="left">$u_r$</td>
<td align="left"><tt>vel_r</tt></td>
<td align="left"><tt>velocity.f90</tt></td>
</tr>
</table>

== Ordering the Fourier modes ==

$A = \sum_{|k|<K}\sum_{|m|<M}\, A_{km} {\mathrm e}^{{\mathrm i}(\alpha kz + m_0 m\theta)},$

Coefficients $A_{km}$ are stored in an order according to the following loop, where
rather than two indices <tt>k</tt> and <tt>m</tt>, the single index <tt>nh</tt> labels the modes:

[[File:km_loop_c.png|400px]]

The index nh is as in the following table for the case K=6, M=4:

[[File:Fourier_km.png|404px]]

NOTES:
* <tt>nh = (2*K-1)*m + k</tt>.
* Loops for <tt>m<0</tt> are skipped. Values of coefficients for <tt>m<0</tt> are inferred from the conjugate symmetric property, $A_{km}=A_{-k,-m}^*$. Similarly for <tt>k<0</tt> when <tt>m==0</tt>.
* This loop is a predefined macro in <tt>parallel.h</tt>. The macro <tt>_loop_km_begin</tt> replaces the double-loop above (or its equivalent for parallel computation).

== 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 defined as though we were working on a single processor. See [[#Modifications_for_the_parallel_implementation|Parallel]] regarding subtle differences in the definitions due to parallelisation.

The most important data types are <tt>type (coll) </tt> and <tt> type (phys)</tt>. The transform between the two types is achieved with calls to <tt>tra_coll2phys(...)</tt>, <tt>tra_phys2coll(...)</tt>.

=== coll, spec — storage of coefficients ===

The collocated data <tt>type (coll)</tt> is the principle type used for mode-independent operations.

<pre> type coll
double precision :: Re(i_N, 0:i_H1)
double precision :: Im(i_N, 0:i_H1)
end type coll</pre>
The element <tt>Re(n,nh)</tt> is the real part of the <tt>nh</tt>$^{\mathrm{th}}$ harmonic evaluated at $r_n$.

The spectral <tt>type (spec)</tt> is a data type for operations independent at each radial point. It is rarely used other than as a transitory data format between the <tt>coll</tt> and <tt>phys</tt> types.

<pre> type spec
double precision :: Re(0:i_H1, i_N)
double precision :: Im(0:i_H1, i_N)
end type spec</pre>

=== phys — real space data ===

The <tt>type (phys)</tt> is used for data in real space.

<pre> type phys
double precision :: Re(0:i_Z-1, 0:i_Th-1, i_N)
end type phys</pre>
The element <tt>Re(k,m,n)</tt> refers to the value at $z_k$, $\theta_m$ and $r_n$.

=== rdom — the radial domain ===

The <tt>type (rdom)</tt> contains information about the radial domain. The public variable <tt>mes_D</tt> is declared in <tt>meshs.f90</tt>. Data includes the number of radial points <tt>mes_D%N</tt>, powers of $r$, $r_n^p$=<tt>mes_D%r(n,p)</tt>, weights for integration
$\int_0^1 f(r)\,r\,{\mathrm d}r$=<tt>dot_product(f(1:i_N),mes_D%intrdr)</tt>, matrices for taking the $p^{\mathrm{th}}$ derivative and more:

<pre> 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</pre>

=== mesh, lumesh — storage of banded matrices ===

The finite-difference stencil has finite width ([[#Finite_differences_.E2.80.94_.24r.24|Finite_Differences]]) and so matrix operations involve matrices that are banded. The <tt>type (mesh)</tt>, defined in <tt>meshs.f90</tt>, is used for matrices on the radial mesh involving $kl$ points to either side of each node.

<pre> type mesh
double precision :: M(2*i_KL+1, i_N)
end type mesh</pre>
For a matrix <tt>A</tt>, the element <tt>M(KL+1+n-j, j) = A(n,j)</tt>. See the man page for the lapack routine <tt>dgbtrf</tt>. The LU-factorisation of a banded matrix is also banded:

<pre> type lumesh
integer :: ipiv(i_N)
double precision :: M(3*i_KL+1, i_N)
end type lumesh</pre>
The element <tt>M(2*KL+1+n-j, j) = A(n,j)</tt>. Pivoting information is stored along with the matrix.

== Modifications for the parallel implementation ==

=== Practical things to remember ===

The code has been parallelised using the Message Passing Interface (MPI) in a way designed to be as unintrusive as possible. 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 <tt>_Np=_Nr*_Ns</tt> set in <tt>parallel.h</tt>. The rank of a process is given by <tt>mpi_rnk</tt> declared in <tt>mpi.f90</tt>, and <tt>mpi_sze</tt>$=$<tt>_Np</tt>.

For modest parallelisations, it is recommended to vary <tt>_Nr</tt> and to keep <tt>_Ns=1</tt> (radial split only). This ensures a couple of conditions are met that could easily be forgotten when setting up a run, namely that <tt>_Ns</tt> must divide both <tt>i_Z</tt> and <tt>i_M</tt>.

=== How the data is distributed ===

The linear parts of the code (evaluating curls, gradients and matrix inversions for the timestepping) involve operations that do not couple Fourier modes. The <tt>type (coll)</tt> holds all radial points for a subset of Fourier modes - the Fourier modes are distributed over the cores.

Products are evaluated in physical space. Here the data is split over <tt>_Nr</tt> sections in $r$ and over <tt>_Ns</tt> sections in $z$. The <tt>type (phys)</tt> holds all azimuthal points for a
given subsection in $r$ and $z$. The following is an example of rank id numbers (<tt>mpi_rnk</tt>) for the case <tt>_Nr=3</tt>, <tt>_Ns=4</tt> ($r$ downwards, $z$ left to right):
0 3 6 9
1 4 7 10
2 5 8 11
In physical space,
* Two cores have data for the same $z$-section if <tt>rank1/_Nr==rank2/_Nr</tt> (integer arithmetic).
* Two cores have data for the same $r$-section if <tt>modulo(rank1,_Nr)==modulo(rank2,_Nr)</tt>.

=== How the code is adapted ===

In the definition of <tt>type (coll)</tt> the parameter <tt>i_H1</tt> is replaced by <tt>i_pH1</tt>.

In the definition of <tt>type (phys)</tt> the parameters <tt>i_N,i_Z</tt> are replaced by <tt>i_pN,i_pZ</tt>.

<tt>i_pH1</tt> is the maximum number of modes on an individual core. The variable <tt>var_H</tt> has elements to determine which harmonics are on which processor. The elements <tt>pH0</tt> and <tt>pH1</tt> are the zeroth index and the number of harmonics for the local process, minus one. The corresponding values for all processors are stored in the arrays <tt>pH0_(rank)</tt> and <tt>pH1_(rank)</tt>. For a variable of <tt>type (coll)</tt>, valid values for the harmonic index in <tt>Re(n,nh)</tt> are <tt>nh=0..pH1</tt> but refer to the indices <tt>nh=pH0..pH0+pH1</tt> of the serial case, see [[#Ordering_the_Fourier_modes]].

In physical space <tt>i_pN</tt> is the maximum number of radial points for each process. The location of the radial subdomain on the local processor is given by additional elements of the variable <tt>mes_D</tt>. The elements <tt>pNi</tt> and <tt>pN</tt> are the location of the first (inner) radial point and the number of points. The values for all processors are stored in the arrays <tt>pNi_(rank)</tt> and <tt>pN_(rank)</tt>. For a variable of <tt>type (spec)</tt>, valid values for the radial index in <tt>Re(nh,n)</tt> are <tt>n</tt>=<tt>1..pN</tt> which refer to the indices of $r_n$, where $n$=<tt>pNi..pNi+pN-1</tt>.

Also in physical space <tt>i_pZ</tt> is exactly the number axial points held by the core. Logically, the local index <tt>k</tt> in <tt>[0,i_pZ-1]</tt>corresponds to index <tt>(mpi_rnk/_Ns)*i_pZ + k</tt> (integer arithmetic) in <tt>[0,i_Z-1]</tt>.

The bulk of communication occurs during the data transposes in the functions <tt>var_coll2spec()</tt> and <tt>var_spec2coll()</tt>. For the case <tt>_Np=_Nr, _Ns=1</tt> the data is transferred in <tt>_Np-1</tt> 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:

<pre> do step = 1, mpi_sze-1
dest = modulo(mpi_rnk+step, mpi_sze)
src = modulo(mpi_rnk-step+mpi_sze, mpi_sze)
...</pre>

=== How the parallelisation works ===

In the Fourier space, all radial points for a particular Fourier mode are located on the same processor, separate modes may be located on separate processes. For the case <tt>K=6, M=4, _Nr=3, _Ns=2,</tt> the modes are split over <tt>_Np=_Nr*_Ns=6</tt> cores (<tt>mpi_rnk=0..5</tt>) indicated by the 6 coloured sections (each contiguous in <tt>nh</tt>) in the following table:

[[File:Fourier_km_par.png|404px]]

Here, the m=0,1,2,3 are split into for <tt>_Ns=2</tt> groups, m=0,1 and m=2,3, then each group is split into <tt>_Nr=3</tt> sections indicated by the colours.

For <tt>_Ns=2</tt>, split the 6 cores into two sets of three, denoted
$\{c_0,c_1,c_2\}$ and $\{c_3,c_4,c_5\}$.
Prior to an FFT, the core $c_0$ gathers all Fourier coefficients from within the first set for a given
subset of radial points -- doing this for each core constitutes a transpose with in the set.
Having gathered these coefficients, sufficient data is present to perform the axial FFTs
(summing over the k for a given m).
Now three sets of cores may be identified, arranged to have data for common
radial points in each set, $\{c_0,c_3\}$, $\{c_1,c_4\}$, $\{c_2,c_5\}$.
Within each set m is distributed, with
m=0,1 is on one core m=2,3 on the other. Gathering the m for a
given subsection of axial points, constitutes a second transpose within each of these
sets. Finally the azimuthal FFTs may be performed (sum over m).
The resulting data is split by both radial and axial position, whilst data for all azimuthal points are present on a given core.

In the double parallelisation two transposes are required, and the total data sent is
doubled. However, the number of messages required is
<tt>(_Ns-1)+(_Nr-1)</tt> where <tt>_Np=_Nr*_Ns</tt>.
Then, for a given $N_p=$<tt>_Np</tt>, there are only approximately $\sqrt{N_p}$ messages when <tt>_Nr</tt> and <tt>_Ns</tt> are approximately equal, versus $N_p$ messages when <tt>_Np=_Nr</tt> (and <tt>_Ns=1</tt>).
This can substantially reduce time lost in latency, the time setting up communications.
There are fewer large messages, rather than many small ones.