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

**URL:** <https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559>\
**Category:** General\
**Tags:** boundary-conditions, dolfinx\
**Created:** [May 5, 2025, 4:12pm UTC](https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559 "2025-05-05T16:12:31Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jo\_A](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jo_a/32/8444_2.png) [@Jo\_A](https://fenicsproject.discourse.group/u/Jo_A)\
**Post date:** [May 5, 2025, 4:12pm UTC](https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559/1 "2025-05-05T16:12:31Z")

</div>

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

```auto
# 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

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 6, 2025, 2:18pm UTC](https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559/2 "2025-05-06T14:18:00Z")

</div>

> [@Jo\_A](#):
>
> ```auto
> solver.report = True
> solver.line_search = True # Enable line search
> 
> ```

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)` ?

---

<div class="post-metadata">

**Author:** ![Jo\_A](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jo_a/32/8444_2.png) [@Jo\_A](https://fenicsproject.discourse.group/u/Jo_A)\
**Post date:** [May 7, 2025, 11:05pm UTC](https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559/3 "2025-05-07T23:05:13Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![Jo\_A](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jo_a/32/8444_2.png) [@Jo\_A](https://fenicsproject.discourse.group/u/Jo_A)\
**Post date:** [May 11, 2025, 9:29pm UTC](https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559/4 "2025-05-11T21:29:00Z")

</div>

@dokken  
In fact the problem comes from the use of the Discontinuous Galerkin in order to impose a constraint through a Lagrange multiplier which is not the right way.  
It seems that I must use **dolfinx\_mpc**.  
However, it is not clear how to implement it for non linear problems.  
There is [this conversation](https://fenicsproject.discourse.group/t/dolfinx-mpc-solve-nonlinear-problems/9854) and [this one](https://fenicsproject.discourse.group/t/unequal-results-for-linear-and-nonlinear-solver-with-dolfinx-mpc/10208), but the file [_test\_nonlinear\_assembly.py_](https://github.com/jorgensd/dolfinx_mpc/blob/main/python/tests/test_nonlinear_assembly.py) seems not to correspond to the latest versions of dolfinx… ?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 12, 2025, 1:38am UTC](https://fenicsproject.discourse.group/t/newton-solver-convergence-of-a-non-linear-problem-for-a-mixed-formulation-lagrange-multiplier-in-1d/17559/5 "2025-05-12T01:38:27Z")

</div>

You mean 0.9.x? See: [dolfinx\_mpc/python/tests/test\_nonlinear\_assembly.py at v0.9.1 · jorgensd/dolfinx\_mpc · GitHub](https://github.com/jorgensd/dolfinx_mpc/blob/v0.9.1/python/tests/test_nonlinear_assembly.py)