Incompressible Linear Elasticity fails in the second time step

1

Hi,

I have an incompressible linear elasticity problem that I need to solve for displacement and pressure. I have the following code which converges once and then gives me a strange error. The code and the error are as follows:

from dolfin import *
import numpy as np

#Solver configs
class Problem(NonlinearProblem):
    def __init__(self, J, F, bcs):
        self.bilinear_form = J
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        assemble(self.linear_form, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self.bilinear_form, tensor=A)
        for bc in self.bcs:
            bc.apply(A)
###############################################################
###############################################################
class CustomSolver1(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)
        PETScOptions.set("ksp_max_it", 1000)
        PETScOptions.set("ksp_initial_guess_nonzero", "true")
        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("ksp_monitor")
        PETScOptions.set("pc_type", "hypre")
        PETScOptions.set("pc_hypre_type", "euclid")

        self.linear_solver().set_from_options()
###############################################################

#Mesh
meshlevel = 5
degree = 1
dim = 2
mesh = UnitDiscMesh.create(MPI.comm_world, meshlevel, degree, dim)
###############################################################

###############################################################
P1_scalar = FiniteElement("P",triangle,1)
P2_scalar = VectorElement("P",triangle,1)
VE = MixedElement([P2_scalar,P1_scalar])
V = FunctionSpace(mesh,VE)
###############################################################

# Define trial and test functions
up = Function(V)
u , p = split(up)
vu , vp = TestFunctions(V)
#######################################################################
#Dimension of the space
d = V.dim()
###############################################################


#time step variables
T = 2000      # final
num_steps= 20000 #number of time steps
dt = T / num_steps # time step
eps = 1       # diffusion coefficient
t=0
k = Constant(dt)
rho = Constant(eps)
###############################################################


#Stress function
def sigma(u,p):
     mu = Constant(5) #tumor is stiffer than regular tissue
     Lambda = Constant(13.867)
     II = Identity(2)            # Identity tensor
     return -p*II+2*mu*(epsilon(u))
###############################################################

#Symmetric part
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
###############################################################

vtkfile_1 = File('test/mesh.pvd')
#######################################################################
t = 0.0
j = int(0)
#######################################################################
vtkfile_1<<(mesh,t)
for n in range(num_steps):

     ##############################################################
     print(t,flush=True)
     t+=dt
     j+=1
     ##############################################################

     F2 =  inner(sigma(u,p),epsilon(vu))*dx+dot((div(u)-(1+t)),vp)*dx
     #######################################################################

     #######################################################################
     J2 = derivative(F2,up)
     problem1 = Problem(J2,F2,[])
     custom_solver1 = CustomSolver1()
     custom_solver1.solve(problem1, up.vector())

     ##############################################################
     u, p = up.split()
     ALE.move(mesh,u)

     vtkfile_1<<(mesh,t)

============= error stack trace ====================
[0] ERROR: zero diagonal in local row 1
iluk_seq_block file= ilu_seq.c line= 428

[0] called from: factor_private file= Euclid_dh.c line= 481
[0] called from: Euclid_dhSetup file= Euclid_dh.c line= 245
[0] called from: HYPRE_EuclidSetup file= HYPRE_parcsr_Euclid.c line= 290

application called MPI_Abort(comm=0x84000007, -1) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=-1
:
system msg for write_line failure : Bad file descriptor
Error in system call pthread_mutex_destroy: Device or resource busy
src/mpi/init/initthread.c:240

The error was really long so I posted the first part of it. IF you run the code you will see that it successfully converges and moves the mesh and then in the next time-step it fails with the above error. Do you know what is happening here and what I should change in my code to get it to work? Any help would be much appreciated.

2

Are you convinced your problem is well posed without boundary conditions?

Also you’re creating PETSc linear algebra objects every step in your loading loop. Euclid is probably complaining because the reference to the preconditioner has changed.

You could debug a lot of this by making sure your problem works with a coarse mesh and a direct solver. If that doesn’t work then you likely have more fundamental problems than the iterative solver.