Issue with the convergence of the Newton method
I a trying to solve a phase change problem governed by the incompressible Navier Stokes equations and an enthalpic energy equation. This problem couples temperature and velocity and it is a non linear one. Thus, I use a Newton solver to get the solution. However, it appears that the Newton solver begins to hang after several iterations : it does not diverge but the residual error oscillates between 10^-4 and 10^-6 whereas my tolerance parameter is at 10^-8. Does someone know what could be the origin of such a problem ?
Here is my code :
from __future__ import print_function from fenics import * from dolfin import * import numpy as np # Define constants and known data Tcold_ = -1. Thot_ = 1. Tcold = Constant(Tcold_) Thot = Constant(Thot_) Pr = Constant(1.0) # Prandtl number Ste = Constant(1.0) # Stefan number Ra = Constant(10**6) # Rayleigh number # Time discretization t = 0.0 T = 0.05 # final time num_steps = 20 # number of time steps dt = (T-t) / num_steps # time step size # Create mesh and define function space N = 50 # Characteristic number of elements mesh = RectangleMesh(Point(0,0),Point(1,1),N,N) # Function space V = VectorElement('P', mesh.ufl_cell(), 2) # velocity Q = FiniteElement('P', mesh.ufl_cell(), 1) # pressure element = MixedElement(V, Q, Q) # velocity, pressure, temperature in this order W = FunctionSpace(mesh, element) # Define boundaries lid = 'near(x, 1.0)' noslip = 'near(x, 0.0) || near(x, 0.0) || near(x, 1.0)' coldwall = 'near(x, 1.0)' hotwall = 'near(x, 0.0)' # Boundary conditions bc_noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), noslip) bc_lid = DirichletBC(W.sub(0), Constant((0.0, 0.0)), lid) bc_hot = DirichletBC(W.sub(2), Thot, hotwall) bc_cold = DirichletBC(W.sub(2), Tcold, coldwall) bc = [bc_noslip, bc_lid, bc_hot, bc_cold] # Variational formulation w = TestFunction(W) (v, q, theta) = split(w) # w.split() does not work here # Create and assign the initial conditions theta_init = 'Th+(Tc-Th)*x'.replace('Tc', str(Tcold_)).replace('Th', str(Thot_)) # theta_init = 'x<=0 ? Th : Tc'.replace('Tc', str(Tcold_)).replace('Th', str(Thot_)) # u_init and p_init are null w_n = interpolate(Expression(("0.", "0.", "0.", theta_init), element=element), W) u_n, p_n, theta_n = w_n.split() """ When you do call interpolate for u_n for example, this returns a new object u_n which replaces the previous one (which was previously returned by w_n.split()) in the namespace. Now later when you update w_n, the updates values aren't actually associated with your variational form. """ # Create the current solution w_ = Function(W) (u_, p_, theta_) = split(w_) # Define symmetric gradient def epsilon(u): return sym(nabla_grad(u)) # Define the viscosity and the phase change curve def phi(temp): Tr = Constant(0.0) r = Constant(0.1) return(1+tanh((Tr-temp)/r)/2) def mu(temp): muL = Constant(1.0) muS = Constant(10**5) return (muL + (muS - muL)*phi(temp)) # Define Boussinesq coupling term fB = Ra/Pr*theta_*Constant((0.0,1.0)) # Weak formulation F == 0 # Equations terms Eq1 = dot(u_,v)/dt*dx - dot(u_n,v)/dt*dx \ + dot(dot(u_, nabla_grad(u_)), v)*dx \ + 2*inner(mu(theta_)*epsilon(u_), grad(v))*dx \ - div(v)*p_*dx - dot(fB,v)*dx Eq2 = q*div(u_)*dx Eq3 = theta_*theta/dt*dx - theta_n*theta/dt*dx \ - 1/Ste*(phi(theta_)-phi(theta_n))*theta*dx \ + 1/Pr*dot(nabla_grad(theta_),nabla_grad(theta))*dx \ - dot(theta_*u_,nabla_grad(theta))*dx # Sum of all contributions F = Eq1 + Eq2 + Eq3 # Create the progress bar progress = Progress('Time-stepping') set_log_level(PROGRESS) # Provide the Jacobian form J = derivative(F,w_) # Create the Newton solver problem = NonlinearVariationalProblem(F, w_, bc, J) solver = NonlinearVariationalSolver(problem) prm = solver.parameters.newton_solver prm.error_on_nonconvergence = False prm = solver.parameters prm['newton_solver']['absolute_tolerance'] = 1E-8 prm['newton_solver']['relative_tolerance'] = 1E-7 prm['newton_solver']['relaxation_parameter'] = 0.8 # Save solution to file in VTK format vtkfile1 = File('Sharp transition/Velocity/velocity.pvd') vtkfile2 = File('Sharp transition/Pressure/pressure.pvd') vtkfile3 = File('Sharp transition/Temperature/temperature.pvd') # Form and solve the linear system for n in range(num_steps): # Update current time t += dt # Solve the variational problem solver.solve() # for use with the NonlinearVariationalProblem structure (u_, p_, theta_) = w_.split() # split(w_) doesn't work # Save data vtkfile1 << (u_, t) vtkfile2 << (p_, t) vtkfile3 << (theta_, t) # Update the previous solution w_n.vector()[:] = w_.vector() # Update the progress bar progress.update(t/T)
I'm mostly copying and pasting my comment from this discussion, because I think it answers your question:
This fully implicit single-domain enthalpy method approach to solving convection-coupled phase-change has some robustness issues, mostly related to convergence of Newton's method, that I discuss in my proceedings paper. I spent the better part of the last year figuring out how to implement this approach in FEniCS. My lessons learned are mostly captured in that paper and in the Phaseflow repository; so I suggest reviewing these materials.
But also more generally:
- The convergence behavior of Newton's method is very much problem dependent. You should try changing numerical parameters such as the grid refinement and the time step size.
- Relaxing Newton's method often helps. Right now you're using "Full Newton".