# Uzawa scheme for stokes bingham flow

**URL:** https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233
**Category:** General Discussion
**Created:** [May 9, 2024, 11:10am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233 "2024-05-09T11:10:13Z")
**Posts on this page:** 20
**Page:** 1

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 9, 2024, 11:10am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/1 "2024-05-09T11:10:13Z")

</div>

Hi all,  
I have written uzawa scheme for Bingham stokes flow on freefem. It is supposed to change when the bn number changes. However when I compile with different bn numbers it does not change. Could anyone help me?  
Here is the code  
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;  
real L = 1; //Pipe length  
real D = 1.; //Pipe height  
real L1=.1;  
real L2=.9;  
real rho=1.; real alpha=1.;  
macro Gradient(u) [dx(u), dy(u)] //  
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //  
macro UgradV(ux,uy,vx,vy) [[ux,uy]'_[dx(vx),dy(vx)] , [ux,uy]'_[dx(vy),dy(vy)] ]//  
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;  
func f2=1.;  
border b1(t=0., 1.){x=t; y=0.;};  
border b2(t=0., 1.){x=L; y=D_t; };  
border b3(t=L, L2){x=t; y=D;};  
border b4(t=L2, L1){x=t; y=D;};  
border b5(t=L1, 0.){x=t; y=D;};  
border b6(t=0., 1.){x=0.; y=D-D_t;};  
int nnL = 50;  
int nnD = 50;  
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6(nnD));  
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);  
fespace Wh(Th,P2);  
fespace Qh(Th,P1);  
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;  
Wh ux,uy;  
Qh p, q, dp,p1;

fespace Lh(Th,P1);  
Lh TmGxx, TmGyy, TmGxy;  
Lh Txx, Tyy, Txy;  
Lh Gdotxx, Gdotyy, Gdotxy;  
Lh dTxx=0., dTyy=0., dTxy=0.;  
Lh TmGxx1, TmGyy1, TmGxy1;  
Lh TmGxx0, TmGyy0, TmGxy0;  
func u5= 1.;  
func u4=1.;  
func u3=1.;

TmGxx=0.;  
TmGyy=0.;

```
TmGxy=0.;

for(int j=0; j<100; ++j)
{  

```

solve bing([ux,uy],[vx,vy])  
=int2d(Th)((Mu)\*(Gradient(ux)‘\*Gradient(vx)  
+ Gradient(uy)’\*Gradient(vy)  
)  
)

```
+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
)
 -int2d(Th)([f1 , f2]' *[vx ,vy]
	)
	
	+ on(b1, ux=0., uy=0.)
			+ on(b2, ux=0., uy=0.)
            + on(b3, ux=u3, uy=0.)
			+ on(b4, ux=u4, uy=0.)
			+ on(b5, ux=u5, uy=0.);

Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
	real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
		 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
		 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
	real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
	real res = TaGnorm;
	if( 1.<res ){
		
		real q = 1./res;
		
		TmGxx[][k] = TaGxx*q;
		TmGyy[][k] = TaGyy*q;
		TmGxy[][k] = TaGxy*q;
	}
	else {
		
		TmGxx[][k] = TaGxx;
		TmGyy[][k] = TaGyy;
		TmGxy[][k] = TaGxy;
	}
}ux0[]=ux[];
uy0[]=uy[];

```

}  
fespace Rh(Th, P2);  
Rh dyg,dyg2;  
solve streamlines (dyg, dyg2,solver=UMFPACK)  
= int2d(Th)(  
dx(dyg)\*dx(dyg2)  
+ dy(dyg)_dy(dyg2)  
)  
+ int2d(Th)(  
- dyg2_(dy(ux) - dx(uy))  
)  
+ on(1,2,3,4,5,6, dyg=0);  
plot(dyg,ps=“ldsre1500.eps”,nbiso=10);  
Here is the uzawa bingham scheme  
laplace u^n+1+divergence lamda^n=f  
lamdaa^n+1=P(lamda^n+D(n^n+1))  
P(q)=q if |q|\<1 p(q)=q/|q| if |q|\>1

