// u_t + b*div(u) = mu*Delta u + rho*g(u)
real mu = 1.0e-2, b1 = -1.0, b2 = -1.0, rho = 10.0;
real[int] viso=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0];
int N = 50;
mesh Th = square(N,N);
fespace Vh(Th,P1); int r=1;
Vh u0,u,v,delta,phi;
bool conv;
int ts = 128;
real t = 0.0,k = 1.0/ts,tol=(1.0/N)^(r+1)/100;
u0 = 16*x*(x-1)*y*(y-1);
real theta = 1.0; // don't use 0 (Newton useless)
macro L(u,v) -mu*dx(u)*dx(v)-mu*dy(u)*dy(v)-(b1*dx(u)+b2*dy(u))*v //
macro g(u) rho*u^2*(1-u) //
macro dg(u) rho*(2*u*(1-u)-u^2) //
varf JADR(delta,phi)=int2d(Th)(delta*phi-k*theta*(L(delta,phi)+dg(u)*delta*phi));
problem ADR(delta,phi,solver=UMFPACK) = JADR
  +int2d(Th)((u-u0)*phi-k*theta*(L(u,phi)+g(u)*phi)-k*(1-theta)*(L(u0,phi)+g(u0)*phi))+on(2,3,delta=0); //+on(1,2,3,4,delta=0)
for (int n=0;n<ts;n++)
  {
    u = u0;
    cout << n << endl;
    ADR;
    conv = (delta[].linfty < tol);
    while (!conv)
      {
        u = u+delta;
	    ADR;
	    cout << delta[].linfty << endl;
	    conv = (delta[].linfty < tol);
      }
    plot(u,wait=true,viso=viso,fill=true,value=true);
    u0 = u;
  }