View source for Newton.f90

{{latexPreamble}}

This page describes usage of the Newton solver for pipe flow.  See also the generic routine [[File:NewtonHook.f90]].

==== Inputs ====

This is the contents of a typical <tt>params.in</tt> file:
  0
  150 100 1 5
  3d-6
  -1d0 1d-7 1d+7
  1d-3 1d-6
You are unlikely to need to change these settings, except perhaps
the value 5, which is the number of leading eigenvalues that
are calculated if newton converges.
150 is the max number of GMRES iterations.  If it is taking this
many iterations (per Newton iteration), convergence is unlikely.
100 is the max number of Newton iterations.

A typical <tt>shifts.in</tt> file:
   10000
 0.  -1.  100.
Respectively: 
* ndts - Number of timesteps to take in period T,
* sth - shift in $\theta$,
* sz - shift in $z$,
* T - period.
A value <tt>-1.</tt> for <tt>sz</tt> informs the code to search for a 
suitable initial guess for the shift.

The timestep size follows from the inputs, T/ndts (=0.01 in the example) which is stretched slightly if 
T changes.  The shifts are shifts ''backwards'', i.e. if a travelling wave
has phase speed c (>0), then after time T we expect sz=c*T (>0).

The value T is arbitrary for a travelling wave, but T=10 usually works 
well.  If the number of GMRES iterations is observed to be large
(more than 50), try doubling T.  If it is small (less than 10), halve
T.

The initial guess should be in <tt>state.cdf.in</tt>.

==== Outputs ====

The output file of most interest is <tt>newton.dat</tt>:
        14679          150       279939
  0     0  0.21542E-02  0.21542E-03  0.25636E-01  0.84033E-01
  1    18  0.46324E-03  0.21542E-03  0.55202E-02  0.83918E-01
  2    19  0.12734E-03  0.83324E-04  0.15176E-02  0.83905E-01
  3    23  0.20499E-04  0.83324E-04  0.24427E-03  0.83920E-01
Along the top: ndts, max GMRES iterations, degrees of freedom.

Then columns:  Newton iteration, GMRES iterations taken,
$||g(x(T))-x(0)||$, $\delta$, 
$||g(x(T))-x(0)||\,/\,||x(0)||$ and $||x(0)||$, where $g(x)$ denotes the spatial shift of $x$ (backwards by <tt>sz</tt>).

The penultimate column, relative error, is the most important one.
If it drops below 3d-6, this is usually considered sufficient.
$\delta$ is the step size being taken (trust-region approach).

'''Other outputs''':

* <tt>shift.dat</tt> records the shift at each iteration.  Each line is like <tt>shifts.in</tt> preceded by the iteration number.

* <tt>state0???.cdf.dat</tt> is the state at each iteration.

* <tt>state1000.cdf.dat</tt> is the converged state.

* <tt>state1001.cdf.dat, state1002.cdf.dat...</tt> are eigenvectors.

Eigenvalues are found in <tt>arnoldi.dat</tt>.
If the descending eigenvalues are e.g.
   real, complex, real, real, complex, real,...
then the respective state numbers for the eigenvectors are
   1001, 1002+3, 1004, 1005, 1006+7, 1008,...
where 1002+3 are real and imaginary components respectively.

==== Continuation of solutions ====

If the file <tt>params_ct.in</tt> is present, the code will attempt
continuation from the input <tt>state.cdf.in</tt>.  A typical <tt>params_ct.in</tt> is:
  2 3 1
   -50.
On the first line, 2 and 3 are the minimum and maximum number of 
newton iterations per continuation step, and 1 indicates the 
parameter to be continued, 1==$Re$, 2==$\alpha$.
On the second line -50 is the size of the first step in the 
chosen parameter from the (assumed to be converged) input 
<tt>state.cdf.in</tt>.

The outputs from each continuation step are moved to directories
<tt>0001</tt>, <tt>0002</tt>, <tt>0003</tt>... and a brief summary
of the data in each is written to <tt>newton.dat_ct</tt>.  Its columns are
 number of continuation step (= directory number), 
 number of newton iterations required, 
 $Re$, $\alpha$, wavespeed $c$.

Continuation restarting from two previous steps is possible.
If <tt>state10.cdf.in</tt> is present, then for example, <tt>params_ct.in</tt>
  2 3 1
    1.5
asks the code to start from a linear extrapolation: if
<tt>state.cdf.in</tt>=$u_0$ and <tt>state10.cdf.in</tt>=$u_1$, then
it would restart with initial guess $u=u_1+ 1.5*(u_1-u_0)$.
The corresponding extrapolation for the initial shifts should be supplied
in <tt>shifts.in</tt>.

Note: actual shifts should be given for the starting point. Do not use the <tt>-1.</tt> flag in <tt>shifts.in</tt>.