Formulation with grad of function squared

variational formulation
1

Hi all,

The problem i’m trying to solve is:

\nabla( \nabla u^2) = 0
u = 1 on \delta \Omega

This is my implementation:

from fenics import *

mesh = UnitIntervalMesh(10)
V = FunctionSpace(mesh, "CG", 1)
u = Function(V)
v = TestFunction(V)

dt = Constant(0.1)

theta = u**2  # scaling

F = dot(grad(theta), grad(v))*dx

bcs = DirichletBC(V, 1, "on_boundary")
solve(F==0, u, bcs)

The above code diverges:

No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.414e+00 (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 "mwe_scaling.py", line 18, in <module>
    solve(F==0, u, bcs)
  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
*** -------------------------------------------------------------------------

Any idea how to solve this problem?

Cheers
Remi

2

I think you can think of this as

\nabla u ^2 = 2 u \nabla u

which leads to the nonlinear Poisson problem

- \nabla \cdot \kappa(u) \nabla u = 0

where \kappa(u) = 2 u. This is only well posed when \kappa \neq 0 which means we need a nonzero initial guess accordingly. So try:

...
u.interpolate(Constant(1.0))
solve(F==0, u, bcs)
...
1 Like
3

Hi @nate thank you very much! It all makes perfect sense.

The following transient problem works:

from fenics import *

num_steps = 3

mesh = UnitIntervalMesh(10)

V = FunctionSpace(mesh, "CG", 1)

u = Function(V)

u_n = Function(V)

v = TestFunction(V)

dt = Constant(0.1)

theta = u**2

theta_n = u_n**2

F = (theta-theta_n)/dt*v*dx

F += dot(grad(theta), grad(v))*dx

bcs = DirichletBC(V, 1, "on_boundary")

u.assign(interpolate(Constant(DOLFIN_EPS), V))

for i in range(num_steps):

    solve(F==0, u, bcs)

    u_n.assign(u)

Although for some reason the first iteration converges slower:

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) = 1.414e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 2.399e+00 (tol = 1.000e-10) r (rel) = 1.696e+00 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 4.505e-01 (tol = 1.000e-10) r (rel) = 3.185e-01 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 2.775e-02 (tol = 1.000e-10) r (rel) = 1.962e-02 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 1.077e-04 (tol = 1.000e-10) r (rel) = 7.616e-05 (tol = 1.000e-09)
  Newton iteration 5: r (abs) = 1.151e-09 (tol = 1.000e-10) r (rel) = 8.136e-10 (tol = 1.000e-09)
  Newton solver finished in 5 iterations and 5 linear solver iterations.
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) = 1.688e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 2.632e-01 (tol = 1.000e-10) r (rel) = 1.560e-01 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 6.272e-03 (tol = 1.000e-10) r (rel) = 3.716e-03 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 3.062e-06 (tol = 1.000e-10) r (rel) = 1.814e-06 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 5.172e-13 (tol = 1.000e-10) r (rel) = 3.065e-13 (tol = 1.000e-09)
  Newton solver finished in 4 iterations and 4 linear solver iterations.
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) = 7.105e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 4.813e-02 (tol = 1.000e-10) r (rel) = 6.774e-02 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 1.709e-04 (tol = 1.000e-10) r (rel) = 2.405e-04 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 1.712e-09 (tol = 1.000e-10) r (rel) = 2.409e-09 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 5.526e-15 (tol = 1.000e-10) r (rel) = 7.777e-15 (tol = 1.000e-09)
  Newton solver finished in 4 iterations and 4 linear solver iterations.

Any idea?

4

After you solve the system for the first time, you have a much better initial guess than a constant value of 1 everywhere. In essence you’re probably deeper in an attractor region. See some hints and tips here.

1 Like