# Problems solving a NoLinearProblem

Hi everyone!

I want to solve a Imcompressible flow, with Navier-Stokes equations.
Its a cavity with a conveyor.

Initial condition at inflow and top

U=(1,0)

Outflow

np-\frac{\partial u}{\partial n} = 0

Cavity

U = (0,0)

Conveyor

U = 0.1*Tangent

The code is:

import dolfinx
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import math
import ufl
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TrialFunction, VectorElement,
as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, MixedElement)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector,
create_vector, create_matrix, set_bc)

from dolfinx.fem import (Constant, Function, FunctionSpace,
assemble_scalar, dirichletbc, form, locate_dofs_topological, set_bc)

##########################################################################################################3

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh_cinta.xdmf", "r") as xdmf:

########################################################################################################

# Definition of Spaces. (dim and polinomial type)

v_cg1 = VectorElement("CG", domain.ufl_cell(), 1)
s_cg1 = FiniteElement("CG", domain.ufl_cell(), 1)
V = dolfinx.fem.FunctionSpace(domain, v_cg1)
Q = dolfinx.fem.FunctionSpace(domain, s_cg1)
Z = MixedElement([v_cg1, s_cg1])
W = dolfinx.fem.FunctionSpace(domain, Z)

########################################################################################################

# Definition of subdomains

x_up  = -40
x_down = 60
y_top = 20
y_floor = 0
x_left = 0
x_right = 10
y_cavity = -5

n = FacetNormal(domain)

def upstream(x):
return np.isclose(x[0], x_up)

def downstream(x):
return np.isclose(x[0], x_down)

def top(x):
return np.isclose(x[1], y_top)

def cavity_y(x):
return np.logical_or(np.logical_and(np.logical_or(np.isclose(x[1], y_floor), np.isclose(x[1], y_cavity)), np.logical_or(np.less(x[0], x_left) , np.greater(x[0], x_right))),np.isclose(x[1], y_cavity))

def cavity_x(x):
return np.logical_and(np.logical_or(np.isclose(x[0], x_left), np.isclose(x[0], x_right)),  np.less(x[1], 0.1))

# Whole cavity

def cavity(x):
return np.logical_or(cavity_x(x) , cavity_y(x))

def cinta(x):
return np.logical_and( np.logical_and(np.greater(x[0], x_left) , np.less(x[0], x_right)) , np.logical_and(np.greater(x[1],y_cavity) , np.less(x[1],5)) )

#############################################################################################################

#Define normal and tangent like functions

def facet_normal_approximation(V, mt: dolfinx.mesh.meshtags, mt_id: int, tangent=False, jit_options: dict = {},
form_compiler_options: dict = {}):

timer = dolfinx.common.Timer("~MPC: Facet normal projection")
comm = V.mesh.comm
n = ufl.FacetNormal(V.mesh)
nh = Function(V)
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
ds = ufl.ds(domain=V.mesh, subdomain_data=mt, subdomain_id=mt_id)
if tangent:
if V.mesh.geometry.dim == 1:
raise ValueError("Tangent not defined for 1D problem")
elif V.mesh.geometry.dim == 2:
a = ufl.inner(u, v) * ds
L = ufl.inner(ufl.as_vector([-n[1], n[0]]), v) * ds
else:
def tangential_proj(u, n):
"""
See for instance:
"""
return (ufl.Identity(u.ufl_shape[0]) - ufl.outer(n, n)) * u
c = dolfinx.fem.Constant(V.mesh, [1, 1, 1])
a = ufl.inner(u, v) * ds
L = ufl.inner(tangential_proj(c, n), v) * ds
else:
a = (ufl.inner(u, v) * ds)
L = ufl.inner(n, v) * ds

# Find all dofs that are not boundary dofs
imap = V.dofmap.index_map
all_blocks = np.arange(imap.size_local, dtype=np.int32)
top_blocks = dolfinx.fem.locate_dofs_topological(V, V.mesh.topology.dim - 1, mt.find(mt_id))
deac_blocks = all_blocks[np.isin(all_blocks, top_blocks, invert=True)]

# Note there should be a better way to do this
# Create sparsity pattern only for constraint + bc
#    bilinear_form = dolfinx.fem.form(a, jit_options=jit_options,
#                              form_compiler_options=form_compiler_options)
bilinear_form = dolfinx.fem.form(a)
pattern = dolfinx.fem.create_sparsity_pattern(bilinear_form)
pattern.insert_diagonal(deac_blocks)
pattern.assemble()
u_0 = dolfinx.fem.Function(V)
u_0.vector.set(0)

bc_deac = dolfinx.fem.dirichletbc(u_0, deac_blocks)
A = dolfinx.cpp.la.petsc.create_matrix(comm, pattern)
A.zeroEntries()

