Cannot solve nonlinear PDE

General
1

Hello,
I have been trying to solve the following, two-dimenional nonlinear PDE with Fenics:

2 grad(u).Hessian(u).grad(u) + (grad(u).grad(u)) * laplacian(u) == -f [in \Omega]
u = u_D [at the boundary of Omega]

where u is a function of two variables and Hessian(u) its Hessian. If I construct the problem by hand in such a way that u = 1 + x^2 + 2y^2 is a solution, then I have f = -8 * (5 * x^2 + 28 * y^2), u_D = u at the boundary of Omega. I gave this problem to fenics with the following code



from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = Function(V)
v = TestFunction(V)
f = Expression('-8*(5*x[0]*x[0]+28*x[1]*x[1])', degree = 2)

# Compute solution
dF = dot(grad(u), grad(u))*dot(grad(u), grad(v))*dx - f*v*dx
solve(dF == 0, u, bc)

# Save solution to file in VTK format
vtkfile = File('poisson/solution.pvd')
vtkfile << u

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')

# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
print('error_L2  =', error_L2)
print('error_max =', error_max)

But when I run it I get some nans

$ python3 ft01_poisson.py 

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 4.057e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 2: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 3: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 4: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 5: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 6: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 7: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 8: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 9: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 10: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 11: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 12: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 13: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 14: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 15: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 16: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 17: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 18: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 19: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 20: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 21: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 22: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 23: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 24: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 25: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 26: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 27: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 28: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 29: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 30: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 31: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 32: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 33: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 34: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 35: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 36: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 37: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 38: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 39: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 40: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 41: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 42: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 43: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 44: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 45: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  Newton iteration 46: r (abs) = nan (tol = 1.000e-10) r (rel) = nan (tol = 1.000e-09)
  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 "ft01_poisson.py", line 35, in <module>
    solve(dF == 0, u, bc)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 220, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 266, 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: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

What am I doing wrong ?

Thank you for your help

2

Besides the rather odd formulation, I suspect your initial guess will be catching you out. Can you try by interpolating (or projecting since this is old DOLFIN), your solution as the initial guess?

See also Default absolute tolerance and relative tolerance - #4 by nate.

3

I copied the code from an online sample program of nonlinear poisson equation, I don’t think that the formulation is odd.
How may I set the initial guess ?

4

See, e.g. Setting initial guess for nonlinear problem