---

<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: [May 9, 2024, 11:45am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/2 "2024-05-09T11:45:16Z")

</div>

Can you please put all your line codes between:  
three backquotes at the begining,  
three backquotes at the end  
so that it appears unaltered to us?

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 9, 2024, 12:17pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/3 "2024-05-09T12:17:30Z")

</div>

Thank you has it changed?  
Here is the code

‘’’  
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;  
real L = 1; //Pipe length  
real D = 1.; //Pipe height  
real L1=.1;  
real L2=.9;  
real rho=1.; real alpha=1.;  
macro Gradient(u) [dx(u), dy(u)] //  
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //  
macro UgradV(ux,uy,vx,vy) [[ux,uy]'_[dx(vx),dy(vx)] , [ux,uy]'_[dx(vy),dy(vy)] ]//  
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;  
func f2=1.;  
border b1(t=0., 1.){x=t; y=0.;};  
border b2(t=0., 1.){x=L; y=D_t; };  
border b3(t=L, L2){x=t; y=D;};  
border b4(t=L2, L1){x=t; y=D;};  
border b5(t=L1, 0.){x=t; y=D;};  
border b6(t=0., 1.){x=0.; y=D-D_t;};  
int nnL = 50;  
int nnD = 50;  
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6(nnD));  
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);  
fespace Wh(Th,P2);  
fespace Qh(Th,P1);  
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;  
Wh ux,uy;  
Qh p, q, dp,p1;

fespace Lh(Th,P1);  
Lh TmGxx, TmGyy, TmGxy;  
Lh Txx, Tyy, Txy;  
Lh Gdotxx, Gdotyy, Gdotxy;  
Lh dTxx=0., dTyy=0., dTxy=0.;  
Lh TmGxx1, TmGyy1, TmGxy1;  
Lh TmGxx0, TmGyy0, TmGxy0;  
func u5= 1.;  
func u4=1.;  
func u3=1.;

TmGxx=0.;  
TmGyy=0.;

```
TmGxy=0.;

for(int j=0; j<100; ++j)
{  

```

solve bing([ux,uy],[vx,vy])  
=int2d(Th)((Mu)\*(Gradient(ux)‘\*Gradient(vx)  
+ Gradient(uy)’\*Gradient(vy)  
)  
)

```
+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
)
 -int2d(Th)([f1 , f2]' *[vx ,vy]
	)
	
	+ on(b1, ux=0., uy=0.)
			+ on(b2, ux=0., uy=0.)
            + on(b3, ux=u3, uy=0.)
			+ on(b4, ux=u4, uy=0.)
			+ on(b5, ux=u5, uy=0.);

Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
	real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
		 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
		 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
	real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
	real res = TaGnorm;
	if( 1.<res ){
		
		real q = 1./res;
		
		TmGxx[][k] = TaGxx*q;
		TmGyy[][k] = TaGyy*q;
		TmGxy[][k] = TaGxy*q;
	}
	else {
		
		TmGxx[][k] = TaGxx;
		TmGyy[][k] = TaGyy;
		TmGxy[][k] = TaGxy;
	}
}ux0[]=ux[];
uy0[]=uy[];

```

}  
fespace Rh(Th, P2);  
Rh dyg,dyg2;  
solve streamlines (dyg, dyg2,solver=UMFPACK)  
= int2d(Th)(  
dx(dyg)\*dx(dyg2)  
+ dy(dyg)_dy(dyg2)  
)  
+ int2d(Th)(  
- dyg2_(dy(ux) - dx(uy))  
)  
+ on(1,2,3,4,5,6, dyg=0);  
plot(dyg,ps=“ldsre1500.eps”,nbiso=10);

‘’’

