I’m working on a 2D Navier-Stokes solver that I want to be able to use to model airflow past an airfoil but I’m wrestling with FEniCS a little bit. I’m trying to use a newton solver with a relaxation factor to help it converge but I don’t think this is working the way that I want it to. This code fails to converge even when parameters are chosen that should almost guarantee a stable solution. I’m still a FEniCS noob so I think there must be some stupid error that I’ve made that I’m hoping those with more experience can pick out easily.

```
#Test and Trial Space
w=Function(W)
w_old=Function(W)
u,p=split(w)
(u_old,p_old)=split(w_old)
# Define functions for solutions at previous and current time steps
(v,q)=TestFunctions(W)
# Define expressions used in variational forms
f = Constant((0, 0))
k = Constant(dt)
R=Constant(1./nu)
tao=Constant(R/dt)
# Define variational problem for step 1
F1= (inner(grad(u),grad(v)))*dx + tao*dot(u,v)*dx - div(v)*p*dx + R*inner(dot(u,grad(u)),v)*dx - tao*dot(u_old,v)*dx
# a(u_,v) b(v,p_) c(u_,u_,v)
F2 = div(u)*q*dx
F=F1+F2
# Time-stepping
t = 0
eps=10**-6
maxit=500
j=derivative(F,w)
problem=NonlinearVariationalProblem(F,w,bcs,J=j)
solver=NonlinearVariationalSolver(problem)
prm=solver.parameters
prm['newton_solver']['relaxation_parameter']=0.0001
n=0
for n in range(num_steps):
print("Elapsed time: {}".format(n*dt))
# Update current time
t += dt
n+=1
solver.solve()
# Update previous solution
w_old.assign(w)
(u_old,p_old)=w.split()
xdmffile_u.write(u_old,t)
xdmffile_p.write(p_old,t)```
```