Hello, I am facing an issue with my eigenmodes. I am attempting to solve a problem in a general case, which requires the use of complex numbers in my matrix. However, for the basic case, I only need real numbers. The code below highlights my problem.
In the first part, the code solves the problem and lists all the calculated eigenvalues for the non-complex case. It then plots the eigenmode associated with the eigenvalue having the largest real part.
In the second part, it repeats the same process but with a complex matrix (equations are the same). While it provides the same eigenvalues, the eigenmodes differ (refer to the pictures below). How can I resolve this issue to ensure that the complex case yields the same eigenmode as the non-complex one (which I know to be true)?
A sincere thanks to those willing to assist me.
load "SLEPc-complex"
load "PETSc-complex"
// Cylinder
real ri = 1;
real re = 50;
real h = 1;
real S=max(re,h)*2;
// Output results in a text file
bool output = 0;
// Number of eigenvalues
int nev=30;
real shift=0.0;
//////////////// Physical parameters /////////////////////////
real Rm=-156;
real eta10=10^5;
real etaGa=1;
func eta0 = 1.0;
func eta1 = eta10;
real aa, qr1, aai, aae, qri, qre, qr0, alt, rad;
real alpha = 0.16*pi;
qr0=cos(alpha);
qri=qr0;
aai=qri/sqrt(1-qri^2);
qre=qr0;
aae=qre/sqrt(1-qre^2);
func qr=qri*(x<1)+qre*(x>=1);
func qphi=sqrt(1-qr^2);
///////////////////////////////////////////////////////////////
// Meshing quality
int M=8;
// Refine function
func f = (S*atan((x-ri))+S*atan(2*y));
func eta0mesh=f;
func eta1mesh=f;
real etavac=10^5;
////////////////////////////////////////////////////////////////////////////////
verbosity=1;
bool debug=false;
////////////////////////////////////////////////////////////////////////////////////
border S1(t=-pi/2,pi/2){x=S*cos(t); y=S*sin(t); label=1;}
border A1(t=0, 1){x=0; y=S-t*(S-h/2); label=2;}
border A2(t=1,-1){x=0; y=h/2*t; label=2;} // Left cylinder
border A3(t=0, 1){x=0; y=-h/2-t*(S-h/2); label=2;}
border C1(t=0,1){x=re*t; y=h/2;label=3;} // Top cylinder
border C2(t=1,-1){x=re; y=h/2*t; label=3;} // Right cylinder
border C3(t=1, 0){x=re*t; y=-h/2; label=3;} // Bot cylinder
//Continue Horizontale
mesh Th=buildmesh(S1(2*M)+A1(2*M)+A2(M)+A3(2*M)+C1(10*M)+C2(M)+C3(10*M));
Th = adaptmesh(Th, [f, f], ratio=1.5, nbvx=100000);
plot(Th,wait=0, cmm="Mesh refined");
// Finiite Element Space
fespace Vh(Th,[P2,P2]);
fespace Vout(Th,P2);
Vh[p,t], [pv,tv];
Vout u,v,Pol,T,Pv,Tv;
// Speed across the cylinder
func ut = 2*Rm*atan(x)/pi;
macro lap(v,u) (dx(v)*dx(u)+dx(v)*u/x+dy(v)*dy(u)) //
varf a([Pol,T],[Pv,Tv]) = int2d(Th,4)(-lap(Pv,Pol)-eta0*lap(Tv,T))
+ int2d(Th,4)(Pv*dy(T)*eta1*qphi*qr/(eta0+eta1*qphi^2))
+ int2d(Th,0)(-lap(Pv,Pol)-etavac*lap(Tv,T))
+ int2d(Th,4)(-dy(Tv)*ut*(Pol/x+dx(Pol))+dx(Tv)*ut*dy(Pol))
+ int2d(Th,4)(-dy(Tv)*dy(T)*eta1*eta0*qr^2/(eta0+eta1*qphi^2))
+ int1d(Th,1)(-Pv*Pol*(2.0/S-1/sqrt(x^2+y^2)))
+ on(1,2,3,T=0)
+ on(2,Pol=0);
varf b([Pol,T],[Pv,Tv]) = int2d(Th,4)(Pv*Pol/(eta0+eta1*qphi^2)+Tv*T)
+ int2d(Th,4)(eta1*qr*qphi/(eta0+eta1*qphi^2)*dy(Tv)*Pol)
+ int2d(Th,0)(Tv*T+Pv*Pol/etavac);
real[int] ev(nev); // to store nev eigein value
real[int] iev(nev); // to store nev eigein value (im part)
Vh[int] [eV,eW](nev);
matrix AA = a(Vh,Vh,solver = sparsesolver);
matrix BB = b(Vh,Vh,solver = sparsesolver);
/////////////////////////////////////////////////////////////////////////
//////////////// EigenValue and Eigenmodes real ////////////////////////
/////////////////////////////////////////////////////////////////////////
real lambda;
int kk;
try{
kk=EigenValue(AA, BB, sigma=shift, which="LM", value=ev, vector=eV, tol=1e-9, maxit=0, ncv=0);
}
catch(...){
cout << "error" << endl ;
}
// Searching for the largest real eigenvalue lambda
kk=min(kk,nev);
lambda=real(ev[0]);
real idxLambda=0;
for (int i=1;i<kk;i++)
{
if (lambda < ev[i]) {
lambda = ev[i];
idxLambda=i;
}
}
// Plot the eignmodes associated with the largest real eigenvalue
plot(eV[idxLambda], cmm="Poloidal", value=true, wait=0);
plot(eW[idxLambda], fill=true, cmm="Toroidal", value=true, wait=0);
/////////////////////////////////////////////////////////////////////////
//////////////// EigenValue and Eigenmodes COMPLEX /////////////////////
/////////////////////////////////////////////////////////////////////////
complex[int] evc(nev); // to store nev eigein value
Vh<complex>[int] [eVc,eWc](nev);
matrix<complex> AAc = a(Vh,Vh,solver = sparsesolver);
matrix<complex> BBc = b(Vh,Vh,solver = sparsesolver);
real lambdac;
int kkc;
try{
kkc=EigenValue(AAc, BBc, sigma=shift, which="LM", value=evc, vector=eVc, tol=1e-9, maxit=0, ncv=0);
}
catch(...){
cout << "error" << endl ;
}
// Searching for the largest real eigenvalue lambda
kkc=min(kkc,nev);
lambda=real(ev[0]);
real idxLambdac=0;
for (int i=1;i<kk;i++)
{
if (lambdac < real(ev[i])) {
lambdac = real(ev[i]);
idxLambdac=i;
}
}
// Plot the eignmodes associated with the largest real eigenvalue
plot(eVc[idxLambdac], cmm="Poloidal", value=true, wait=0);
plot(eWc[idxLambdac], fill=true, cmm="Toroidal", value=true, wait=0);
As we can see the global shape is similar but the values are not the same, we are not symmetric in the non complex case while we are in the second one. the eigenvalues are respectively : +0.385965+0i (non-complex case) and +0.385965+2.8916e-16i in the complex case.
