Kernel appears to be dead when tried to execute solve() command

General
1

Hi, I am having difficulty in executing the following block of my code -

# Compute solution
u = Function(V, name = "Displacement")
solve(a == L, u, bcs, solver_parameters={'linear_solver':'mumps'})

MWE is -

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

tol = 1e-5

#Making a geometry
cylinder = Cylinder(Point(0.5, 0.5, 0), Point(0.5, 0.5, 2), 0.25, 0.25)
cube = Box(Point(0, 0, 0),Point(1, 1, 1))

domain = cube - cylinder

# Making Mesh (40 corresponds to the mesh density)
mesh = generate_mesh(domain, 40)

V = VectorFunctionSpace(mesh, 'CG', 2)
tol = 1E-8

def bottom(x, on_boundary):
    return near(x[0], 0., tol)

def top(x, on_boundary):
    return near(x[0], 1., tol)

bcbot = DirichletBC(V,  Constant((0, 0, 0)), bottom)
bctop = DirichletBC(V.sub(0),  Constant(0.1), top)
bcs = [bcbot, bctop]

E = Constant(200e3) #E = 200 GPa 
nu = Constant(0.3)
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

# Define strain and stress

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return lmbda*tr(epsilon(u))*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V, name = "Displacement")
solve(a == L, u, bcs, solver_parameters={'linear_solver':'mumps'})

It’s is working fine for a unit cube mesh but when I try to change the domain like here a cube with a cylindrical cavity in it, it takes too long before showing kernel died.
I hope to know what mistake I am making and looking forward to solve it.
Thank You

2

Have you tried on a coarser (less dense mesh) generated by mshr?

I would guess you run out of memory.

3

Thanks @dokken for this information, it works fine till mesh density = 35 and beyond that the kernel dies, but what alternative should I use if I don’t want to make my mesh coarser and execute as it is ?

4

Use an iterative solver.
As you are solving what appears to be an elasticity problem, using CG with a GAMG preconditioner is usually efficient.

You could alternatively use the setup from Bitbucket

5

Hi @dokken, thanks for the suggestion, I tried to use the code from Bitbucket after some modifications but ran into error while executing this cell -

# Assemble system, applying boundary conditions and preserving
# symmetry)
A, b = assemble_system(a, L, bcs)

Error -

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
/tmp/ipykernel_739/3299545883.py in <module>
      1 # Assemble system, applying boundary conditions and preserving
      2 # symmetry)
----> 3 A, b = assemble_system(a, L, bcs)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py in assemble_system(A_form, b_form, bcs, x0, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, A_tensor, b_tensor, backend)
    428         assembler.assemble(A_tensor, b_tensor, x0)
    429     else:
--> 430         assembler.assemble(A_tensor, b_tensor)
    431 
    432     return A_tensor, b_tensor

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to match boundary condition to function space.
*** Reason:  Function space on axis 0 does not contain BC space.
*** Where:   This error was encountered inside SystemAssembler.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

Minimum working example which I’m using -

from dolfin import *
import matplotlib.pyplot as plt


# Test for PETSc
if not has_linear_algebra_backend("PETSc"):
    print("DOLFIN has not been configured with PETSc. Exiting.")
    exit()

# Set backend to PETSC
parameters["linear_algebra_backend"] = "PETSc"

def build_nullspace(V, x):
    """Function to build null space for 3D elasticity"""

    # Create list of vectors for null space
    nullspace_basis = [x.copy() for i in range(6)]

    # Build translational null space basis
    V.sub(0).dofmap().set(nullspace_basis[0], 1.0);
    V.sub(1).dofmap().set(nullspace_basis[1], 1.0);
    V.sub(2).dofmap().set(nullspace_basis[2], 1.0);

    for x in nullspace_basis:
        x.apply("insert")

    # Create vector space basis and orthogonalize
    basis = VectorSpaceBasis(nullspace_basis)
    basis.orthonormalize()

    return basis

from dolfin import *
from mshr import *
import matplotlib.pyplot as plt

#Making a cylindrical geometry
cylinder = Cylinder(Point(0.5, 0.5, 0), Point(0.5, 0.5, 2), 0.25, 0.25)
cube = Box(Point(0, 0, 0),Point(1, 1, 1))

domain = cube - cylinder

# Making Mesh (40 corresponds to the mesh density)
mesh = generate_mesh(domain, 40)

# Stress and strain computation
E = Constant(200e3) #E = 200 GPa 
nu = Constant(0.3)
model = "plane_stress"

mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

# Define strain and stress

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return lmbda*tr(epsilon(u))*Identity(3) + 2*mu*epsilon(u)

# Create function space
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = inner(sigma(u), epsilon(v))*dx
L = inner(f, v)*dx

# Set up boundary condition 
class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0., 0.1)
class right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1., 0.1)

boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)

left().mark(boundaries, 1)
right().mark(boundaries, 2)

V = VectorFunctionSpace(mesh, "Lagrange", 1)
bc_left = DirichletBC(V, Constant((0., 0., 0.)), boundaries, 1)
bc_right = DirichletBC(V.sub(0), Constant(.1) , boundaries, 2)

bcs = [bc_left, bc_right]

# Assemble system, applying boundary conditions and preserving
# symmetry)
A, b = assemble_system(a, L, bcs)