Issue with the convergence of the Newton method

3 months ago by
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], 1.0)'
noslip = 'near(x[1], 0.0) || near(x[0], 0.0) || near(x[0], 1.0)'
coldwall = 'near(x[0], 1.0)'
hotwall = 'near(x[0], 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[0]'.replace('Tc', str(Tcold_)).replace('Th', str(Thot_))
# theta_init = 'x[0]<=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)

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')

# 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
Community: FEniCS Project
Since you're not actually using any of these stabilization terms, would you clean up this working example, making it more minimal? There's a bunch of noise here not related to the question, mostly just the "Stabilization terms" section which can be deleted.
written 3 months ago by Alexander G. Zimmerman  

1 Answer

3 months ago by

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".
Please login to add an answer/comment or follow this question.

Similar posts:
Search »