# Divergence Issue with Nonlinear Solver and Mixed elements in FEniCS

Hello FEniCSx community,

I am encountering divergence issues with a nonlinear solver in FEniCS while attempting to solve a 2D problem with mixed elements. My goal is to solve a 2D mechanical flow field problem. The problem involves solving a system of nonlinear partial differential equations governing the deformation and stress distribution in a domain with specified boundary conditions. I have implemented the problem formulation using mixed elements, with deformation u represented as a vector and pressure pi as a scalar.

However, when attempting to solve the problem using a Newton solver in FEniCS, I encounter divergence issues. The residuals grow rapidly after the second iteration, r (abs) is already 3.00366e+32, and after the tenth iteration r (abs) is inf, causing the solver to diverge.

The field equations are as following written:

and the weak statement of the equations are as following written:

and my code is as following:

``````import numpy as np
import ufl
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx import nls, io
from dolfinx.mesh import *
from dolfinx.mesh import CellType, create_interval, locate_entities
from dolfinx.fem import (Constant, dirichletbc, Function, functionspace, assemble_scalar)
from dolfinx import fem, mesh, plot
from mpi4py import MPI
from basix.ufl import element, mixed_element
from petsc4py import PETSc
from ufl import dot, dx, grad, inner

############################# Mesh Part ############################

domain = create_rectangle(MPI.COMM_WORLD, np.array([[0, 0], [0.75, 0.75]]),
[150, 150], cell_type=CellType.triangle)
# Mixed element
P2 = element("Lagrange", domain.basix_cell(), 1, shape = (2,))    # Mixed element for tensor
P1 = element("Lagrange", domain.basix_cell(), 1)   # Mixed element for scalar without shape = (1,)

TH = ufl.MixedElement([P2, P1])
ME = functionspace(domain, TH)

v_u, v_pi = ufl.TestFunctions(ME)
Variable = Function(ME)
u, pi = ufl.split(Variable)

############################# Boundary Conditions ############################

def left_boundary(x):
return np.isclose(x[0], 0.0)
def right_boundary(x):
return np.isclose(x[0], 0.75)

boundary_u = PETSc.ScalarType((0, 0))
boundary_pi = PETSc.ScalarType(1)

fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left_boundary)
right_facets = mesh.locate_entities_boundary(domain, fdim, right_boundary)

marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack([np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

left_dofs = fem.locate_dofs_topological(ME.sub(0).collapse()[0], facet_tag.dim, facet_tag.find(1))
l_bc_u = fem.dirichletbc(boundary_u, left_dofs, ME.sub(0).collapse()[0])

right_dofs = fem.locate_dofs_topological(ME.sub(1).collapse()[0], facet_tag.dim, facet_tag.find(2))
r_bc_pi = fem.dirichletbc(boundary_pi, right_dofs, ME.sub(1).collapse()[0])

bcs = [l_bc_u,  r_bc_pi]

############################# Parameters ############################

B = fem.Constant(domain, default_scalar_type((0, 0)))
T_neumann = fem.Constant(domain, default_scalar_type((0, 1)))
pi_neumann = fem.Constant(domain, default_scalar_type(0))
# Spatial dimension
d = 2

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Elasticity parameters
E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2

# Hyper-elasticity
P = ufl.diff(psi, F)
T = P - pi * I

############################# Weak Statements ############################

F_u = inner(grad(v_u), T) * ufl.dx - inner(v_u, B) * ufl.dx - inner(v_u, T_neumann) * ds(2)
F = F_u + F_pi

problem = NonlinearProblem(F, Variable, bcs)
solver = NewtonSolver(domain.comm, problem)

############################# Post Processing ############################

# Set Newton solver options
solver.atol = 1e-8
solver.convergence_criterion = "incremental"
solver.rtol = 5.0e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
log.set_log_level(log.LogLevel.INFO)
r = solver.solve(Variable)
``````

Then I got the error as:
RuntimeError: Newton solver did not converge because maximum number of iterations reached

I think the reason could be setting of Neumann Conditions but I’m not certainly sure. Please reply if I didn’t explain enough and any suggestions, advice, recommendations or alternative solver strategies would be greatly appreciated.