// example33.edp, from chapt3/NSNewton.edp
// Original Author: F. Hecht  
// Jan 2012 Stationary incompressible Navier-Stokes Equation with Newton method.
// around a 3d Cylinder 

verbosity=0;

// build the  Mesh
real R = 5, L=15;
border cc(t=0, 2*pi){ x=cos(t)/2; y=sin(t)/2; label=1;}
border ce(t=pi/2, 3*pi/2) { x=cos(t)*R; y=sin(t)*R; label=1;}
border beb(tt=0,1) { real t=tt^1.2; x= t*L; y= -R; label = 1;}
border beu(tt=1,0) { real t=tt^1.2; x= t*L; y= R; label = 1;}
border beo(t=-R,R) {  x= L; y= t; label = 0;}
border bei(t=-R/4,R/4) {  x= L/2; y= t; label = 0;}

mesh Th=buildmesh(cc(-50) + ce(30) + beb(20) + beu(20) + beo(10) + bei(10));
plot(Th, wait=true);

//  macros
macro Grad(u1,u2) [ dx(u1),dy(u1), dx(u2),dy(u2)]// 
macro UgradV(u1,u2,v1,v2) [ [u1,u2]' * [dx(v1),dy(v1)] , 
                            [u1,u2]' * [dx(v2),dy(v2)] ]// 
macro div(u1,u2)  (dx(u1) + dy(u2))//

//  FE Spaces
fespace Xh(Th,P2);
fespace Mh(Th,P1);
Xh u1,u2,v1,v2,du1,du2,u1p,u2p;
Mh p,q,dp,pp;

//  Physical parameter 
real nu= 1./50, nufinal=1/200., cnu=0.5; 

// stop test for Newton 
real eps=1e-6;

// intial guess with B.C. 
u1 = ( x^2 + y^2 ) > 2;
u2=0;

while(true){  //  Loop on viscosity 
  int n;
  real err=0;	
  for( n=0; n< 15; n++){ // Newton Loop 
    solve Oseen([du1,du2,dp],[v1,v2,q]) =
      int2d(Th) (  nu*( Grad(du1, du2)' * Grad(v1, v2) )
                 + UgradV(du1, du2, u1, u2)' * [v1, v2]
                 + UgradV( u1, u2, du1, du2)' * [v1, v2]
                 - div(du1, du2)*q - div(v1, v2)*dp 
                 - 1e-8*dp*q // stabilization term 
                )
      - int2d(Th) (  nu*(Grad(u1, u2)' * Grad(v1, v2) )
                     + UgradV(u1,u2, u1, u2)' * [v1, v2]
                     - div(u1, u2)*q - div(v1, v2)*p 
                     - 1e-8*p*q 
                  )
      + on(1,du1=0,du2=0)  
      ;

    u1[] -= du1[];
    u2[] -= du2[];
    p[]  -= dp[];
    real Lu1=u1[].linfty,  Lu2 = u2[].linfty , Lp = p[].linfty;
    err= du1[].linfty/Lu1 + du2[].linfty/Lu2 + dp[].linfty/Lp;
	    
    cout << n << " err = " << err << " " << eps << " Re  =" << 1./nu << endl;
    if(err < eps) break; // converge 
    if( n>3 && err > 10.) break; //  Blowup ????           
  }

  if(err < eps){  // if converge  decrease nu (more difficult)
    plot([u1,u2], p, wait=1, cmm=" Re = " + 1./nu , coef=0.3);

    if( nu == nufinal) break; 

    if( n < 4) cnu=cnu^1.5; // fast converge => change faster 
    nu = max(nufinal, nu * cnu); // new viscosity 
    u1p=u1;  u2p=u2;  pp=p; //  save correct solution ...

  } else {  // if blowup,  increase nu (more simple)
    assert(cnu< 0.95); // final blowup ...  
    nu = nu/cnu; //  get previous value of viscosity 
    cnu= cnu^(1./1.5); // no conv. => change lower 
    nu = nu * cnu;  // new vicosity
    cout << " restart nu = " << nu << " Re= "<< 1./nu << 
            "  (cnu = " << cnu << " )" << endl;
    // restore correct solution ..
    u1=u1p;
    u2=u2p;
    p=pp; 
  }   
} // end while loop