Newton

load "gmsh"
int Verbosity = verbosity;
verbosity = 0;
mesh Th = gmshload("Mesh/teslaSmall02.msh");

bool backward = 1;

int INLET = 2-backward;
int OUTLET = 1+backward;
int WALL = 3;

real Remax = 500;
int steps = 10;
real k = 6e3;
int maxIters = 50;
real tol = 1e-8;
real Re = Remax/steps;

fespace Vh(Th, P2);
Vh u1, u2, v1, v2, du1, du2, w;

fespace Ph(Th, P1);
Ph p, q, dp;

macro grad(u) [dx(u),dy(u)]//
macro Grad(u1, u2) [grad(u1), grad(u2)]//
macro div(u1, u2) (dx(u1)+dy(u2))//
macro S(u1, u2) (0.5*(Grad(u1,u2)'+Grad(u1,u2)))//
macro coord(x) ((x+1.333018)*sqrt(2))//

func gbackward = k*coord(x)*(coord(x)-0.1);
func gforward = k*y*(y+0.1);
solve Stokes(u1,u2,p,v1,v2,q) =
	int2d(Th) (
		q*div(u1,u2)
		-p*div(v1,v2)
		+2*(S(u1,u2):S(v1,v2))
	)
	+on(INLET, u1=(backward ? gbackward : gforward)*N.x, u2=(backward ? gbackward : gforward)*N.y)
	+on(WALL, u1=0, u2=0)
;

plot(p, fill=1, value=1, wait=1);

problem NewtonNS(du1,du2,dp,v1,v2,q) =
	int2d(Th) (
		q*div(du1,du2)
		+Re*([du1, du2]'*Grad(u1, u2))*[v1, v2]
		+Re*([u1, u2]'*Grad(du1, du2))*[v1, v2]
		-dp*div(v1,v2)
		+2*(S(du1,du2):S(v1,v2))
		+dp*q*1e-8
	)
	+int2d(Th) (
		q*div(u1,u2)
		+Re*([u1,u2]'*Grad(u1,u2))*[v1,v2]
		-p*div(v1,v2)
		+2*(S(u1,u2):S(v1,v2))
	)
	+on(INLET, OUTLET, WALL, du1=0, du2=0)
;

for (int i = 1; i <= maxIters; i++) {
	if (Re > Remax) break;
	NewtonNS;
	real err = int2d(Th) (sqrt(du1^2+du2^2));
	u1 = u1+du1;
	u2 = u2+du2;
	p = p+dp;
	real norm = int2d(Th) (sqrt(u1^2+u2^2));
	cout << i << ": " << err/norm << endl;
	if (err/norm < tol) {
		Vh w;
		w = sqrt(u1^2+u2^2);
//		plot(w, fill=1, value=1, cmm=Re);
		plot(p, fill=1, value=1, cmm=Re);
		cout << "SOLUTION REYNOLDS " << Re << endl;
		real inflow = int1d(Th, INLET) (u1*N.x+u2*N.y);
		real outflow = int1d(Th, OUTLET) (u1*N.x+u2*N.y);
		cout << "Inflow: " << inflow << endl;
		cout << "Outflow: " << outflow << endl;
		real pin = int1d(Th, INLET) (p);
		real pout = int1d(Th, OUTLET) (p);
		cout << "PRESSURE DROP: " << pin-pout << endl;
		Re+=Remax/steps;
		i = 0;
	}
}