Uzawa scheme for stokes bingham flow

Thank you ,Professor.

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.

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

.
It stops to compile.
How can I fix?

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.

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