HyperElastic Example with second order vectorfunctionspace failed

Question

Hi all,

I am running the hyperelastic demo. with first order functionspace, it was running perfectly,
but when I changed to second order, it gave error on out of memory, which does not make sense to me.

I am running on a machine with i7-4770K and 16 GB memory, mint 17.1. A double check using ***top*** in the terminal did not reveal any memory issue.

Any thoughts will be appreciated. Thanks.

Here's error infor:

 Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 6.943e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)

UMFPACK V5.6.2 (Apr 25, 2013): ERROR: out of memory

Traceback (most recent call last):
  File "/home/haibo/Research_Projects/Projects/Mech_Atria/Test_Poisson/Test_Mesh_degree.py", line 67, in <module>
    form_compiler_parameters=ffc_options)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 297, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 345, in _solve_varproblem
    solver.solve()
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@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76.
*** Where:   This error was encountered inside /build/dolfin-k_QrtL/dolfin-1.6.0/dolfin/la/PETScLUSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.6.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

[Finished in 895.5s with exit code 1]
[shell_cmd: python -u "/home/haibo/Research_Projects/Projects/Mech_Atria/Test_Poisson/Test_Mesh_degree.py"]
[dir: /home/haibo/Research_Projects/Projects/Mech_Atria/Test_Poisson]
[path: /usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games]

And here's the source code:

from dolfin import *

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
                # "quadrature_degree":2,\
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True
               }

# Create mesh and define function space
# mesh = UnitCubeMesh(24, 16, 16)
mesh = UnitCubeMesh(24,16,16)
V = VectorFunctionSpace(mesh, "Lagrange", 2)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"))
r = Expression(("scale*0.0",
                "scale*(y0 + (x[1] - y0)*cos(theta) - (x[2] - z0)*sin(theta) - x[1])",
                "scale*(z0 + (x[1] - y0)*sin(theta) + (x[2] - z0)*cos(theta) - x[2])"),
                scale = 0.5, y0 = 0.5, z0 = 0.5, theta = pi/3)

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, -0.5, 0.0))  # Body force per unit volume
T  = Constant((0.1,  0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)

# Save solution in VTK format
file = File("displacement_odr2.pvd");
file << u;

# Plot and hold solution
plot(u, mode = "displacement", interactive = True)

hernan_mella · Answer 1 · Nov 23, 2015

The increase in the memory usage when you use a direct solver and P2 elements is very high.
May be you can use a combination of these alternatives:

Use an iterative method (e.g. gmres) to solve the linear system in each iteration of the newton algorithm.
Reduce the number of quadrature points.
Run your code in parallel (use mpirun -n np python your_code.py, where np is the number of processes).

The above suggestions result in a modification of your code as follows:

# form compiler parameters
parameters["form_compiler"]["representation"] = "quadrature"
parameters["form_compiler"]["quadrature_degree"] = 5

# Defining non linear problem
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
#prm['newton_solver']['absolute_tolerance'] = 1E-8
#prm['newton_solver']['relative_tolerance'] = 1E-7
#prm['newton_solver']['maximum_iterations'] = 25
#prm['newton_solver']['relaxation_parameter'] = 1.0
prm['newton_solver']['linear_solver'] = 'gmres'

# Solve variational problem
solver.solve()

HyperElastic Example with second order vectorfunctionspace failed

1 Answer