Here is the uzawa bingham scheme  
laplace u^n+1+divergence lamda^n=f  
lamdaa^n+1=P(lamda^n+D(n^n+1))  
P(q)=q if |q|\<1 p(q)=q/|q| if |q|\>1

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 9, 2024, 12:20pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/4 "2024-05-09T12:20:53Z")

</div>

======Code=====================  
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;  
real L = 1; //Pipe length  
real D = 1.; //Pipe height  
real L1=.1;  
real L2=.9;  
real rho=1.; real alpha=1.;  
macro Gradient(u) [dx(u), dy(u)] //  
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //  
macro UgradV(ux,uy,vx,vy) [[ux,uy]'_[dx(vx),dy(vx)] , [ux,uy]'_[dx(vy),dy(vy)] ]//  
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;  
func f2=1.;  
border b1(t=0., 1.){x=t; y=0.;};  
border b2(t=0., 1.){x=L; y=D_t; };  
border b3(t=L, L2){x=t; y=D;};  
border b4(t=L2, L1){x=t; y=D;};  
border b5(t=L1, 0.){x=t; y=D;};  
border b6(t=0., 1.){x=0.; y=D-D_t;};  
int nnL = 50;  
int nnD = 50;  
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6( nnD));  
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);  
fespace Wh(Th,P2);  
fespace Qh(Th,P1);  
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;  
Wh ux,uy;  
Qh p, q, dp,p1;

fespace Lh(Th,P1);  
Lh TmGxx, TmGyy, TmGxy;  
Lh Txx, Tyy, Txy;  
Lh Gdotxx, Gdotyy, Gdotxy;  
Lh dTxx=0., dTyy=0., dTxy=0.;  
Lh TmGxx1, TmGyy1, TmGxy1;  
Lh TmGxx0, TmGyy0, TmGxy0;  
func u5= 1.;  
func u4=1.;  
func u3=1.;

TmGxx=0.;  
TmGyy=0.;

TmGxy=0.;

for(int j=0; j\<100; ++j)  
{  
solve bing([ux,uy],[vx,vy])  
=int2d(Th)((Mu)\*(Gradient(ux)'\*Gradient(vx)

- Gradient(uy)'\*Gradient(vy)  
)  
)

+int2d(Th)( Bn\*[TmGxx,TmGyy,TmGxy]‘\* StrainRate(vx,vy)  
)  
-int2d(Th)([f1 , f2]’ \*[vx ,vy]  
)

- on(b1, ux=0., uy=0.)
- on(b2, ux=0., uy=0.)
- on(b3, ux=u3, uy=0.)
- on(b4, ux=u4, uy=0.)
- on(b5, ux=u5, uy=0.);

Gdotxx = 2.\*dx(ux);  
Gdotyy = 2.\*dy(uy);  
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k\<Lh.ndof; ++k){  
real TaGxx=TmGxx[k]+r2_Gdotxx[][k],  
TaGyy=TmGyy[][k]+r2_Gdotyy[k],  
TaGxy=TmGxy[k]+r2_Gdotxy[][k];  
real TaGnorm = sqrt( .5_(TaGxx^2+TaGyy^2+2.\*TaGxy^2) );  
real res = TaGnorm;  
if( 1.\<res ){

real q = 1./res;

TmGxx[k] = TaGxx_q;  
TmGyy[x][k] = TaGyy_q;  
TmGxy[k] = TaGxy\*q;  
}  
else {

TmGxx[k] = TaGxx;  
TmGyy[k] = TaGyy;  
TmGxy[k] = TaGxy;  
}  
}ux0=ux;  
uy0=uy;

}  
fespace Rh(Th, P2);  
Rh dyg,dyg2;  
solve streamlines (dyg, dyg2,solver=UMFPACK)  
= int2d(Th)(  
dx(dyg)\*dx(dyg2)

- dy(dyg)\*dy(dyg2)  
)
- int2d(Th)(

- dyg2\*(dy(ux) - dx(uy))  
)

