# Error convergence on Stokes problem

**URL:** <https://community.freefem.org/t/error-convergence-on-stokes-problem/3550>\
**Category:** General Discussion\
**Created:** [October 17, 2024, 6:06am UTC](https://community.freefem.org/t/error-convergence-on-stokes-problem/3550 "2024-10-17T06:06:26Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![ggeedorah](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ggeedorah/32/2678_2.png) [@ggeedorah](https://community.freefem.org/u/ggeedorah)\
**Post date:** [October 17, 2024, 6:06am UTC](https://community.freefem.org/t/error-convergence-on-stokes-problem/3550/1 "2024-10-17T06:06:26Z")

</div>

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(pi_x)_sin(pi_y) - pi_sin(pi_x)_cos(pi_y);  
func f2 = 2\*(pi^2)_sin(pi_x)_sin(pi_y) - pi_cos(pi_x)_sin(pi_y);

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(pi_x)sin(piy);  
func ue2 = -sin(pi_x)_sin(pi_y);  
func pe = cos(pi\*x)_cos(pi_y);

func ue1x = pi_cos(pi_x)_sin(pi_y);  
func ue1y = pi_sin(pi_x)_cos(pi_y);  
func ue2x = -pi_cos(pi_x)_sin(pi_y);  
func ue2y = -pi_sin(pi_x)_cos(pi_y);

real nu = 1;

func f1 = -2\*(pi^2)_sin(pi_x)_sin(pi_y) - pi_sin(pi_x)_cos(pi_y);  
func f2 = 2\*(pi^2)_sin(pi_x)_sin(pi_y) - pi_cos(pi_x)_sin(pi_y);

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

---

<div class="post-metadata">

**Author:** ![fb77](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/fb77/32/3796_2.png) [@fb77](https://community.freefem.org/u/fb77)\
**Post date:** [October 17, 2024, 8:58am UTC](https://community.freefem.org/t/error-convergence-on-stokes-problem/3550/2 "2024-10-17T08:58:58Z")

</div>

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

> **[FreeFEM-documentation.pdf](https://doc.freefem.org/pdf/FreeFEM-documentation.pdf)**
>
> 33.81 MB

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)))`