# 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:** 10
**Page:** 2

<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:24am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/21 "2024-05-14T11:24:57Z")

</div>

Thank you ,Professor.

---

<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 16, 2024, 3:33pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/22 "2024-05-16T15:33:04Z")

</div>

Sorry, about the incompressible version I forgot one line to determine the additive constant in the pressure  
`+int2d(Th)(-1.e-10*p*q)`  
With that the scheme converges rather nicely! Here is the full code (with other minor modifications).  
I think that Mu=sqrt(2.) corresponds to the test of Olshanskii.

```auto
verbosity=0;
real Bn = 1.; real Mu=sqrt(2.);

real L = 1; //Pipe length
real D = 1.; //Pipe height
real L1=.1;
real L2=.9;
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=0.;
func f2=0.;
border b1(t=0., 1.){x=t*L; 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 = 64;

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

Wh ux,uy,ux0,uy0,vx,vy;
Qh p,q;
Qh normsigma,normsigmaplot;

fespace Lh(Th,P1dc);
Lh TmGxx, TmGyy, TmGxy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh TmGxx0, TmGyy0, TmGxy0;
Lh Tauxx,Tauyy,Tauxy,Taunorm;

func u3=1.;

        TmGxx=0.; 
        TmGyy=0.;
        TmGxy=0.;
        TmGxx0=TmGxx;TmGyy0=TmGyy;TmGxy0=TmGxy;
        ux0=0.; uy0=0.;

        real errold=1.e10;
        real errstop = 0.1/8.*L/nn;
        for(int j=0; j<500; ++j){  
solve bing([ux,uy,p],[vx,vy,q],solver=UMFPACK)
=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])
        -int2d(Th)(Mu*p*Divergence(vx, vy)+Mu*q*Divergence(ux, uy))
        +int2d(Th)(-1.e-10*Mu*p*q)
        + 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));

        Tauxx=TmGxx+r2*Gdotxx;
        Tauyy=TmGyy+r2*Gdotyy;
        Tauxy=TmGxy+r2*Gdotxy;
        Taunorm=sqrt(.5*(Tauxx^2+Tauyy^2+Tauxy^2));
        TmGxx=Tauxx/max(1.,Taunorm);
        TmGyy=Tauyy/max(1.,Taunorm);
        TmGxy=Tauxy/max(1.,Taunorm);

   real normsigest=sqrt(int2d(Th)((TmGxx+2.*Mu/Bn*dx(ux))^2+(TmGyy+2.*Mu/Bn*dy(uy))^2
               +(TmGxy+2.*Mu/Bn*(dx(uy)+dy(ux))/sqrt(2.))^2));
   real err = sqrt(int2d(Th)(square(ux - ux0)+square(uy - uy0)));
   real errsigma = sqrt(int2d(Th)(square(TmGxx-TmGxx0)+square(TmGyy-TmGyy0)+square(TmGxy-TmGxy0)));
   cout << err << " " << errsigma << " --iteration error u,sigma /iteration step --" << j<< endl;
   if (r2*max(err,errold)<normsigest*errstop) break;
   ux0=ux;
   uy0=uy;
   TmGxx0=TmGxx;TmGyy0=TmGyy;TmGxy0=TmGxy;
   errold=err;
}// end of iteration loop on j

normsigma=sqrt((TmGxx+2.*Mu/Bn*dx(ux))^2+(TmGyy+2.*Mu/Bn*dy(uy))^2
               +(TmGxy+2.*Mu/Bn*(dx(uy)+dy(ux))/sqrt(2.))^2);
        plot(ux,wait=1,fill=1,value=true);
        plot([ux,uy],wait=1,nbiso=10,value=true);

normsigmaplot=max(normsigma,1.37);
real[int] valiso(15);
for (int i=0; i<valiso.n;i++)
    valiso[i]=1.37+i*0.016/3.;
real[int] colorhsv=[
4./6., 1 , 0.5,
4./6., 1 , 1,
5./6., 1 , 1,
1, 1. , 1,
1, 0.5 , 1
];
plot(normsigmaplot, ps = "normsigma_solution.eps",viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);

load "Element_P3"
fespace Rh(Th, P3);
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)
    ;
    real errdiv=int2d(Th)(abs(Divergence(ux, uy)));
    cout << "errdiv = " << errdiv << endl;
    real errstream=sqrt(int2d(Th)((dx(dyg)-uy)^2+(dy(dyg)+ux)^2));
    cout << "errstream = " << errstream << endl;
    plot(dyg,ps="ldsre1500.eps",nbiso=10,value=true);

```