- on(1,2,3,4,5,6, dyg=0);  
================end of the code=====================

---

<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: [May 9, 2024, 1:10pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/5 "2024-05-09T13:10:40Z")

</div>

Sorry it appears still modified by the text interpreter. You have to use only back quotes (you can copy/paste this one ` if you are not sure of the character on your keyboard).  
Normally while composing your post you can see the result after the text interpreter in a separate right window, in which you can check what we will see.

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 9, 2024, 1:51pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/6 "2024-05-09T13:51:24Z")

</div>

I hope it can be seen.  
Thanks

```auto

real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;
real L = 1;	//Pipe length
real D = 1.;	//Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [[ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-D*t;};
int nnL = 50;
int nnD = 50;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6(nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy;
Qh p, q, dp,p1;

fespace Lh(Th,P1);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 1.;
func u4=1.;
func u3=1.;

	
TmGxx=0.; 
	TmGyy=0.;
	
	TmGxy=0.;
	
	for(int j=0; j<100; ++j)
	{  
solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu)*(Gradient(ux)'*Gradient(vx)
               + Gradient(uy)'*Gradient(vy)
			   )
			  )
		
	+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
	)
	 -int2d(Th)([f1 , f2]' *[vx ,vy]
		)
		
		+ on(b1, ux=0., uy=0.)
				+ on(b2, ux=0., uy=0.)
                + on(b3, ux=u3, uy=0.)
				+ on(b4, ux=u4, uy=0.)
				+ on(b5, ux=u5, uy=0.);

	Gdotxx = 2.*dx(ux);
	Gdotyy = 2.*dy(uy);
	Gdotxy = dx(uy)+dy(ux);

	for(int k=0; k<Lh.ndof; ++k){
		real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
			 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
			 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
		real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
		real res = TaGnorm;
		if( 1.<res ){
			
			real q = 1./res;
			
			TmGxx[][k] = TaGxx*q;
			TmGyy[][k] = TaGyy*q;
			TmGxy[][k] = TaGxy*q;
		}
		else {
			
			TmGxx[][k] = TaGxx;
			TmGyy[][k] = TaGyy;
			TmGxy[][k] = TaGxy;
		}
	}ux0[]=ux[];
	uy0[]=uy[];
	
	
	
}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + on(1,2,3,4,5,6, dyg=0);
	plot(dyg,ps="ldsre1500.eps",nbiso=10);

