Newton Solver Convergence of a non linear problem for a mixed formulation Lagrange multiplier in 1d

I am new in fenicsx and I’m trying to solve in 1 dimension a non-linear problem with a mixed formulation in the context of a Lagrange multiplier.

However, the Newton Solver does not converge…

Would you have some advices please ?
Just in case, here is the code with some comments

# Import necessary packages
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot, nls
from dolfinx.mesh import locate_entities_boundary, meshtags
from ufl import dx, ds, grad
from mpi4py import MPI

# --- Mesh and function space setup ---
from basix.ufl import element,mixed_element

# Create a 1D mesh over the interval [0, 1] with 50 elements
msh = mesh.create_interval(comm=MPI.COMM_WORLD, points=(0.0, 1.0), nx=100)

# Define the function space V for r(z), the shape function of the bridge
# Using first-order (linear) Lagrange elements
PV=element("Lagrange",  "interval", 1)
PQ=element("DG",  "interval", 0)
V = fem.functionspace(msh, PV)
# Function space for lambda (scalar field — constant function over domain)
Q = fem.functionspace(msh, PQ)  # Discontinuous Galerkin of degree 0

#Mixed function space
W = fem.functionspace(msh, mixed_element([PV, PQ]))


# --- Define trial/test functions and coordinates ---

w=fem.Function(W)
dw=ufl.TrialFunction(W)
# Create the function r(z), representing the profile to optimize
r, lambda_ = ufl.split(w)#fem.Function(V),fem.Function(Q)

# Define the test function v(z) used in the weak formulation
v,mu = ufl.TestFunctions(W)#ufl.TestFunction(V), ufl.TestFunction(Q)

# Get the spatial coordinate (z) in the mesh
x = ufl.SpatialCoordinate(msh)

dr,dlambda = ufl.TrialFunctions(W)#ufl.TrialFunction(V)  ,ufl.TrialFunction(Q)     # Used to define Jacobian

# -------- initial condition setup---------------

w.sub(0).interpolate(lambda x: np.sqrt(V_target/np.pi)*np.ones(x.shape[0]))
w.sub(1).interpolate(lambda x: 1.0*np.ones(x.shape[0]))  

# --- Boundary condition setup ---

# Determine topological dimension of facets (0D in 1D mesh)
fdim = msh.topology.dim - 1

# Locate the mesh boundary facets at z = 0 and z = 1
facets = locate_entities_boundary(msh, fdim, 
                                  lambda x: np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0))

# Assign tags: 1 for z=0 (left), 2 for z=1 (right)
facet_values = np.array([1 if np.isclose(x[0], 0.0) else 2 for x in msh.geometry.x[facets]])

# Create a meshtag object associating facet indices to their boundary tags
facet_tags = meshtags(msh, fdim, facets, facet_values)

# Define a boundary integral measure ds that uses the tags
ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tags)

# --- Define physical parameters ---

# Adhesion coefficient at the contact lines
alpha = 1e-6

# Target volume of the bridge (not enforced directly here)
V_target = 1

#--- Define energy functional and weak form ---

# Compute r'(z), the derivative of r with respect to z
r_prime = grad(r)[0]
   
# Common denominator appearing in area and energy terms
denom = ufl.sqrt(1 + r_prime**2)

# Bulk energy integrand: surface area minus Lagrange multiplier term
bulk_term = (r * denom - 2 * lambda_ * np.pi * r) * v * dx    

# Variation of the area functional (arising from integration by parts)
variation_term = -(r * r_prime / denom) * grad(v)[0] * dx     

# Boundary term at z = 0 (tag 1), includes adhesion and curvature
bc_term_0 = (2 * alpha * r - r * r_prime / denom) * v * ds(1)

# Boundary term at z = H = 1 (tag 2), includes adhesion and curvature
bc_term_H = (2 * alpha * r + r * r_prime / denom) * v * ds(2)

# Volume constraint equation: residual in lambda direction
volume_constraint = -(np.pi * r**2 - V_target) * mu * dx

# Combine all contributions into the full residual (first variation)
F = bulk_term + variation_term + volume_constraint+ bc_term_0 + bc_term_H 

#-- Nonlinear problem setup and solution ---

# Compute the Jacobian of F with respect to (r, lambda)
J = ufl.derivative(F, w, dw)

# ----------------Define the nonlinear problem and solver---------------------------------------

problem = fem.petsc.NonlinearProblem(F, w, J=J)
solver = nls.petsc.NewtonSolver(msh.comm, problem)

# Set solver options
solver.rtol = 1e-5
solver.report = True
solver.line_search = True  # Enable line search
solver.max_it = 500  # Increase the maximum iterations

# Solve
n, converged = solver.solve(w)

assert converged

Heads up, these arguments does not do anything.
Secondly, what is the output of running this code with
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO) ?

Thank you for your answer.
I corrected one error that was to write x.shape[0] instead of x.shape[1]

However, it is still not converging.
Writing dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO) leads to:

[2025-05-08 01:00:28.387] [info] Requesting connectivity (0, 0) - (1, 0) ... [2025-05-08 01:00:28.391] [info] Column ghost size increased from 0 to 0 [2025-05-08 01:00:28.392] [info] Newton iteration 0: r (abs) = 0.020917390938103513 (tol = 1e-10), r (rel) = inf (tol = 0.0001) [2025-05-08 01:00:28.392] [info] PETSc Krylov solver starting to solve system. .... [2025-05-08 01:00:28.393] [info] Newton iteration 5: r (abs) = -nan (tol = 1e-10), r (rel) = -nan (tol = 0.0001)

Is that what you mean ?
From what it seems, at step 0, the absolute residual norm is ok, but the relative residual is not ?