Description
Axisymmetric, hydromagnetic, modulated Couette flow between infinite or finite
cylinders in the small Prandtl number limit. Timestepping is via 2nd order
accurate implicit CrankNicolson for the linear terms and 2nd order accurate
explicit AdamsBashforth for the nonlinear terms. For the diffusive equations
the code uses operator factorisation to allow a tridiagonal system to be
solved; the linear algebra package LAPACK is used to solve the Poisson
equations associated with the stream function, current and magnetic field. The
spatial discretisation is via 2nd order accurate centred finite differences.
It includes a 'homotopy' parameter, tau
, to continuously deform the
boundaries from the infinite cylinder case (tau = 0
) to the finite cylinder
case (tau = 1
). A basic particle path subroutine to track the trajectory of
a particle is also included. Spatial modulation of the inner and/or outer
Reynolds numbers in the axial direction is possible to mimic a wavy cylinder
boundary.
Requirements
 LAPACK  the linear algebra package.
 BLAS  the basic linear algebra subprograms.
 A Fortran 90/95 compiler.
Files
 couette_mod.f90  Main program file.
 current.f90  Routines to solve the Poisson equation for the azimuthal current.
 derivs.f90  Routines to calculate derivatives.
 ic_bc.f90  Routines to set up initial and boundary conditions.
 io.f90  Routines to do with input/output.
 linear.f90  Routines to set up the linear parts of the RHS of the azimuthal velocity and vorticity fields.
 magnetic.f90  Routines to solve the Poisson equation for the azimuthal magnetic field.
 Makefile  Makefile to compile the code (see below).
 matrices.f90  Routines to set up matrices involved in timestepping and solving Poisson equations.
 nonlinear.f90  As linear.f90 but for the nonlinear terms.
 parameters.f90  Parameters to set (see below).
 README.md  This file.
 setup  Setup script to compile and set up ready to run in separate directory.
 solve.f90  Routines to solve for the azimuthal velocity and vorticity fields, as well as the Thomas algorithm.
 stream.f90  Routines to solve Poisson equation for stream function.
 variables.f90  Derived types and routines to do with variables.
Main parameters

alpha  Wavenumber (
2pi/wavelength
) in the infinite case. Set equal zero if using finite cylinders. 
gamma  Aspect ratio (
height/gap
) for the finite case. Set equal to2pi/alpha
if using infinite cylinders. 
eta  Radius ratio (
R1/R2
).  Q  Chandrasekhar number. Measure of the imposed magnetic field.

Re1,2  Reynolds number in
Re1,2(t)=Re1,2+Re1,2mod*cos(om1,2*t)
.  Re1,2_mod  Modulation amplitude as above.
 om1,2  Frequency of modulation as above.

Re_incr  How much
Re1
should be incremented or decremented when searching for critical values. 
growth_tol  How close two successive growth rates should be before
Re1
is altered in searching for critical values. 
dt  Time step. In general
10^4
is a good choice. Once the Reynolds numbers are large (>700, say) and/or the spatial resolution is increased significantly (>80 points in z) thendt = 10^5
is a better choice for stability. 
seed  Initial seed for initial conditions. Set small
O(10^10)
for calculating linear growth rates. SetO(10^3)
for nonlinear saturation. In practice this can be set to zero, since the boundary conditions of the azimuthal velocity at the cylinder wall(s) can start the flow.  end_time  Final viscous diffusion time.

tau_init  Initial value of homotopy parameter, tau.
0<=tau<=1
.tau = 0
=> infinite cylinders,tau = 1
=> finite cylinders.  tau_step  For steady case, how much tau should be increased after saturation at each tau.
 tau_end  For steady case, the final value of tau desired.
 nx_init  Number of radial grid points.
 nt_init  Number of azimuthal grid points for use when a 3D OpenDX isosurface is desired. This is not actually used in any computation in the code other than for the isosurface plots.

nz_init  Number of axial grid points. For infinite
2*nx
is sufficient. For finitegamma*nx
is usually required for full resolution.  save_rate  After how many time steps should velocities be saved?
 save_rate_2  After how many time steps should crosssections be saved?
 xsect_save  Should crosssections (surfaces) of fields be saved?
 save3d  Should a 3D isosurface be saved (OpenDX)?