```

Here is the algoritm

 ![uzawa](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/c17e477a19e300cdfe3fb585c8ed83cea4adb5fe.png)

---

<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: [May 9, 2024, 2:13pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/7 "2024-05-09T14:13:53Z")

</div>

Great, it is clear now.

I find two minor problems:

1. The border b4 should be discretized as `b4(40)` instead of `b4(5)`

2. You strain rate is not normalized correctly, it should be  
`macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt(2.)]`  
in order for the scalar product  
`[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)`  
to give the correct value, and accordingly  
`Gdotxy = 2.*(dx(uy)+dy(ux))/sqrt(2.);`  
`real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+TaGxy^2) );`

I think that also you should take the strain rate space Lh as P1dc instead of P1, so that Grad(u) belongs to Lh when u is in P2.

But I think that the main problem is that your parameter `r2` is much too small. This implies that the convergence of `u` in the iteration is very bad (`u` appears in the strain rate update only behind the factor `r2`).  
You should take `r2` of order unity. But even in that case I am not sure that the iteration converges, in any case it could take thousands of iterations. You could try to add a relaxation term in (u^{k+1}-u^k) in the velocity equation (k is the iteration index), with a factor of order r2/h^2, see the PhD thesis of Duc Nguyen Hoai Theorem 4.2.6 p.71  
[https://theses.hal.science/tel-03078670/](https://theses.hal.science/tel-03078670/)

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 9, 2024, 6:55pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/8 "2024-05-09T18:55:37Z")

</div>

Thank you so much. I will correct the code what you have written.

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 10, 2024, 6:07am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/9 "2024-05-10T06:07:41Z")

</div>

Dear François Bouchut  
I did the things you wrote but the streamlines don’t change when the Bn number changes. I even added the relaxation term.

---

<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: [May 10, 2024, 10:56am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/10 "2024-05-10T10:56:27Z")

</div>

I think that the iteration is far from converged after your 100 iterations. You shoud check the difference ||u^{k+1}-u^k|| for each level k to have an idea. The method is known to be very slow. Later I will propose you a slightly different scheme with faster convergence.

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 10, 2024, 11:21am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/11 "2024-05-10T11:21:36Z")

</div>

> [@fb77](#):
>
> I think that the iteration is far from converged after your 100 iterations. You shoud check the difference ||u^{k+1}-u^k||. The method is known to be very slow. Later I will propose you a slightly different scheme with faster convergence.

Dear Professor I checked the difference it starts with 9e-5 ends as 1e-3 at the 100th iteration.

---

<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: [May 10, 2024, 12:10pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/12 "2024-05-10T12:10:38Z")

</div>

The first iteration is meaningless, but 1e-3 at the end looks not bad. Can you send me your code, I will look at it if I have time this afternoon.

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 10, 2024, 12:21pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/13 "2024-05-10T12:21:24Z")

</div>

Dear Professor ,  
Here is the latest code.

```auto
real Bn = 1.;real Mu=1.;

real r = 1.; 
real L = 1;	//Pipe length
real D = 1.;	//Pipe height
real L1=.1;
real L2=.9;
real rho=.1; real alpha=1.;real r2=10.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [[ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt(2.)] //

func f1=1e-1;
func f2=1e-1;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-D*t;};
int nnL = 32;
int nnD = 32;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(4)+b4(24)+b5(4)+b6(nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy,ux0,uy0;
Qh p, q, dp,p1;

fespace Lh(Th,P1dc);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 16.;
func u4=16.;
func u3=16.;

	
TmGxx=0.; 
	TmGyy=0.;
	
	TmGxy=0.;
	ux0=0.; uy0=0.;
	for(int j=0; j<100; ++j)
	{  
solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu+1./rho)*(Gradient(ux)'*Gradient(vx)
               + Gradient(uy)'*Gradient(vy)
			   )
			  )
		-int2d(Th)((1./rho)*(Gradient(ux0)'*Gradient(vx)
               + Gradient(uy0)'*Gradient(vy)
			   )
			  )  
	+int2d(Th)( [TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
	)
	 -int2d(Th)([f1 , f2]' *[vx ,vy]
		)
		
		+ on(b1, ux=0., uy=0.)
				+ on(b2, ux=0., uy=0.)
                + on(b3, ux=u3, uy=0.)
				+ on(b4, ux=u4, uy=0.)
				+ on(b5, ux=u5, uy=0.);

	Gdotxx = 2.*dx(ux);
	Gdotyy = 2.*dy(uy);
	Gdotxy = sqrt(2.)*(dx(uy)+dy(ux));

	for(int k=0; k<Lh.ndof; ++k){
		real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
			 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
			 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
		real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+TaGxy^2) );
		real res = TaGnorm;
		if( 1.<res ){
			
			real q = 1./res;
			
			TmGxx[][k] = TaGxx*q;
			TmGyy[][k] = TaGyy*q;
			TmGxy[][k] = TaGxy*q;
		}
		else {
			
			TmGxx[][k] = TaGxx;
			TmGyy[][k] = TaGyy;
			TmGxy[][k] = TaGxy;
		}
	}
	dux[] = ux[] - ux0[];
   duy[] = uy[] - uy0[];
   real err = sqrt(int2d(Th)(square(dux)+square(duy)));
   cout << err << "--iteration error /iteration step --" << j<< endl;
   ux0[]=ux[];
	uy0[]=uy[];
}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + on(1,2,3,4,5,6, dyg=0);
	plot(dyg,ps="ldsre1500.eps",nbiso=10);

