Error convergence on Stokes problem

General Discussion
1

Hi, i’m trying to implement a code for a Stokes problem, but there is two things that i’m having trouble. First of all is the introduction of the quotient space L_0^2(\Omega) via Lagrange Multipliers, I tried the following

macro grad(f) [dx(f), dy(f)]//

mesh Th = square(20,20);

fespace Vh(Th,P2);
fespace Qh(Th,P1);
fespace Lh(Th,P1);

Vh u1, u2, v1, v2;
Qh p, q;
Lh mu, lambda;

real nu = 1;

func f1 = -2*(pi^2)sin(pix)sin(piy) - pisin(pix)cos(piy);
func f2 = 2*(pi^2)sin(pix)sin(piy) - picos(pix)sin(piy);

solve stokes([u1, u2, p, lambda],[v1, v2, q, mu]) =
int2d(Th)( nu * ( grad(u1)‘*grad(v1) + grad(u2)’*grad(v2) )
- p * (dx(v1) + dy(v2))
+ q * (dx(u1) + dy(u2))
)
- int2d(Th)( f1 * v1 + f2 * v2 )
+ int2d(Th)( lambda * q)
- int2d(Th)( mu * p)
+ on(1, 2, 3, 4, u1 = 0, u2 = 0)
;

plot([u1,u2], p, wait = true, value = true);

It seems okay but i wanted to double check. Also I’m trying to compute the seminorm H^1 error and L^2 for velocity and pressure, and the results are not the expected, here is the code

func ue1 = sin(pix)sin(piy);
func ue2 = -sin(pix)sin(piy);
func pe = cos(pi*x)cos(piy);

func ue1x = picos(pix)sin(piy);
func ue1y = pisin(pix)cos(piy);
func ue2x = -picos(pix)sin(piy);
func ue2y = -pisin(pix)cos(piy);

real nu = 1;

func f1 = -2*(pi^2)sin(pix)sin(piy) - pisin(pix)cos(piy);
func f2 = 2*(pi^2)sin(pix)sin(piy) - picos(pix)sin(piy);

ofstream datafile(“convergence_data_stokes.txt”);

for (int i = 0; i < 5; ++i){
int n = 10*(2^i);

  mesh Th = square(n,n);      

  fespace Vh(Th,P2);          
  fespace Qh(Th,P1);         
  fespace Lh(Th,P1);          

  Vh u1, u2, v1, v2;         
  Qh p, q;                  
  Lh mu, lambda;              

  solve stokes([u1, u2, p, lambda], [v1, v2, q, mu]) =
      int2d(Th)( nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2) ) 
               - p * (dx(v1) + dy(v2))                           
               - q * (dx(u1) + dy(u2))                            
               )
      - int2d(Th)( f1 * v1 + f2 * v2 )                           
      - int2d(Th)( lambda * q)                              
      + int2d(Th)( mu * p) 
      + on(1, 2, 3, 4, u1 = 0, u2 = 0)                       
      ;

  int Ndofs = Vh.ndof + Qh.ndof + Lh.ndof;           

  real errgraduL2 = sqrt(int2d(Th)((dx(u1)-ue1x)^2) + int2d(Th)((dy(u1)-ue1y)^2) + int2d(Th)((dx(u2)-ue2x)^2) + int2d(Th)((dy(u2)-ue2y)^2 ));  
  real errpL2 = sqrt(int2d(Th)((p-pe)^2));

  datafile << Ndofs << " " << errgraduL2 << " " << errpL2 << endl;

}

Any help would be appreciated

2

Hello,
Your Stokes formulation is not correct, because lambda and mu should be constant, instead of arbitrary P1 functions (space Lh).
One way of doing it is to build the matrix directly, adding just one degree of freedom to it. You can look at the Lagrange multiplier example, p.671 in the FreeFem documentation

In this example the condition \int u=0 is imposed.

From your code, another way could be to force lambda to be constant by adding a penalty term
+ int2d(Th)(tgv*(dx(lambda)*dx(mu)+dy(lambda)*dy(mu)))