---

<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 17, 2024, 6:48am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/23 "2024-05-17T06:48:53Z")

</div>

> [@fb77](#):
>
> ```auto
> verbosity=0;
> real Bn = 1.; real Mu=sqrt(2.);
> 
> real L = 1; //Pipe length
> real D = 1.; //Pipe height
> real L1=.1;
> real L2=.9;
> 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=0.;
> func f2=0.;
> border b1(t=0., 1.){x=t*L; 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 = 64;
> 
> mesh Th=buildmesh(b1(nn)+b2(nn)+b3(nn)+b4(nn));
> //Th = adaptmesh(Th, 1./64., IsMetric=1, nbvx=10000);
> plot(Th,wait=1);
> fespace Wh(Th,P2);
> fespace Qh(Th,P1);
> 
> Wh ux,uy,ux0,uy0,vx,vy;
> Qh p,q;
> Qh normsigma,normsigmaplot;
> 
> fespace Lh(Th,P1dc);
> Lh TmGxx, TmGyy, TmGxy;
> Lh Gdotxx, Gdotyy, Gdotxy;
> Lh TmGxx0, TmGyy0, TmGxy0;
> Lh Tauxx,Tauyy,Tauxy,Taunorm;
> 
> func u3=1.;
> 
> TmGxx=0.; 
> TmGyy=0.;
> TmGxy=0.;
> TmGxx0=TmGxx;TmGyy0=TmGyy;TmGxy0=TmGxy;
> ux0=0.; uy0=0.;
> 
> real errold=1.e10;
> real errstop = 0.1/8.*L/nn;
> for(int j=0; j<500; ++j){  
> solve bing([ux,uy,p],[vx,vy,q],solver=UMFPACK)
> =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])
> -int2d(Th)(Mu*p*Divergence(vx, vy)+Mu*q*Divergence(ux, uy))
> +int2d(Th)(-1.e-10*Mu*p*q)
> + 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));
> 
> Tauxx=TmGxx+r2*Gdotxx;
> Tauyy=TmGyy+r2*Gdotyy;
> Tauxy=TmGxy+r2*Gdotxy;
> Taunorm=sqrt(.5*(Tauxx^2+Tauyy^2+Tauxy^2));
> TmGxx=Tauxx/max(1.,Taunorm);
> TmGyy=Tauyy/max(1.,Taunorm);
> TmGxy=Tauxy/max(1.,Taunorm);
> 
> real normsigest=sqrt(int2d(Th)((TmGxx+2.*Mu/Bn*dx(ux))^2+(TmGyy+2.*Mu/Bn*dy(uy))^2
> +(TmGxy+2.*Mu/Bn*(dx(uy)+dy(ux))/sqrt(2.))^2));
> real err = sqrt(int2d(Th)(square(ux - ux0)+square(uy - uy0)));
> real errsigma = sqrt(int2d(Th)(square(TmGxx-TmGxx0)+square(TmGyy-TmGyy0)+square(TmGxy-TmGxy0)));
> cout << err << " " << errsigma << " --iteration error u,sigma /iteration step --" << j<< endl;
> if (r2*max(err,errold)<normsigest*errstop) break;
> ux0=ux;
> uy0=uy;
> TmGxx0=TmGxx;TmGyy0=TmGyy;TmGxy0=TmGxy;
> errold=err;
> }// end of iteration loop on j
> 
> normsigma=sqrt((TmGxx+2.*Mu/Bn*dx(ux))^2+(TmGyy+2.*Mu/Bn*dy(uy))^2
> +(TmGxy+2.*Mu/Bn*(dx(uy)+dy(ux))/sqrt(2.))^2);
> plot(ux,wait=1,fill=1,value=true);
> plot([ux,uy],wait=1,nbiso=10,value=true);
> 
> normsigmaplot=max(normsigma,1.37);
> real[int] valiso(15);
> for (int i=0; i<valiso.n;i++)
> valiso[i]=1.37+i*0.016/3.;
> real[int] colorhsv=[
> 4./6., 1 , 0.5,
> 4./6., 1 , 1,
> 5./6., 1 , 1,
> 1, 1. , 1,
> 1, 0.5 , 1
> ];
> plot(normsigmaplot, ps = "normsigma_solution.eps",viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);
> 
> load "Element_P3"
> fespace Rh(Th, P3);
> 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)
> ;
> real errdiv=int2d(Th)(abs(Divergence(ux, uy)));
> cout << "errdiv = " << errdiv << endl;
> real errstream=sqrt(int2d(Th)((dx(dyg)-uy)^2+(dy(dyg)+ux)^2));
> cout << "errstream = " << errstream << endl;
> plot(dyg,ps="ldsre1500.eps",nbiso=10,value=true);
> 
> ```

Dear Professor,  
Thank you so much for this code. Right now I am trying to run the code but it stops at the 3rd iteration like this

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/4/4724f1fc37d47726551e8de9a8cfbdafdc4d4601.png).  
It stops to compile.  
How can I fix?

