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