Newton solver did not converge

Question

Hello. I tried different coding for solving nonlinear variational problem for weeks, however, all do not work. Please help me clarify the problem. Briefly, I am working on forced convection flow at vertical plate. Thanks in advance. I really appreciate it. Here is my coding:

from __future__ import print_function
from fenics import *

#set parameter values
rho = 997.1 #density water
mu = 0.00089 #0.890*10^-3 #dynamic viscosity water
sigma = 0.05501 #electrical conductivity water
C = 4179 #heat capacitance water
c = 0.591 #thermal conductivity water
B_0 = 0

#create mesh
mesh = UnitSquareMesh(8,8)

#define function space
V = VectorElement('P', mesh.ufl_cell() , 2)
Q = FiniteElement('P', mesh.ufl_cell(), 2)
element = MixedElement([V, Q, Q])
W = FunctionSpace(mesh, element)

#define boundaries
walls = 'near(x[1],0) || near(x[1],1)'
inflow = 'near(x[0],0)'
outflow = 'near(x[0],1)'

#define boundary condition

#walls
bcu_walls = DirichletBC(W.sub(0), Constant((0,0)), walls)

#bcp_walls = DirichletBC(W.sub(1), Constant(0), walls) #pressure
bcT_walls = DirichletBC(W.sub(2), Constant(290), walls) #T_w

#inflow
bcT_inflow = DirichletBC(W.sub(2), Constant(295), inflow) 

#outflow
bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow) #p constant
bcT_outflow = DirichletBC(W.sub(2), Constant(292), outflow) #uniform temperature

bcu = [bcu_walls]
bcT = [bcT_walls, bcT_inflow, bcT_outflow] 
bcp = [bcp_outflow]

bcs = bcu + bcT + bcp

#define test functions
(v, q, s) = TestFunctions(W)

#define functions
u = Function(W)
p = Function(W)
T = Function(W)

#split functions
w = Function(W)
(u, p, T) = split(w)

#define expression used in variational forms
mu = Constant(mu)
rho = Constant(rho)
sigma = Constant(sigma)
C = Constant(C)
c = Constant(c)
B_0 = Constant(B_0)

#define variational problem
F1 = div(u)*q*dx 

F2 = rho*dot(dot(u,nabla_grad(u)),v)*dx + div(v)*p*dx + mu*inner(nabla_grad(u),nabla_grad(v))*dx \
   + sigma*pow(B_0,2)*dot(u,v)*dx 

F3 = rho*C*dot(u,grad(T))*s*dx-c*div(grad(T))*s*dx

F = F1 + F2 + F3

# Create VTK files for visualization output
vtkfile_u = File('steady_system/u.pvd')
vtkfile_p = File('steady_system/p.pvd')
vtkfile_T = File('steady_system/T.pvd')

# Solve nonlinear variational problem
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

# Save solution to file (VTK)
(u, p, T) = w.split()
vtkfile_u << u
vtkfile_p << p
vtkfile_T << T

# Plot solution
plot(u, title='Velocity')
plot(p, title='Pressure')
plot(T, title='Temperature')    

# Hold plot
interactive()

and here is the error:

File "steadyforced.py", line 114, in
solver.solve() RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

---

*** fenics-support@googlegroups.com

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0

---

*** DOLFIN version: 2016.2.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Aborted (core dumped)

Answer 1

This are just the steady, incompressible, isothermal Navier-Stokes equations, which don't depend on $T$, and the heat equation. Why don't you solve the flow problem and the heat transport (with $u$ given) separately?

Anyways, I assume the Newton method does not converge because your discretization of the flow problem is unstable. Try using P2 for the velocity and P1 for the pressure space, check the tutorial: https://fenicsproject.org/pub/tutorial/html/._ftut1009.html#ftut1:NS

Also, Newton's method has a small radius of convergence (your initial value needs to be "close" to the solution) which depends on the Reynolds number. You might need to perform some Picard iterations (or fixed point) first, in order to achieve convergence with Newton's method.

Hope this helps

Comments

I try to omit the temperature and try to solve the flow problem. This is my coding:

from __future__ import print_function
from fenics import *

#set parameter values
rho = 997.1 #density water
mu = 1 #0.890*10^-3 #dynamic viscosity water
sigma = 0.05501 #electrical conductivity water
C = 4179 #heat capacitance water
c = 0.591 #thermal conductivity water
B_0 = 0

#create mesh
mesh = UnitSquareMesh(8,8)

#define function space
V = VectorElement('P', mesh.ufl_cell() , 2)
Q = FiniteElement('P', mesh.ufl_cell(), 1)
element = MixedElement([V, Q])
W = FunctionSpace(mesh, element)


#define boundaries
walls = 'near(x[1],0)'
inflow = 'near(x[0],0)'
outflow = 'near(x[0],1)'

#define boundary condition

#walls
bcu_walls = DirichletBC(W.sub(0), Constant((0,0)), walls)

#inflow
#bcu_inflow = DirichletBC(W.sub(0), Constant((1,0)), inflow) #u_e=1
#bcp_inflow = DirichletBC(W.sub(1), Constant(2), inflow)

#outflow
bcp_outflow = DirichletBC(W.sub(1), Constant(0), outflow) #p constant

bcu = [bcu_walls]
bcp = [bcp_outflow]

bcs = bcu + bcp

#define test functions
(v, q) = TestFunctions(W)

#define functions
u = Function(W)
p = Function(W)

#split functions
w =Function(W)
(u, p) = split(w)


#define expression used in variational forms
n = FacetNormal(mesh)
mu = Constant(mu)
rho = Constant(rho)
sigma = Constant(sigma)
C = Constant(C)
c = Constant(c)
B_0 = Constant(B_0)

#define variational problem
F1 = div(u)*q*dx 

F2 = rho*dot(dot(u,nabla_grad(u)),v)*dx + div(v)*p*dx + mu*inner(nabla_grad(u),nabla_grad(v))*dx \
   + sigma*pow(B_0,2)*dot(u,v)*dx 

F = F1 + F2 

# Create VTK files for visualization output
vtkfile_u = File('steady_system/u.pvd')
vtkfile_p = File('steady_system/p.pvd')


# Solve variational problem for time step
J = derivative(F, w)
problem = NonlinearVariationalProblem(F, w, bcs, J)
solver  = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 25
prm['newton_solver']['relaxation_parameter'] = 1.0
solver.solve()

# Save solution to file (VTK)
(u, p) = w.split()
vtkfile_u << u
vtkfile_p << p

# Plot solution
plot(u, title='Velocity')
plot(p, title='Pressure')    

# Hold plot
interactive()

However, this occurs in terminal. I cant spot the mistakes. Thanks for your help.

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-08) r (rel) = -nan (tol = 1.000e-07)
  Newton solver finished in 0 iterations and 0 linear solver iterations.

---

In this last example I think that the solution is (0,). Try a non null rhs to debug your code...
Hope this helps

---

Sorry raraya, I didnt get what you meant. Could you explain? Thank you so much.

---

I mean that your system (last code) have a null solution (u,p) =(0,0)!! Probably the Newton scheme is dividing for the first residual which is zero, may be that explains "r (rel) = -nan"... Put a different condition (bc o rhs) to have a non zero solution to test your scheme