---

<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 17, 2024, 8:40am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/24 "2024-05-17T08:40:03Z")

</div>

It looks like the iteration is just finished, just type “enter” in the plot window to see the diagnostics!

For convenience you can add a line  
`cout << "end of computations, type return in the plot window to plot " << endl;`  
just after the end of the iteration loop.

---

<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 17, 2024, 9:36am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/25 "2024-05-17T09:36:08Z")

</div>

Dear Professor,  
It looks great. Thank you so much.

---

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

</div>

Dear Professor,  
When I plot the streamlines ,I also want to draw rigid zones as in the examples;gray shaded areas. I couldn’t figure out how to include rigid zones into the streamlines plot. Could you give me any ideas please?

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/4/40e46a19da898fca6761f512ac3991d7f8aed10d.jpeg)

---

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

</div>

I think it’s not possible to draw two different type of plots on the same window in FF++. One thing you could do is post treatment with an image freeware in order to merge two output files.

But a try is to rescale one function to be in a range close to the other. For the version of the code I have posted, the plotting section becomes

```auto
normsigmaplot=max(normsigma,1.37);
real[int] valiso(15);
for (int i=0; i<valiso.n;i++)
    valiso[i]=1.37-0.0746666+2*i*0.016/3.;// doubled range
real[int] colorhsv=[
4./6., 1 , 0.5,
4./6., 1 , 1,
5./6., 1 , 1,
1, 1. , 1,
1, 0.5 , 1
];
plot(normsigmaplot, ps = "normsigma_solution.eps",viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);

load "Element_P3"
fespace Rh(Th, P3);
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)
    ;
    real errdiv=int2d(Th)(abs(Divergence(ux, uy)));
    cout << "errdiv = " << errdiv << endl;
    real errstream=sqrt(int2d(Th)((dx(dyg)-uy)^2+(dy(dyg)+ux)^2));
    cout << "errstream = " << errstream << endl;
    plot(dyg,ps="ldsre1500.eps",nbiso=10,value=true,wait=1);

dyg=min(dyg,0.07);// to avoid larger values
dyg=1.37-0.0746666+(dyg+0.005)/0.093*14*0.016/3.;//rescaling
plot(dyg,normsigmaplot,viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);

```

Here the values \<1.37 are for the stream function, and the values \>1.37 are for \|\sigma\| (values between 1.37 and 1.44 because we focus around the value \sqrt{2}).

 ![normsigma-stream_solution](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/2X/9/93ab4a2caa1f5e474ed1f94e6747612bd93463c7.jpeg)

---

<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 30, 2024, 1:54pm UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/28 "2024-05-30T13:54:42Z")

</div>

Dear Professor ,  
Thank you so much.

---

<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: [June 3, 2024, 10:11am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/30 "2024-06-03T10:11:17Z")

</div>

Your result looks correct. You can take for plot commands

```auto
plot(p, cmm="Pressure",ps="pressurepoiseuille.eps",wait=1,value=true);
plot([ux, uy], cmm="Velocity",ps="velpoiseuille.eps",value=true);

```

Please open a new topic for a new subject (your question is not related to the subject “Uzawa scheme for stokes Bingham flow”). To open a new topic press the “+” in the bottom right of the FF window.

---

<div class="post-metadata">

### Author: ![ozil101](https://yyz2.discourse-cdn.com/flex030/user_avatar/community.freefem.org/ozil101/32/1820_2.png) [@ozil101](https://community.freefem.org/u/ozil101)
#### Post date: [June 3, 2024, 10:13am UTC](https://community.freefem.org/t/uzawa-scheme-for-stokes-bingham-flow/3233/31 "2024-06-03T10:13:48Z")

</div>

Thank you for your kind reply sir.

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