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])
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
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))
# Deformation gradient
F = ufl.variable(I + grad(u))
print(grad(u)._dim)
# 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_pi = inner(ufl.grad(v_pi), ufl.grad(pi)) * ufl.dx - inner(v_pi, pi_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.
Thank you for your help!