// u_t + b1*div(u) = mu*Delta u + rho*g(u) with exact bc // travelling wave real mu= 1.0e-2, b1 = -1.0, rho = 100.0; real a = sqrt (rho / (4 * mu)); real b = 2 * b1 + sqrt (rho * mu); real c = a * (b - 1); 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,phi,delta,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 aL(u,phi) -mu*dx(u)*dx(phi)-mu*dy(u)*dy(phi)-(b1*dx(u)+b1*dy(u))*phi // macro g(u) rho * u ^ 2 * (1 - u) // macro dg(u) rho * (2 * u * (1 - u) - u ^ 2) // varf boundary(u,phi) = on(1,2,3,4,u=1); real[int] onboundary = boundary(0,Vh); //cout << onboundary << endl; varf F(dummy,phi) = int2d(Th)((u-u0)*phi -k*theta*(aL(u,phi)+g(u)*phi) -k*(1-theta)*(aL(u0,phi)+g(u0)*phi)); varf JF(delta,phi) = int2d(Th)(delta*phi -k*theta*(aL(delta,phi)+dg(u)*delta*phi)); problem ADR(delta,phi,solver = GMRES,eps = tol / 100) = JF + F +on(1,2,3,4,delta=0); for (int n = 0;n < ts;n++) { bc = 1.0 / (1 + exp (a * (x + y - b * (t + k)) + c)); u[] = onboundary ? bc[] : u0[]; // enforce bc //u = 16*x*(1-x)*y*(1-y)*u0+(1-16*x*(1-x)*y*(1-y))*bc; // better Newton cout << "time step: " << n << endl; ADR; conv = (delta[].linfty < tol); cout << "max delta: " << delta[].linfty << endl; while (!conv) { u = u + delta; ADR; conv = (delta[].linfty < tol); cout << "max delta: " << delta[].linfty << endl; } u = u + delta; // extra Newton iteration t = t + k; plot(u, wait = true, fill = true, value = true); u0 = u; } cout << "L2 error: " << sqrt (int2d(Th)((1.0 / (1 + exp (a * (x + y - b * t) + c)) - u) ^ 2)) << endl;