Is PetSc unable to handle multiple basic nonlinear ufl expressions? - PETSc error code is: 77, Petsc has generated inconsistent data

Errors
1

Hello,

I am working on a (nonlinear) contact problem which uses the euclidean norm of the displacement vectors in its penalty term + several trigonometric functions.

However, PetSc somehow did not accept expressions in the weak formulation and fails with:

Traceback (most recent call last):
  File "/home/fenics/shared/mwe_rec_cropped.py", line 122, in <module>
    n, converged = solver.solve(u)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/nls/petsc.py", line 41, in solve
    n, converged = super().solve(u.vector)
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 77, Petsc has generated inconsistent data

I tried different versions of this problem in order to find out when these kind of errors occour. The less nonlinear basic ufl functions I am using the more likely it is for PetSc to compile.

MWE:

import gmsh
import numpy as np
from dolfinx import fem, mesh, log
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from petsc4py.PETSc import ScalarType, Options
from ufl import dx, inner, derivative, tr, Identity, sym, grad, TrialFunction, TestFunction, Measure, SpatialCoordinate, dot, sqrt, cos, atan_2

gmsh.initialize()
gmsh.model.add("rectangle")
lc=0.5e-1

# generate points for inner contour
gmsh.model.geo.addPoint(1, 1, 0, lc, 1)
gmsh.model.geo.addPoint(-1, 1, 0, lc, 2)
gmsh.model.geo.addPoint(-1, -1, 0, lc, 3)
gmsh.model.geo.addPoint(1, -1, 0, lc, 4)
# generate inner contour
gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
# generate points for outer contour
gmsh.model.geo.addPoint(2, 2, 0, lc, 5)
gmsh.model.geo.addPoint(-2, 2, 0, lc, 6)
gmsh.model.geo.addPoint(-2, -2, 0, lc, 7)
gmsh.model.geo.addPoint(2, -2, 0, lc, 8)
# generate outer contour
gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 5, 8)
# generate curveloops
gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 1)
gmsh.model.geo.addCurveLoop([5, 6, 7, 8], 2)
# generate surface
gmsh.model.geo.addPlaneSurface([2, 1], 1)
gmsh.model.geo.synchronize()
# generate physical groups for boundary conditions
inside, outside, surface = 1, 2, 3
gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4], inside, name="Inner Contour")
gmsh.model.addPhysicalGroup(1, [5, 6, 7, 8], outside, name="Outer Contour")
gmsh.model.addPhysicalGroup(2, [1], surface, name="Surface")
# finish gmsh
gmsh.model.mesh.generate(2)

# create a DOLFINx mesh 
dim = 2
gmsh_model_rank = 0
domain_com = MPI.COMM_SELF
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, domain_com, gmsh_model_rank, dim)

gmsh.finalize()

tdim = domain.topology.dim
fdim = tdim - 1

# define function space and model parameters
V = fem.VectorFunctionSpace(domain, ("CG",1))
x = SpatialCoordinate(domain)
d = 2.4
r = d/2

# returns negative scalar if deformed configuration is still within a a circle r = D/2 acc. definition
def gapp_x(u, x):
     R = sqrt( pow(x[0], 2) + pow(x[1], 2)) 
     angle = atan_2( x[1], x[0] )
     ures = sqrt( pow(u[0], 2) + pow(u[1], 2))
     ux = cos(angle) * ures
     return  R - r - ux

metadata = {"quadrature_degree": 4}
ds = Measure('ds', domain=domain, subdomain_data=facet_markers, metadata=metadata)

# test and trial functions
u = fem.Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)

# define outer contour as homogenous boundary condition
u_zero = np.array((0,)*domain.geometry.dim, dtype=ScalarType)
bc = fem.dirichletbc(u_zero, fem.locate_dofs_topological(V, fdim, facet_markers.find(outside)), V)

# constituive relations - linear elasticity
E = fem.Constant(domain, ScalarType(10.))
nu = fem.Constant(domain, ScalarType(0.3))
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

# maculay bracket
def maculay(x): 
        return (x+abs(x))/2

pen = fem.Constant(domain, ScalarType(1e3))
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[0], maculay(-gapp_x(u, x)))*ds(1) 
J = derivative(form, u, du)

problem = fem.petsc.NonlinearProblem(form, u, bcs=[bc], J=J)
from dolfinx import nls
solver = nls.petsc.NewtonSolver(domain.comm, problem)

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

ksp = solver.krylov_solver
opts = Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(u)
assert(converged)

How is it possible to use multiple nonlinear basic functions in my weak formulation without getting a PetSc error?

Many thanks in advance!

2

Start by changing your iterative solver to be “preonly” and “lu”, then you will see that the residual blows up.

This is likely due tot the fact that you have an initial guess for u that is 0.

Setting an initial guess for u and reducing the size of your penalty parameter:

u.x.array[:] = 1
pen = fem.Constant(domain, ScalarType(1e1))

yields a converged solution.

If you increase the penalty parameter, a naive newton method is not sufficient, you would need adaptive descent direction steps to avoid oscillation between two solutions.

1 Like