(* PPF-LinAn.m
   Pb. 2.1 Linear stability analysis of Plane Poiseuille Flow with a spectral method
   Partial solution by E Plaut
   Some ... or .... have to be replaced by instructions in Mathematica Language
 *)

<<commands.m

(************** Q1 **************)
(* Nb of modes *) Nz= 17; Print["Nz= ",Nz]

(* Base functions for Psi *)

F[n_,z_]= (z^2-1)^2 ChebyshevT[2 n-2,z]

Plot[{F[1,z],F[2,z],F[3,z],F[4,z]},{z,-1,1},
PlotStyle->{Black,Red,Green,Blue},
 Frame->True,FrameLabel->{"z","F"}]

(*Comments:
- seem to form a basis of the functions that vanish,
  with their first derivative, at z=+-1,
  i.e., the functions that fulfill the BC of the problem
- higher n functions oscillate more*)

(* Gauss-Lobatto Points... The . matters ! *)

z[m_]= Cos[m Pi/(2 Nz+1.)]

(* Operators *)

lz= Table[{0,z[m]},{m,Nz}]
gz= ListPlot[lz,PlotMarkers->{Graphics[Disk[{0,0},1]],.04}]

(*Comments:
- all points in the interval [0,1] ok since all terms in OS eq.
   are even, no need to test non positive values of z !
- correspond to a kind of mesh which is more dense
  near z = 1 which is were high gradients occur => good
*)

(************** Q3.a **************)
(* Base flow *)
U[z_]= 1-z^2

(* Operators BEWARE THE SAME form of the variable - a function -
   should appear on the left AND the right !
   Thus Dz2[f_]:= D[f[z],z,z] is a NONSENSE *)
   
Dz2[f_]:= D[f,z,z];
Del[f_]:= -k^2 f + Dz2[f]
Dop[f_]:= -Del[f]
Lop[f_]:= -1/R Del[Del[f]] + I k (U[z] Del[f] - U''[z] f)

MatD = MatL = IdentityMatrix[Nz];

Do[
   Do[
      MatD[[m,n]] = ReplaceAll[Dop[F[n,z]], z->z[m]];
      MatL[[m,n]] = ReplaceAll[Lop[F[n,z]], z->z[m]] ;
     ,{m,1,Nz}]
  ,{n,1,Nz}]

(************** Q3.b **************)

MD[k_] = MatD;  ML[k_,R_] = MatL;

(************** Q4 **************)

spectrum[k_,R_]:= Eigenvalues[{ML[k,R], MD[k]}]

sigma1[k_?NumericQ,R_?NumericQ]:= Last[Sort[spectrum[k,R]]]

(* Evidence the instability at high R *)
print["sigma1[1,1000]"]
print["sigma1[1,2000]"]
print["sigma1[1,....]"]
print["sigma1[1,....]"]
print["sigma1[1,....]"]
print["sigma1[1,....]"]

(************** Q5 **************)

fr= FindRoot[Re[sigma1[1, R]]==0, {R, ...., ....}]

R0[k_?NumericQ]:= R/.FindRoot[Re[sigma1[k, R]]==0, {R, ...., ....}]

lkR0= Table[{k, R0[k]},{k,.95,1.09,.01}]

fs= Directive[Thickness[.004],Black];

(*taille des ticks*) tt= .019;
tk= Table[{k,Tost[k],{tt,0}},{k,.96,1.1,.02}];
tr= Table[{r,Tost[r],{tt,0}},{r,5000,7000,500}];

tksans= Map[nolabel,tk];
trsans= Map[nolabel,tr];
ft= {tk,tr,tksans,trsans};

pr= {{.95,1.1},{5000,7000}};
gkR0= ListLinePlot[lkR0,
 PlotStyle->{Thick,Black},PlotRange->pr];

fl= {Style["k",Italic,Black,16],Style["R",Italic,Black,16]};
GkR0= Show[gkR0,
 Graphics[{Text[Style["stable",Bold,15],{1.07,5500}],
  Text[Style["unstable",Bold,Red,15],{1,6500}]}],
 Frame->True,FrameStyle->fs,FrameTicks->ft,
 FrameLabel->fl,RotateLabel->False]

tk= {{1,Tost[1],{tt,0}}};
tksans= Map[nolabel,tk];
tr= {{6000,Tost[6000],{tt,0}}};
trsans= Map[nolabel,tr];
ft= {tk,tr,tksans,trsans};
gkR00= ListLinePlot[lkR0,PlotStyle->White,PlotRange->pr];

GkR00= Show[gkR00,
 Frame->True,FrameStyle->fs,FrameTicks->ft,
 FrameLabel->fl,RotateLabel->False]

Export["images/PPF-LA-gkR0-nosol.pdf",GkR00,ImageSize->300]

(************** Q6.a **************)

fm= FindMinimum[R0[k],{k,....,....}]

Rc= fm[[1]]
kc= k/.fm[[2]]

Print["kc= ",kc," , Rc= ",Rc]

...
