The geometric nonlinear problem

Question

I want to solve the Problem

Newton iteration 47: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 48: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 49: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton iteration 50: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Traceback (most recent call last):
  File "demo_leaf.py", line 83, in <module>
    solve(F == 0, u_total)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 300, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 349, 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-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 solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2016.1.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Comments

Quoting the error message: "include a minimal running example to reproduce the error."

Answer 1

Since you have not shared your problem, consider the following suggestions to be generic in nature.

• The most common reason behind this is wrongly chosen/ physically incompatible or unsound parameters. Check if your physical parameters are realistic for the problem you are trying to solve. Write the units as comments in your code, so you are not using just any random values.

• Check the time step. Try playing around with it. Reduce it to a level after which you are sure that it cannot be a reason for non-convergence.

• Check your boundary conditions. The problem may not converge because of incorrect BC that do not make any sense physically.

• Increase the nonlinear solver iteration parameter to a high value above 50

  nonlinear_solver.parameters["maximum_iterations"] = 500;

• Try changing the value of your forcing function (source). It may be too much for your model parameters. I don't know how to put this in words technically but hope you get what I mean.

One of these corrections or a combination of them should get your code working considering that you have programmed it correctly. Play around with the values a little bit.

Comments

from dolfin import *
import numpy as np
import scipy as sp

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

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
    # Left boundary is "target domain" G
    def inside(self, x, on_boundary):
        return bool(x[0] < DOLFIN_EPS and x[0] > -DOLFIN_EPS and on_boundary)
    # Map right boundary (H) to left boundary (G)
    def map(self, x, y):
        y[0] = x[0] - 20.0
        y[1] = x[1]

# computation domain
L = 5.0
W = 1.0
h = 0.02

# growth epsilon parameters
n = 10
beta = 0.065
epsilon_g = Expression("0.065*pow(x[1]/1, 10)")

# create mesh 
mesh = RectangleMesh(Point(0.0, 0.0), Point(L, W), 20, 4, "right/left")
print("Plotting a RectangleMesh")
plot(mesh, title="Rectangle")

# define function space
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)

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

# define functions
du_total = TrialFunction(W)
u_total = Function(W)
v = TestFunction(W) 
u, zeta = split(u_total)

# Kinematics
epsilon = sym(grad(u)) + grad(grad(zeta))/2
epsilon_gg = as_matrix([[epsilon_g, 0], [0, 0]])
epsilon_e = epsilon - epsilon_gg

# Elasticity parameters
E = 1          # the Young 's modulus
mu = 0.3        # the Poisson 's ratio

# the bending energy density
F_b = E*pow(h, 3)/12/(1 - pow(mu, 2))*(pow(tr(grad(grad(zeta))), 2) + 2*(1 - mu)*det(tr(grad(grad(zeta)))))*dx

# the stretching energy density
F_s = E*h/(1 - pow(mu, 2))*(pow(tr(epsilon_e), 2) + 2*(1 - mu)*det(epsilon_e))*dx

# total potential energy
Psi = F_s + F_b

# compute first variation of Pi
F = derivative(Psi, u_total, v)

# computer jacobian of F
# J = derivative(F, u_total, du_total)

# solver variational problem
solve(F == 0, u_total)

# save solution in VTK format
file = File("output/leaf.pvd");
file << zeta

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

this is my code, the problem cannot solved!