Newton Solver does not change guess

Hello,
I have run into a strange error, or at least I think its strange.
For some reason the output of the Newton solver does not change and therefore it cant converge.
The output always remains:

Newton iteration 0: r (abs) = 1.161e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

It remains this for the whole 50 tries and even when rerunning it the numbers do not change, even when I give a random number as an initial guess the numbers remain the same.

from dolfin import *
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

tol=1E-12

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 2)

# ft01_poisson.py Example
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Laplace=plot(u, title='show u')
plt.colorbar(Laplace)
#plt.show()

######################################################
ho0 = Function(V)
s=TestFunction(V)
s_lm = TestFunction(V)

F0 = dot((-1*(ho0**2)/0.001),s)* dx + dot((3.5-u),s_lm)*dx

ho0J=derivative(F0,ho0)

Ho0problem = NonlinearVariationalProblem(F0, ho0, None, ho0J)
Ho0solver = NonlinearVariationalSolver(Ho0problem)

Ho0prm = Ho0solver.parameters
Ho0prm['newton_solver']['maximum_iterations'] = 5 
Ho0solver.solve()

As you can see the first part of the code is build on the ft01 Poisson example code but it only functions as an input for the second part where its output u is supposed to be part of a lagrange multiplier in F0 which is the formula that has been giving me and the group I’m working with quiet a few problems.
(which is why you might have seen it somewhere already. I’m trying a different approach and find the result encouraging but also strange)

Anyways if anybody got any insight into why this is happening and want to share, Thanks!

@Mister_Worst care to share the part where you actually use the Newton solver ? There could be a subtlety there…

The last line in the code is what uses should start the solver.
After that there is just the standard plot function for the result.
But that part I haven’t reached with this code yet, since I mostly either get error messages or the solver not converging.

ho0plot=plot(ho0, title='show ho0')
plt.colorbar(ho0plot)
plt.show()
print ('ho0 calculated')

Interestingly enough I ran the same code I had yesterday giving me the above output in the consol again today and it changed to being

Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)

after the first iteration instead, with no changes in the code

Of course sorry I hadn’t realised your code block was so large I had to scroll to get to the end…

I’m confused with your F0. If I was a Newton Solver and I was fed that, I think I would set ho0 to zero and then end up with a finite residual equal to 3.5-u. Are you sure that using different test functions for the first and second term is the way to go ?

no problem.
As far as I have seen this or similar has been the way a Lagrange Multiplier has been implemented in different demo/tutorial codes.
For example the Poisson Equation with pure Neumann Border Conditions
the thing is that they normally use a mixed Functions Space for this and then have s,s_lm=TrialFunctions(W) instead of two seperated commands.
But I havent been able to get that setup to work for the derivate function, which is needed for the NonLinearVariationalSolver as far as I know.
But yes I’m aware that this part is a bit iffy so I have been working on that as well.

I think you need to supply a test function to derivative in case of ambiguity. Something like : ho0J=derivative(F0,ho0,s).

I also agree with you that creating two TestFunction on the same space is extremely fishy. I think you should be able to go back to the MixedElement at little cost. Consider something like :

P2 = FunctionSpace(mesh, 'Lagrange', 2)
R = FunctionSpace(mesh, 'Real', 0)
V=P2*R

u,_=TrialFunction(V)
v,s = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u,_ = Function(V)
solve(a == L, u, bc)

Laplace=plot(u, title='show u')
plt.colorbar(Laplace)
#plt.show()

######################################################
ho0,_ = Function(V)

F0 = dot(-ho0**2/0.001,v)* dx + dot(3.5-u,s)*dx

ho0J=derivative(F0,ho0,v)

Also be advised that supplying ho0J to the Newton Solver is optional. I’m not 100% sure it’s necessary to declare u twice either.