// u_t + alpha*div(u) = epsilon*Delta u + gamma*g(u) with exact bc
real epsilon= 1.0e-2, alpha = -1.0, gamma = 100.0;
real a = sqrt(gamma/(4*epsilon));
real b = 2*alpha+sqrt(gamma*epsilon);
real c = a*(b-1);
// uexact = 1.0/(1+exp(a*(x+y-b*t)+c));
func u0fun = 1.0/(1+exp(a*(x+y-b*0)+c));
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,bc;
bool conv;
int ts = 64;
real t = 0.0,k = 1.0/ts,tol=(1.0/N)^(r+1)/100;
u0 = u0fun;
real theta = 1.0; // don't use 0 (Newton useless)
macro L(u,v) -epsilon*dx(u)*dx(v)-epsilon*dy(u)*dy(v)-(alpha*dx(u)+alpha*dy(u))*v //
macro g(u) gamma*u^2*(1-u) //
macro dg(u) gamma*(2*u*(1-u)-u^2) //
  varf boundary(u,v) = on(1,2,3,4,u=1);
real[int] onboundary = boundary(0,Vh);
varf JADR(delta,phi)=int2d(Th)(delta*phi-k*theta*(L(delta,phi)+dg(u)*delta*phi));
problem ADR(delta,phi,solver=GMRES,eps=tol/100) = JADR
  +int2d(Th)((u-u0)*phi-k*theta*(L(u,phi)+g(u)*phi)-k*(1-theta)*(L(u0,phi)+g(u0)*phi))+on(1,2,3,4,delta=0); //+on(1,2,3,4,delta=0)
for (int n=0;n<ts;n++)
  {
    u = u0;
    bc = 1.0/(1+exp(a*(x+y-b*(t+k))+c));
// enforce bc
    u[] = onboundary ? bc[] : u[];
    cout << n << endl;
    ADR;
    conv = (delta[].linfty < tol);
    while (!conv)
      {
    	u = u+delta;    
	    ADR;    
    	cout << delta[].linfty << endl;
    	conv = (delta[].linfty < tol);
      } 
    t = t+k;    
    plot(u,wait=true,fill=true,value=true);
    u0 = u;
  }
cout << sqrt(int2d(Th)((1.0/(1+exp(a*(x+y-b*t)+c))-u)^2)) << endl;