SLEPc different values with -np 1 vs -np 2+

General Discussion
1

Dear forum,

I am using FreeFEM++ v4.8 for Windows to perform a stability analysis with the SLEPc library. If I set -np 1 and I repeat the simulations without changing anything else I get everytime the same results. When I change -np 2 and I repeat the simulations without changing anything else I get different results and I don’t understand why. I attach the code without the equations:

mesh Th=…
fespace XXXXXFES(Th,[P1,P1,P1,P1]);
XXXXXFES[u,v,w,p];
XXXXXFES[ue,ve,we,pe];
XXXXXFES[hy,hx,hz,q];
func E = [P1, P1, P1, P1];
varf LNSE ([ue,ve,we,pe],[hy,hx,hz,q])
=int2d(Th)…
varf b([ue,ve,we,pe],[hy,hx,hz,q])
=int2d(Th)…
Mat M;
macro def(j) [j, j#A, j#B, j#C]
macro init(j) [j, j, j, j]
createMat(Th, M, E);
M = LNSE(XXXXXFES, XXXXXFES);
matrix B = b(XXXXXFES, XXXXXFES);
Mat J(M, B);
int nev = 100;
int ncv = 2*nev;
real tol = 1e-12;
complex[int] EigenVal(nev);
XXXXXFES[int] [EigenVecu,EigenVecv,EigenVecw,EigenVecp] (nev);
complex Shift = getARGV(“-shift_real”, 0) + getARGV(“-shift_imag”, 10) * 1i; /
string EVSParam =
“-eps_tol " + tol +” "+
“-eps_nev " + nev +” "+
"-eps_type krylovschur " +
"-st_type sinvert " +
"-eps_monitor_conv " +
"-eps_target " + real(Shift) + “+” + imag(Shift) + “i”;
int i, k;
k = EPSSolve(M, J, values=EigenVal, vectors=EigenVecu, sparams = EVSParam);

Did someone encounter this problem before?
Regards,
Robert

2

Near impossible to help you out without a runnable piece of code. What do you mean by “different results”? Consider upgrading your FreeFEM installation as well.

3

Also note that your finite element space consisting of P1-P1 elements is not inf-sup stable for Stokes problems. Try Taylor-Hood (P2-P1) or MINI (P1b-P1) elements.

4

Hello,
Thank you both for suggestions. I did some more tests and the problem might come from the boundary conditions as I have an axi-symmetrical problem and I use the weak form of the N-S equations written in the cylindrical components. By different results I mean that I just run the simulation to compute 100 eigenvalues. Instead of getting the exactly 100 eigenvalue I get different values each time, not even close to the previous simulation. In another case where I use the weak form of the N-S equation written in the cartesian components I always get the same eigenvalues.

I have to put as a boundary condition on the symmetry-axis the partial derivative of a flow component in respect to a direction equal to 0, du/dx=0. Is it correct to write it as it follows?

+int1d(Mesh, SymmetryAxis)((dx(u)*TestFunction)*1e30)

Regards,
Robert

5

No, it does not look like a proper way to impose such a boundary condition.

6

I recommend to check out the StabFem library. The StabFem gitlab repository contains various examples which can be very helpful to guide the correct implementation of the boundary conditions.

7

Dear @aszaboa ,
Thank you for the link, indeed is a very nice repository. Actually in this example, on line 148 I saw that kind of boundary condition.
Regards,
Robert