```

---

<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: [May 10, 2024, 6:11pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/14 "2024-05-10T18:11:03Z")

</div>

Alice,  
The results obtained with your code look correct (but take care that in the last version you sent, you forget the factor Bn in front of `[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)`). Indeed for Bn=1. or Bn=10. it converges quite fast, this is because the solution is smooth, there is no plug zone (zone where strainrate=0).  
For Bn=100 it converges very slowly (and 100 iterations are not enough), this is because there is a plug zone.  
I remarked that you put Neumann BC on boundary b6. Is this intentional or is this that you forgot to put ux=uy=0?

Otherwise there is a way to have faster convergence. It works if you consider in your equation 2\mathop{\rm div}(\mathop{\rm StrainRate})(u) instead of \Delta u. It follows that you get  
`int2d(Th)(2.*Mu*(StrainRate(ux,uy)'*StrainRate(vx,vy)))`  
instead of `Grad'*Grad` previously.  
Then if you do not put the relaxation terms and if you choose the optimal coefficient `r2=Mu/Bn`, you get faster convergence.

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 11, 2024, 11:32am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/15 "2024-05-11T11:32:32Z")

</div>

Dear Professor,  
Thank you so much. I updated the code according to what you mentioned. Here is the latest code.  
The results for lid driven cavity squeezes to upward as bn increases as what is required.  
However I when bn is lower such as 1 and the force changes like in this paper the result does not match with the paper “Analysis of semi-staggered finite-difference method with application  
to Bingham flows ,Maxim A. Olshanskii”. Here mu=1,force=0.,u=1 when y=1.  
How can I fix?

```auto
real Bn = 1.; real Mu=1.;

real r = 1.; 
real L = 1;	//Pipe length
real D = 1.;	//Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;real r2=Mu/Bn;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [[ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt(2.)] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, 0.){x=t; y=D;};

border b4(t=0., 1.){x=0.; y=D-D*t;};
int nn = 50;

mesh Th=buildmesh(b1(nn)+b2(nn)+b3(nn)+b4(nn));
Th = adaptmesh(Th, 1./64., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy,ux0,uy0;
Qh p, q, dp,p1;

fespace Lh(Th,P1dc);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;

func u3=1.;

	
TmGxx=0.; 
	TmGyy=0.;
	
	TmGxy=0.;
	ux0=0.; uy0=0.;
	for(int j=0; j<10; ++j)
	{  
solve bing([ux,uy],[vx,vy])
=int2d(Th)(2.*Mu*(StrainRate(ux,uy)'*StrainRate(vx,vy)))

	+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
	)
	 -int2d(Th)([f1 , f2]' *[vx ,vy]
		)
		
		+ on(b1, ux=0., uy=0.)
				+ on(b2, ux=0., uy=0.)
                + on(b3, ux=u3, uy=0.)
				+ on(b4, ux=0., uy=0.);
				

	Gdotxx = 2.*dx(ux);
	Gdotyy = 2.*dy(uy);
	Gdotxy = sqrt(2.)*(dx(uy)+dy(ux));

	for(int k=0; k<Lh.ndof; ++k){
		real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
			 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
			 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
		real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+TaGxy^2) );
		real res = TaGnorm;
		if( 1.<res ){
			
			real q = 1./res;
			
			TmGxx[][k] = TaGxx*q;
			TmGyy[][k] = TaGyy*q;
			TmGxy[][k] = TaGxy*q;
		}
		else {
			
			TmGxx[][k] = TaGxx;
			TmGyy[][k] = TaGyy;
			TmGxy[][k] = TaGxy;
		}
	}
	dux[] = ux[] - ux0[];
   duy[] = uy[] - uy0[];
   
   real err = sqrt(int2d(Th)(square(dux)+square(duy)));
   cout << err << "--iteration error /iteration step --" << j<< endl;
   ux0[]=ux[];
	uy0[]=uy[];
}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + on(1,2,3,4,5,6, dyg=0);
	plot(dyg,ps="ldsre1500.eps",nbiso=10);

```

---

<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: [May 13, 2024, 1:15pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/16 "2024-05-13T13:15:19Z")

</div>

In that paper the flow is taken incompressible, whereas in your code you use compressible. This is why the solution is different (but looks correct).  
You can modify easily your code to do the incompressible case, by using the extra variables p,q as  
`solve bing([ux,uy,p],[vx,vy,q])`  
and adding  
`-int2d(Th)(p*Divergence(vx, vy)+q*Divergence(ux, uy))`  
But then it happens that the scheme converges extremely slowly (you need thousands of iterations). In order to improve it you should use the Augmented Lagrangian method as in that paper.

* * *

Another point is that your computation of the stream function could be misleading. You look for \psi satisfying  
-\Delta\psi=\partial\_y u\_x-\partial\_x u\_y.  
An important property that one would like that if \partial\_x u\_x+\partial\_y u\_y=0 then \partial\_x\psi=u\_y and -\partial\_y\psi=u\_x.  
You see then that the Dirichlet condition on \psi is not appropriate. It is better to take  
\frac{\partial\psi}{\partial n}=u\_y n\_x-u\_x n\_y on \partial\Omega.  
We are thus led to the problem  
-\Delta \psi=f in \Omega,\qquad\frac{\partial\psi}{\partial n}=g on \partial\Omega.  
To solve this system you need the compatibility condition (to see that it is necessary, integrate the equation over \Omega)  
\int\_\Omega f=-\int\_{\partial\Omega}g  
Assuming this, the system is not invertible since \psi is defined up to a constant. You can determine it by imposing \int\_{\partial\Omega}\psi=0.  
In order to get the solution you can then solve  
-\Delta \psi=f in \Omega,\qquad\frac{\partial\psi}{\partial n}+\varepsilon\psi=g on \partial\Omega,  
for some small \varepsilon. This system is well-posed. This leads to the code

```auto
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + int1d(Th)(
        - dyg2*(uy*N.x-ux*N.y)
    )
    + int1d(Th)(1e-8*dyg*dyg2)
    ;

```

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 14, 2024, 8:54am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/17 "2024-05-14T08:54:10Z")