iso_hel  If
save3d = .true.
should the isosurface be of the helicity? Ifiso_hel = .false.
then the stream function is saved. 
restart  Should we restart from a previous run? If
.true.
, fileend_state.dat
should be copied to restart directory. 
auto_tau  For steady case, should tau be automatically stepped after saturation at each tau. Set in conjunction with
tau_step
andtau_end
. 
auto_Re  For steady case, should
Re1
be stepped automatically? 
dec_Re  Can the critical value for
Re1
only be found by a quasistatic decrease ofRe1
? 
hyst_Re  Is the critical value of
Re1
in a hysteresis region? (Specifically for 1 and 2cell flows at very small aspect ratio).
Rarely used parameters

eps1,2  Amplitude of oscillation of
Re_1,2(t,z)
in axial direction. 
freq1,2  Frequency of oscillation of
Re_1,2(t,z)
in axial direction.  x_par_pos  Initial radial position of a particle in the fluid as a fraction of gap width.

z_par_pos  Initial axial position of a particle in the fluid as a fraction of
gamma
(finite) oralpha
(infinite).  save_part  Should a particle path be saved?
Makefile
Use the makefile to compile the code and handle module dependencies.
 OBJECT  The compiled executable name.
 OBJS  The object files that should be linked.
 FC  The Fortran compiler.
 FFLAGS  Compiler flags.
 LDFLAGS  Any extra flags required by the linker.
To run
Set parameters.f90
then run
source setup <directory>
which compiles code and copies parameters.f90
and moves OBJECT
to
<directory>
. ./couette_mod
runs code.
If restarting from a previous run then be sure to set restart = .true.
in
parameters.f90
, and have the file end_state.dat
from the preceding run in
the run/submit directory.
Errors are output if either:

end_state.dat
exists butrestart = .false.
or 
restart = .true.
butend_state.dat
does not exist.
If, at any time during a run, the file SAVE
is present in the run directory
(e.g. by giving the command touch SAVE
), the program will note this and
output data files for contours and surfaces at that time. Once the data files
are output SAVE
is removed. This can be done at any time and as many times
as desired during the run.
Any run should be cleanly stopped by removing the empty file RUNNING
which is
present in the run directory once the job has started. This avoids having to
use CtrlC
to interrupt the program. Using this procedure will ensure that
all runtime files are cleaned up and unneeded directories removed.
Output files
 end_state.dat  Saves time index, p, timestep, dt, and fields u, Z, psi, current and magnetic field for restart.
 energy.dat  Saves the kinetic energy due to CCF and the total kinetic energy including CCF at each time.
 max_psi.dat  Saves maximum value of psi (and vr, vz) over whole field. Mainly for use with IDL plots.
 particle.dat  Saves the radial and axial position of a particle at each time.

time_tau.dat  If
auto_tau = .true.
saves time at which each step in tau occured.  torque.dat  Saves the torques on the inner (G1) and outer cylinders (G2), the torque due to CCF (Gc), and the ratio G1/Gc. For steady states G1=G2.

u_growth.dat  Saves time, radial velocity (outflow, inflow), growth rate, axial velocity, stream function, azimuthal velocity, vorticity, azimuthal current and magnetic field, and Reynolds number. Velocities are saved at the points defined in
io.f90
 subroutinesave_growth
. 
p1234567.dat  Stream function field saved at the time defined by
(timestep)*(number following 'p')
in filename. (Only ifxsect_save == .true.
).  u1234567.dat  As above but for azimuthal velocity field.
 z1234567.dat  As above but for vorticity field.
 vr1234567.dat  As above but for radial velocity field.
 vz1234567.dat  As above but for axial velocity field.
 j1234567.dat  As above but for azimuthal current field.
 b1234567.dat  As above but for azimuthal magnetic field.

xsect1234567.dat  Crosssections of all fields except vorticity saved at the time defined as above. (Only if
xsect_save == .true.
).