// u_t + alpha*div(u) = epsilon*Delta u + gamma*g(u) with exact bc // travelling wave 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 = 32; real t = 0.0,k = 1.0/ts,tol=(1.0/N)^(r+1)/100; u0 = u0fun; real theta = .5; // 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 F(dummy,phi) = int2d(Th)((u-u0)*phi-k*theta*(L(u,phi)+g(u)*phi)-k*(1-theta)*(L(u0,phi)+g(u0)*phi)); 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 + F +on(1,2,3,4,delta=0); //+on(1,2,3,4,delta=0) for (int n=0;n