load "PETSc-complex"
load "SLEPc"
load "Element_P3"
macro dimension()3L//
include "macro_ddm.idp"
load "lapack"
func Pb2 = [P2, P2];
func Pc = [P2, P2, P2,P1];
real ymin = 0.;
real ymax = 35.;
border yline(t=ymin,log(ymax+1)){x=0.;y=exp(t)-1;}
int np =2000;
if(!mpirank)
meshL thL=buildmeshL(yline(np));
broadcast(processor(0), thL);
meshL thLzero = thL;
fespace VhC(thL,Pc);
fespace Bh(thL,P2);
fespace VhBzero(thLzero,Pb2);
VhBzero<complex> [fG1,fG2];
/* ################################################################################## */
/* ######################## Solution of the Blasius equation ######################## */
/* ################################################################################## */
{
meshL thL2 = thL;
fespace VhB(thL2,Pb2);
VhB<complex> [f1, f2];
varf Resid([u1, u2], [v1, v2]) = int1d(thL2)(
-dy(f1)*v1 + f2*v1
+0.5*dy(f2)*f1*v2 - dy(f2)*dy(v2))
+ on(1, u1 = f1-0.) /* wall */
+ on(2, u2 = f2-1.) /* free stream */
+ on(1, u2 = f2-0.); /* wall */
varf Jacob([u1, u2], [v1, v2]) = int1d(thL2)(
-1.*dy(u1)*v1 + u2*v1
+0.5*dy(f2)*u1*v2 - dy(u2)*dy(v2) +0.5*dy(u2)*f1*v2)
+ on(1, u1 = f1-0.) /* wall */
+ on(2, u2 = f2-1.) /* free stream */
+ on(1, u2 = f2-0.); /* wall */
int[int] n2othL2;
macro thL2N2O()n2othL2//
Mat<complex> dJ;
{
macro def(u)[u, u#B]//EOM
macro init(i)[i, i]//EOM
MatCreate(thL2, dJ, Pb2);
}
set(dJ, sparams = "-ksp_type preonly -pc_type lu");
[f1, f2] = [y, 1];
func complex[int] funcRes(complex[int]& inPETSc) {
ChangeNumbering(dJ, f1[], inPETSc, inverse = true, exchange = true);
complex[int] out(VhB.ndof);
out = Resid(0, VhB, tgv = -1);
complex[int] outPETSc;
ChangeNumbering(dJ, out, outPETSc);
return outPETSc;
}
func int funcJ(complex[int]& inPETSc) {
ChangeNumbering(dJ, f1[], inPETSc, inverse = true, exchange = true);
dJ = Jacob(VhB, VhB, tgv = -1);
return 0;
}
complex[int] xPETSc;
ChangeNumbering(dJ, f1[], xPETSc);
SNESSolve(dJ, funcJ, funcRes, xPETSc, sparams = "-snes_monitor -snes_max_it 20 -snes_linesearch_monitor -snes_linesearch_order 2 -snes_atol 1e-10 -snes_rtol 1e-10 -snes_stol 1e-10 -snes_converged_reason");
ChangeNumbering(dJ, f1[], xPETSc, inverse = true, exchange = true);
int[int] rest = restrict(VhB, VhBzero, n2othL2);
f1[].re .*= dJ.D;
f1[].im .*= dJ.D;
VhBzero<complex> [fReduce1,fReduce2];
for[i, v : rest] fReduce1[][v] = f1[][i];
mpiAllReduce(fReduce1[], fG1[], mpiCommWorld, mpiSUM);
/* scaling: displacement thickness */
real Cm = 1.72089;
[fG1, fG2] = [fG1(x,Cm*y,z), fG2(x,Cm*y,z)];
}
/*#####################################################################################################################################################################################################*/
real Rey = getARGV("-Rey", 500.);
real kx = getARGV("-kx", 0.);
real kz = getARGV("-kz", 0.2);
real shiftr = getARGV("-shiftr", 0.4);
real shifti = getARGV("-shifti", 0.1);
complex shift = shiftr + 1i*shifti;
int nEV = getARGV("-nEV", 1);
int nKryl = getARGV("-nKryl", 300);
real nu = 1./Rey;
real k2= kx*kx+kz*kz;
varf A([u,v,w,p],[a1,a2,a3,a4])=int1d(thL)(a1*1i*kx*fG2*u+a1*dy(fG2)*v+ a1*1i*kx*p +nu*(a1*k2*u + dy(a1)*dy(u))
+a2*1i*kx*fG2*v +a2*dy(p) +nu*( a2*k2*v + dy(a2)*dy(v))
+a3*1i*kz*fG2*w + a3*1i*kz*p +nu*( a3*k2*w +dy(a3)*dy(w) )
+a4*(1i*kx*u + dy(v)+ 1i*kz*w))+on(1,2, u=0,v=0,w=0,p=0);
varf D([u,v,w,p],[a1,a2,a3,a4])=int1d(thL)(a1*u+a2*v+a3*w)+on(1,2, u=0,v=0,w=0,p=0);
macro def(u)[u, u#B, u#C,u#D]//EOM
macro init(i)[i, i, i,i]//EOM
matrix <complex> dA,dB;
dA = A(VhC, VhC, tgv = -2);
dB = D(VhC, VhC, tgv = -1);
matrix <complex> dAT=dA';
Mat <complex> J=dA;
Mat <complex> JT=dAT;
Mat<complex> Prod;
MatMatMult(J, JT, Prod);
real tol = 1e-16;
complex[int] nevEPS(nEV);
string ssparams = " -eps_type krylovschur" +
" -eps_nev" + nEV +
"-eps_tol " + tol +
" -eps_smallest_magnitude"+
" -ksp_type preonly -pc_type lu"
;
VhC<complex>[int] [ut,vt,wt,trt](nEV);
int epsK = EPSSolve(Prod, vectors = ut, values = nevEPS, sparams =ssparams);
Bh <complex> Ut,Vt,Wt,Tt;
Ut=ut[0];
Vt=vt[0];
Wt=wt[0];
I have tried to implement the method that you have showed in the screenshot; to find the gain using the (iw-L)(iw-L)*V=V\sigma^{-2} for Blasius Boundary Layer. But I am not getting the right values for gains and the forcings. Can you point out if it is due to wrong parameters in “ssparams” ? Also \omega for my case is 0 so I have not taken it.
