Problems solving a NoLinearProblem

General
1

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

    # Read the mesh

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

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

    # 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:
                https://link.springer.com/content/pdf/10.1023/A:1022235512626.pdf
                """
                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]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    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)
    nh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    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):
    return sym(nabla_grad(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:

Thank you in advance!

2

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

3

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 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.

I obtain the same error with the custom Newton solver