Iterative solver for incompressible solid

Linear Algebra
1

Hello,

I am trying to obtain the solution for a unit cube with an incompressible hyperelastic density energy function. The cube is loaded only with displacement boundary conditions in two opposed faces. If I use a direct solver as “mumps”, results are obtained successfully. However, if I use an iterative solver as “cg” with any kind of preconditioner, code fails, either after the first iteration or when the maximum number of iterations is reached.

Can somebody help me? I don’t know the reason for this kind of failure just at the beginning of the solution and even when I impose very small displacements. Thank you

# Clear all variables
import dolfin
import time

t1 = time.time()

# coarse deformation gradient (target lambdas)
lam1t = 1.01
lam2t = 1.0
lam3t = 1.0 

# Elasticity parameters
E, nu = 5.0, 0.4999

# paths
results_path  = './'

dolfin.set_log_level(20)
dolfin.parameters['form_compiler']['cpp_optimize'] = True
dolfin.parameters['form_compiler']['representation'] = "uflacs"
dolfin.parameters["form_compiler"]["quadrature_degree"] = 2


# ABOUT MESH

# create mesh
nelem_side = 15
mesh = dolfin.UnitCubeMesh.create(nelem_side, nelem_side, nelem_side, dolfin.CellType.Type.tetrahedron)
outer_boundary_markers = dolfin.MeshFunction('size_t',mesh, mesh.topology().dim() - 1)
inner_boundary_markers = dolfin.MeshFunction('size_t',mesh, mesh.topology().dim() - 1)

tol = 0.01
class BoundaryX0(dolfin.SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary and dolfin.near(x[0],0.0,tol)

class BoundaryX1(dolfin.SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary and dolfin.near(x[0],1.0,tol)

# initial setting
outer_boundary_markers.set_all(0)

# create instances of boundaries
bX0 = BoundaryX0()
bX1 = BoundaryX1()

# mark boundaries
bX0.mark(outer_boundary_markers, 1)
bX1.mark(outer_boundary_markers, 4)

print('Mesh loaded')



# ABOUT ELEMENTS


# modeling parameters
mu  = dolfin.Constant(E/(2.0*(1.0 + nu)))


# DEFINE ELEMENT DOMAIN!
P2 = dolfin.VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = dolfin.FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = dolfin.MixedElement([P2, P1]) 
W  = dolfin.FunctionSpace(mesh, TH)



# ABOUT FUNCTION SPACES

w = dolfin.Function(W)
v = dolfin.TestFunction(W)
du = dolfin.TrialFunction(W)

(u, p) = dolfin.split(w)
(v, delta_p) = dolfin.split(v)



# ABOUT BOUNDARY CONDITIONS


# imposed displacements on outer boundary
zero = dolfin.Constant((0.0,0.0,0.0))
disp_value = dolfin.Expression(('(lmbda1-1.0)*x[0]', '(lmbda2-1.0)*x[1]', '(lmbda3-1.0)*x[2]'), 
                               lmbda1 = lam1t, lmbda2 = lam2t, lmbda3 = lam3t, degree=1, domain=mesh)

bc_xmin = dolfin.DirichletBC(W.sub(0), zero, outer_boundary_markers, 1)
bc_xmax = dolfin.DirichletBC(W.sub(0), disp_value, outer_boundary_markers, 4)
DBC = [bc_xmin, bc_xmax]

    

# ABOUT FORMULATING THE PROBLEM

I = dolfin.Identity(3)             # Identity tensor
F = I + dolfin.grad(u)             # Deformation gradient
C = F.T*F                          # Right Cauchy-Green tensor
Ic = dolfin.tr(C)
J  = dolfin.det(F)

# Piola-Kirchhoff stress tensor for simply incompressible Neo-Hookean energy
# psi = mu/2*(I3-1) - p(J-1)
Finv_t = dolfin.inv(F).T
P = mu*F  - p*J*Finv_t

# nonlinear variational equation
F = dolfin.inner(P, dolfin.grad(v) )*dolfin.dx + delta_p*(J-1.0)*dolfin.dx
Jac = dolfin.derivative(F, w, du)


# ABOUT SOLUTION      

# parameters
newton_solver_dict = {'convergence_criterion': 'residual',
                      'relative_tolerance': 1E-12,
                      'absolute_tolerance':1e-10,
                      'maximum_iterations':16,
                      'report':True,
                      'relaxation_parameter': 1.0,
#                      'linear_solver':'mumps',
                      'linear_solver':'cg',
                      'preconditioner':'default'
                      }

#solution
dolfin.solve(F == 0, w, DBC, J=Jac, solver_parameters={"newton_solver": newton_solver_dict}) 

#projection solution
u = dolfin.split(w)[0]
U_proj = dolfin.VectorFunctionSpace(mesh, 'CG', 1)
u_proj = dolfin.project(u, U_proj, solver_type='mumps')

#RVE.solve(RVE.solver, nsteps, disp_value, lam1t, lam2t, lam3t)
file_ob_u = dolfin.File(results_path+'defmap.pvd')
file_ob_u << u_proj
2

Hello,
preconditioning for incompressible materials is an issue since the system becomes indefinite. You should check preconditioning solutions for Stokes flow for instance:
https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/stokes-iterative/python/documentation.html for a simple strategy
or https://fenicsproject.org/qa/12856/fieldsplit-petscpreconditioner_set_fieldsplit-arguments/ with PETSc fieldsplit

1 Like