# Assemble the matrix with all entries
form_coeffs = dolfinx.cpp.fem.pack_coefficients(bilinear_form)
form_consts = dolfinx.cpp.fem.pack_constants(bilinear_form)
dolfinx.cpp.fem.petsc.assemble_matrix(A, bilinear_form, form_consts, form_coeffs, [bc_deac])
if bilinear_form.function_spaces[0] is bilinear_form.function_spaces[1]:
A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)
A.assemblyEnd(PETSc.Mat.AssemblyType.FLUSH)
dolfinx.cpp.fem.petsc.insert_diagonal(A, bilinear_form.function_spaces[0], [bc_deac], 1.0)
A.assemble()
#linear_form = dolfinx.fem.form(L, jit_options=jit_options,
#                        form_compiler_options=form_compiler_options)
linear_form = dolfinx.fem.form(L)
b = dolfinx.fem.petsc.assemble_vector(linear_form)

dolfinx.fem.petsc.apply_lifting(b, [bilinear_form], [[bc_deac]])
dolfinx.fem.petsc.set_bc(b, [bc_deac])

# Solve Linear problem
solver = PETSc.KSP().create(MPI.COMM_WORLD)
solver.setType("cg")
solver.rtol = 1e-8
solver.setOperators(A)
solver.solve(b, nh.vector)
timer.stop()
return nh

############################################################################################################

# Boundary conditions

ucavity_condition = Function(V)
uup_condition = Function(V)
pd_condition = Function(Q)

def u_exact(x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 1
return values

uup_condition.interpolate(u_exact)

ucavity_condition.x.array[:] = 0
pd_condition.x.array[:] = 0

with dolfinx.io.XDMFFile(domain.comm, "Plot/u_initial.xdmf", "w") as file:
file.write_mesh(domain)
file.write_function(uup_condition)

fdim = domain.topology.dim -1

up_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, upstream)
down_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, downstream)
cavity_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, cavity)
top_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, top)
cinta_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, cinta)

up_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, up_facets)
t_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, top_facets)
cinta_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, cinta_facets)
cavity_f = dolfinx.fem.locate_dofs_topological(W.sub(0), fdim, cavity_facets)
d_f = dolfinx.fem.locate_dofs_topological(W.sub(1), fdim, down_facets)

boundaries = [(1,up_facets),(2,down_facets),(3,top_facets),(4,cavity_facets),(5,cinta_facets)]

facet_indices, facet_markers = [], []
for (marker, face ) in boundaries:
facet_indices.append(face)
facet_markers.append(np.full_like(face, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

########################################################################################################################

# Normal an tangent around the conveyoy

nh = facet_normal_approximation(V, facet_tag, 5, tangent=False)
th = facet_normal_approximation(V, facet_tag, 5, tangent=True)

# Velocity at cinta

module = 0.1
th.x.array[:] = module*th.x.array[:]

#############################################################################################################

# Define ds measure

ds = Measure("ds", domain=domain, subdomain_data=facet_tag)

#######################################################################################################################

# Bounday conditions

bc_u = dirichletbc(uup_condition, up_f)
bc_t = dirichletbc(uup_condition, t_f)
bc_cavity = dirichletbc(ucavity_condition, cavity_f)
bc_cinta = dirichletbc(th, cinta_f)
bc_p = dirichletbc(pd_condition, d_f)

bc = [bc_u, bc_t, bc_cavity, bc_cinta, bc_p]

#########################################################################################################################

# Navier-Stokes Equations

mu = Constant(domain, PETSc.ScalarType(0.001))

def epsilon(u):

# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))

w = Function(W)
u, p = ufl.split(w)

m = TestFunction(W)

v, q = ufl.split(m)

a1 = dot(dot(u,nabla_grad(u)),v)*dx + div(u)*q*dx + inner(sigma(u,p),epsilon(v))*dx + dot(p*n, v)*ds - dot(mu*nabla_grad(u)*n, v)*ds

problem = dolfinx.fem.petsc.NonlinearProblem(a1, w, bcs=bc)
solver = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
n, converged = solver.solve(w)
assert(converged)
print(f"Number of interations: {n:d}")



when I try to solve, the kernel dies instantaneously.

Does anyone know where this problem is or how I can locate it?

The mesh is in:

Try running the code as a python script instead of inside a jupyter notebook (as it will hopefully give you some more output).

You could also use a custom newton solver to pinpoint where it crashes, see: Custom Newton solvers — FEniCSx tutorial

I tried run like script and obtain:

corrupted double-linked list (not small)
[carlos-PROX15-AMD:10783] *** Process received signal ***
[carlos-PROX15-AMD:10783] Signal: Aborted (6)
[carlos-PROX15-AMD:10783] Signal code:  (-6)
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind
[0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac OS X to find memory corruption errors
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run