Obtain Eigenvectors using SLEPc parallely

Dear FreeFEM Developer,
I try to solve the following eigenvalue problem using SLPEc parallelly.

The FE space for [u, v, p] is

fespace Xh(Th, [P2, P2, P1]); 

Note that the variables U and V are independent of spatial variables, so they can not be defined through fespace.

I solve the eigenvalue using

EPSSolve(A, B, vectors = vec, values = val, ....);

My question is that U and V are not defined in any finite element space, how can we compose vec in vectors = vec?
It seems impossible. I checked all the SLEPc examples provided by FreeFEM (link), but I didn’t find similar cases.

Can we store the eigenvector in an array, just like this

For example:
real[int, int] EigVec(Xh.ndof+2, num);
where num denotes the number of eigenvalues.

Then we get eigenvectors through

Xh [eu, ev, ep]; 
eu[] = EigVec[0:Xh.ndof-1][0];

where eu denotes the leading eigenvector.

Yes, you can use the array parameter instead of the vectors parameter.
Then, with one column of the array, you can use ChangeNumbering to go back to the fespace function you defined.

1 Like

Dear Dr. Pierre Jolivet,
Your answer is really helpful. Thanks again.

Best

Dear Dr. Pierre Jolivet,
According to your suggestion, I have replaced vectors with array:

real[int, int] EigVec(Xh.ndof+2, nev);
EPSSolve(A, B, array=EigVec, values=EigVal, sparams = params); 

The code runs well. Then I try to transfer the data using ChangeNumbering

real[int] EigSig(Xh.ndof+2);
ChangeNumbering(OpJ, EigSig, EigVec(:, 0), inverse=true, exchange = true);
fespace Vh(Th, P2);  fespace Ph(Th, P1);
Vh uh, vh; 
Ph ph;
if(mpirank == 0) EigSig.resize(EigSig.n - 2);
Xh [eu, ev, ep] = [uh, vh, ph];
eu[] = EigSig;
uh = eu; vh = ev; ph = ep; 

I find that the eigenvalues are correct but the eigenvectors are almost zero (1e-39).
Is there something wrong with the above codes?

Best,

I cannot tell for sure without running the actual code, but it looks wrong indeed:

  1. EigVec should not have the dimension of the local fespace, but the local dimension of the Mat (without ghost unknowns);
  2. ChangeNumbering should be called with a vector of dimension Xh.ndof first, and OpJ.n then. You are calling it with a vector of dimension Xh.ndof + 2 first, and Xh.ndof + 2 then;
  3. you should not resize EigSig after a call to ChangeNumbering.

See FreeFem-sources/laplace-lagrange-PETSc.edp at develop · FreeFem/FreeFem-sources · GitHub for an similar example.

Dear Dr. Pierre Jolivet,
Thanks for your kind reply, I will revise my code according to your suggestion and the attached example.