USE cbmat
USE gnuplot
USE symbolic
kTsurhoni2d=1.508E-8
INTEGER ny=1200; ymax=30
ARRAY(-1..ny+1) OF REAL y
DO y(i)=ymax∗i/ny FOR ALL i
dy=ymax/ny
INTEGER CONSTANT vpos=1,upos=2,wpos=3,ppos=4,nvars=4,
contpos=1,umompos=2,wmompos=3,vmompos=4
ARRAY(-1..ny) OF REAL Psi=0,oldPsi
ARRAY(1+nvars..ny∗nvars) OF COMPLEX perturbation,oldpert,adjoint,adjdA
#define u(iy) perturbation[upos+(iy)∗nvars]
#define v(iy) perturbation[vpos+(iy)∗nvars]
#define w(iy) perturbation[wpos+(iy)∗nvars]
#define p(iy) perturbation[ppos+(iy)∗nvars]
#define dym(f,iy) [f(iy)-f(iy-1)]/dy
#define dyp(f,iy) [f(iy+1)-f(iy)]/dy
#define dyc(f,iy) [f(iy+1)-f(iy-1)]/2/dy
#define dy2(f,iy) [f(iy+1)-2∗f(iy)+f(iy-1)]/dy^2
#define laplacian(f) dy2(f,iy)-(alpha^2+beta^2)∗f(iy)
#define mum(f) [f(iy)+f(iy-1)]/2
#define mup(f) [f(iy+1)+f(iy)]/2
#define U(iy) dym(Psi,iy)
#define V(iy) [oldPsi(iy)-Psi(iy)]
#define d2Psi(iy) dy2(Psi,iy)
cont==dym(v,iy)-I∗alpha∗u(iy)-I∗beta∗w(iy)
umom==I∗omega∗u(iy)-I∗alpha∗U(iy)∗u(iy)+mum(d2Psi)∗mum(v)-I∗alpha∗p(iy)-ni∗laplacian(u)
vmom==I∗omega∗v(iy)-I∗alpha∗dyc(Psi,iy)∗v(iy)+dyp(p,iy)-ni∗laplacian(v)
wmom==I∗omega∗w(iy)-I∗alpha∗U(iy)∗w(iy)-I∗beta∗p(iy)-ni∗laplacian(w)
D_eps_umom==-dym(V,iy)∗u(iy)+mum(V)∗dyc(u,iy)
D_eps_vmom== V(iy)∗dyc(v,iy)+dyc(V,iy)∗v(iy)
D_eps_wmom== mum(V)∗dyc(w,iy)
REAL omega,beta,ni,x
COMPLEX alpha
ARRAY(1+nvars..ny∗nvars,-(2∗nvars-1)..2∗nvars-1) OF COMPLEX A,dA
INLINE FUNCTION AA(pos1,pos2)=A(pos1+nvars∗(∗),pos2-pos1+nvars∗(∗))
INLINE FUNCTION dAA(pos1,pos2)=dA(pos1+nvars∗(∗),pos2-pos1+nvars∗(∗))
SUBROUTINE BaseFlow()
REAL U0(0..ny)
Psi(0)=0; U0(0)=0; U0(1)=0.332/SQRT(x)∗y(1); U0(ny)=1
DO
U0(1)=U0(1)∗U0(ny)**(-1.5)
LOOP FOR i=1 TO ny-1
d1=2/[y(i+1)-y(i-1)]; d2p=d1/[y(i+1)-y(i)]; d2m=d1/[y(i)-y(i-1)]; d20=d2p+d2m
Psi(i)=Psi(i-1)+[U0(i-1)+U0(i)]∗[y(i)-y(i-1)]/2
U0(i+1)={d20∗U0(i)-[d2m-d1∗Psi(i)/4/x]∗U0(i-1)}/[d2p+d1∗Psi(i)/4/x]
REPEAT LOOP
FOR 2 TIMES
Psi(ny)=Psi(ny-1)+[U0(ny-1)+U0(ny)]∗[y(ny)-y(ny-1)]/2
END BaseFlow
SUBROUTINE BuildMats()
A=0; dA=0
LOOP FOR iy=1 TO ny-1 AND j=-1 TO 1 AND jv=1 TO nvars
eps=0
AA(contpos,jv,iy,j)=dy∗D(cont,perturbation[(iy+j)∗nvars+jv])
AA(umompos,jv,iy,j)=dy∗D(umom,perturbation[(iy+j)∗nvars+jv])
AA(vmompos,jv,iy,j)=dy∗D(vmom,perturbation[(iy+j)∗nvars+jv])
AA(wmompos,jv,iy,j)=dy∗D(wmom,perturbation[(iy+j)∗nvars+jv])
dAA(contpos,jv,iy,j)=dy∗D[D_alpha(cont),perturbation[(iy+j)∗nvars+jv]]
dAA(umompos,jv,iy,j)=dy∗D[D_alpha(umom),perturbation[(iy+j)∗nvars+jv]]
dAA(vmompos,jv,iy,j)=dy∗D[D_alpha(vmom),perturbation[(iy+j)∗nvars+jv]]
dAA(wmompos,jv,iy,j)=dy∗D[D_alpha(wmom),perturbation[(iy+j)∗nvars+jv]]
REPEAT
END BuildMats
COMPLEX FUNCTION CMAX(COMPLEX FUNCTION(INTEGER i) f; INTEGER imin..imax)
RESULT=0
LOOP FOR i=imin TO imax
IF ABS(RESULT)<ABS(f(i)) THEN RESULT=f(i)
REPEAT
END CMAX
SUBROUTINE InvIt()
DO
BuildMats
LUdecomp A
TYPEOF(perturbation) dAp=dA∗perturbation; perturbation=A\dAp
adjoint=adjdA/A; adjdA=adjoint∗dA
savepert=perturbation
dAp=dA∗perturbation; perturbation=A\dAp
adjoint=adjdA/A; adjdA=adjoint∗dA
dalpha=-(adjdA∗savepert)/(adjdA∗perturbation)
alpha=alpha+dalpha
perturbation=perturbation/CMAX[u(i),1..ny-1]
normfact=1/(adjdA∗perturbation)
adjoint=adjoint∗normfact
adjdA=adjdA∗normfact
WRITE ":" ./.
WHILE ABS(dalpha)>1.E-8
WRITE "." ./.
END InvIt
#define uadj(iy) adjoint[umompos+(iy)∗nvars]
#define vadj(iy) adjoint[vmompos+(iy)∗nvars]
#define wadj(iy) adjoint[wmompos+(iy)∗nvars]
COMPLEX FUNCTION hfa()
RESULT=adjdA∗(perturbation-oldpert)
LOOP FOR iy=2 TO ny-2
RESULT=~+dy∗[uadj(iy)∗D_eps_umom+vadj(iy)∗D_eps_vmom+wadj(iy)∗D_eps_wmom]
REPEAT
END hfa
REAL FUNCTION SpectRec()=dy∗SUM
4∗NORM[alpha∗uadj(iy)]+4∗NORM[dym(vadj,iy)]+4∗NORM[beta∗wadj(iy)]+
2∗NORM[dyc(uadj,iy)-I∗alpha∗mum(vadj)]+2∗NORM[dyc(wadj,iy)-I∗beta∗mum(vadj)]+
2∗NORM[I∗beta∗uadj(iy)+I∗alpha∗wadj(iy)]
FOR iy=2 TO ny-2
ni=1/2700
xmin=0.090909; xmax=1; nx=40
omin=0.02; omax=0.08; nomg=20
bmax=0.12; nbeta=10
dx=(xmax-xmin)/nx
domega=(omax-omin)/nomg
dbeta=bmax/(nbeta+0.5); bmin=0.5∗dbeta
REAL ums(0..nx,1..ny-1)=0
FILE spectr=CREATE("spectr.out")
COMPLEX Na(0..nomg,0..nx)
LOOP FOR omega=omin TO omax+domega/2 BY domega
LOOP FOR beta=bmin TO bmax+dbeta/2 BY dbeta
WRITE BY NAME omega,beta
perturbation=0; u(1)=1
adjdA=0; adjdA(umompos+nvars)=1
alpha=omega/0.3
x=xmin; BaseFlow; InvIt
COMPLEX h(0..nx); h(0)=1
Na(ROUND[(omega-omin)/domega],0)=1
REAL ampl=0
LOOP FOR ix=1 TO nx
x=x+dx
oldPsi=Psi
BaseFlow
oldalpha=alpha
oldpert=perturbation
oldadjoint=adjoint
oldadjdA=adjdA
InvIt
h(ix)=h(ix-1)∗EXP(-I∗(alpha+oldalpha)/2/ni∗dx-hfa())
IF beta=bmin THEN
Na(ROUND[(omega-omin)/domega],ix)=Na(ROUND[(omega-omin)/domega],ix-1)∗EXP(-I∗(alpha+oldalpha)/2/ni∗dx)
IF IMAG(oldalpha)<0 AND IMAG(alpha)>0 THEN
c=CMAX(1/Na(ROUND[(omega-omin)/domega],i),ix-1..ix)
Na(ROUND[(omega-omin)/domega],0..ix)=~∗c
END IF
END IF
ampl=~+dx/NORM(h(ix))∗SpectRec()
IF IMAG(oldalpha)<0 AND IMAG(alpha)>0 THEN
c=CMAX(1/h(i),ix-1..ix)
h(0..ix)=~∗c
ampl=ampl/NORM(c)
END IF
WRITE x,alpha,kTsurhoni2d∗ni∗ni∗ampl∗NORM(h(ix))
DO ums(ix,iy)=~+kTsurhoni2d∗ni∗ni∗ampl∗NORM(h(ix))∗NORM(u(iy))∗domega/PI∗dbeta/PI FOR ALL iy
REPEAT
WRITE TO spectr omega,beta,kTsurhoni2d∗ni∗ni∗ampl∗NORM(h(nx)),ampl,LOG(MOD(h(nx)))
REPEAT
WRITE TO spectr
REPEAT
FILE rmsfile=CREATE("rmsfile.out")
LOOP FOR ix=0 TO nx
DO WRITE TO rmsfile xmin+ix∗dx,y(iy),SQRT(ums(ix,iy)) FOR iy=1 TO ny-1
WRITE TO rmsfile
REPEAT
FILE maxfile=CREATE("maxfile.out")
LOOP FOR ix=0 TO nx
Reh=SQRT(xmin+ix∗dx)/ni
WRITE TO maxfile Reh,SQRT(MAX(ums(ix,∗))),4.6+7∗(Reh-1200)/1000,EXP[4.6+7∗(Reh-1200)/1000],MAXABS(Na(∗,ix))
REPEAT