Newton solver does not converge due to non-linear term

1

Hello,

I am very new to FEniCS and trying to solve the following time-dependent, non-linear, 1D equation:
\frac{\partial h}{\partial t} + \frac{\partial}{\partial x}\left(h^3\frac{\partial^3h}{\partial x^3} + \frac 1h \frac{\partial h}{\partial x} \right)=0
with boundary conditions:
\frac{\partial h}{\partial x}=\frac{\partial^3 h}{\partial x^3}=0 ~~\text{at}~~ x=0,\frac{\lambda}{2}
and initial condition:
h(x,t=0)=\frac{10}{9} \left(1-\frac{1}{10}\cos(\frac{2\pi x}{\lambda})\right)

Introducing p=\frac{\partial^2 h}{\partial x^2} it is straightforward to obtain the variational formulation, with homogeneous Neumann boundary conditions for h and p. For now I only use a simple forward Euler scheme for the time.

I can’t do anything with that equation in FEniCS because of the term \frac{1}{h} \frac{\partial h}{\partial x}: at the first time-iteration, the newton algorithm does not converge and in fact the residual is already NaN at iteration 0. Everything seems fine when I remove it, or even replace the 1/h with some other functions of h.

I do not understand what could cause this. I do expect 1/h to be a issue eventually and lead to a singularity: h(x=0,t) \rightarrow 0 and I will have to deal with that later, but it should be OK at early times since the initial condition is nowhere close to 0.

I searched a bit in this forum but did not find tricks that could work, except using a random initial guess: I do not understand why but this seemingly improves things a bit: the Newton solver doesn’t converge, but at least the residual starts with a reasonable, non-NaN value (I tried to adjust the relaxation parameter, tolerances, … but couldn’t get it to converge).

Thanks for any help. My code is below

Note: This is a problem taken from the following reference
Zhang, W.W. and Lister, J.R., 1999. Similarity solutions for van der Waals rupture of a thin film on a solid substrate. Physics of Fluids , 11 (9), pp.2454-2462.

from fenics import *
import numpy as np
from numpy.random import rand

t_max = 20
DT = 0.001
dt = Constant(DT)

lambd=10
h0=Expression('10./9. * (1-0.1*cos(2*3.1415926535*x[0]/L))',degree=1,L=lambd)

Xmax=lambd/2.
ddx = 0.01

# mesh
mesh = IntervalMesh(int(Xmax/ddx), 0, Xmax)

# basis functions
element = FiniteElement('Lagrange', interval, 1) #Lagrange polynomial of degree 1?

# mixed elements and associated function space
mixed_elements = MixedElement([element, element])
V = FunctionSpace(mesh, mixed_elements)

# Variables
hp = Function(V)
h, p = split(hp)
h_old = interpolate(h0, V.sub(0).collapse())

# Trial functions
u, v = split(TestFunction(V))

# Weak form of the equation to solve
eq1 = p*v*dx + inner(h.dx(0), v.dx(0))*dx
eq2 = h*u*dx - h_old*u*dx -dt*h**3*inner(p.dx(0),u.dx(0))*dx-dt/h*inner(h.dx(0),u.dx(0))*dx
eq = eq1+eq2


t=0
while t<t_max:
    t+=DT
    J = derivative(eq,hp)
    problem = NonlinearVariationalProblem(eq,hp,J=J)
    solver  = NonlinearVariationalSolver(problem)
    prm = solver.parameters
    #prm["newton_solver"]["absolute_tolerance"] = 1E-6
    #prm["newton_solver"]["relative_tolerance"] = 1E-6
    #prm["newton_solver"]["maximum_iterations"] = 10000
    #prm["newton_solver"]["relaxation_parameter"] = 0.01
    
    # RANDOM INITIAL GUESS
    #hp.vector().set_local(rand(hp.vector().size()))
    #hp.vector().apply  
    
    solver.solve()
    
    height, height2 = hp.split(deepcopy=True)
    h_old.assign(height)
2

In your problem, if you do not start with a non-zero hp, the right hand side of your newton problem (Where F=eq1+eq2)
dF/d(hp)\Delta hp=-F(hp,uv)
Will be infinity, since eq(hp)=infinity.

I do not currently have the time to run your example. If I find time for it, i will keep you posted.

3

I just searched more info taking into account your comment and apparently the default initial guess of the Newton’s solver is zero. I somehow overlooked that fact before so thanks, this is already quite helpful.

So I’m assuming i should either find a trick to remove the singularity in the equation, come up with a smarter initial guess, or use something else than Newton’s method. I’ll look into it.

edit: just taking another initial guess seems to be acceptable already, it was quite stupid… Good learning experience I guess. I’ll still need to deal with the singularity but that’s another issue

1 Like