</div>

Dear Professor,  
Thank you so much for your valuable precious support. When I changed the streamfunction as you wrote at the last message,the graph looks like this.

 ![ldsre1500](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/1/1f323aaeceaabd61a25796b9f6f694a8d7ee2daa.jpeg)  
Previously it looked like this  
 ![ldc_bing](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/c/cf06c2db3a2956995c7031fdd5165c0b3f5377ca.jpeg).What do you think that need to change?  
Regards

---

<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: [May 14, 2024, 10:32am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/18 "2024-05-14T10:32:57Z")

</div>

I think that the first (with non homogeneous Neumann BC) is the correct one.  
But as far as div(u) is non zero, this is not much meaningful.

---

<div class="post-metadata">

### Author: ![Alice](https://avatars.discourse-cdn.com/v4/letter/a/779978/32.png) [@Alice](https://community.freefem.org/u/Alice)
#### Post date: [May 14, 2024, 11:00am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/19 "2024-05-14T11:00:47Z")

</div>

Dear Professor ,  
Thank you so much. So for the incompressible case will the stream function algoritm change ?

---

<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: [May 14, 2024, 11:16am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/20 "2024-05-14T11:16:12Z")

</div>

The formulation with Neumann BC is especially good for the incompressible case, since as I said it is not much meaningfull when div u is not zero.

[Next page](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233.md?page=2)