Newton solver does not converge for nonlinear problem with neumann boundary conditions

It is hard for us to give specific answers if you do not provide any more explicit information.

However, in general, it can very often be that a differential equation does not have a unique solution when only Neumann conditions are subscribed. Check whether you get a sensible answer if you only consider the third equation with the Neumann conditions only.

If that does not clarify, try posting the _smallest possible_ (but running) piece of code that reproduces your problem.

Thanks Marie!

I tried considering only the third equation, however the built-in Newton solver still does not converge.

Here is the 'smallest possible' running code (i know, it's still kind of long;-):

# START OF CODE
import numpy
from dolfin import *

# Create mesh and define function space
meshsize = 32
timesteps = 10
mesh = UnitSquare(meshsize,meshsize)
V = FunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V,V,V])
Q = VectorFunctionSpace(mesh, "CG", 1)

# Define Constants
dt = 10**(-4) # Time step
eps = 10**(-3)
lamda = 10**(-2)
Exstart = 0.0
Eystart = 0.0

################################################################################
# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS \
        or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

################################################################################
class InitialConditions(Expression):
 def eval(self, values, x):
     if (x[0] < 1 - DOLFIN_EPS and x[0] > 0 + DOLFIN_EPS) and (x[1] < 1 - DOLFIN_EPS and x[1] > 0 + DOLFIN_EPS):
              values[0] = 1/sqrt(2*pi)*exp(-1/2*((x[0]-0.5)/(1)*4)**2)*exp(-1/2*((x[1]-0.5)/(1)*4)**2) + 1.0 + 0.01
              values[1] = 0.0
              # values[2] = sqrt(sin(2*pi*x[0]))*sqrt(sin(2*pi*x[1]))
              values[2] = 0.0
         elif x[0] > 1 - DOLFIN_EPS:
              values[0] = 1.0
              values[1] = 0.0
              values[2] = 0.0
         elif x[1] > 1 - DOLFIN_EPS:
              values[0] = 1.0
              values[1] = 0.0
              values[2] = 0.0
         elif x[0] < 0 + DOLFIN_EPS:
              values[0] = 1.0
              values[1] = 0.0
              values[2] = 0.0
         elif x[1] < 0 + DOLFIN_EPS:
              values[0] = 1.0
              values[1] = 0.0
              values[2] = 0.0
     def value_shape(self):
         return (3,)

################################################################################
class State(object):
 def __init__(self):
  self.uk = Function(W)
  self.u0 = Function(W)
  self.E = Function(Q)
  self.problem = None

 def solve(self, u0, E):

  self.E = E
  self.u0 = u0
  self.uk = Function(W, self.u0)
  self.du = TrialFunction(W)
  self.F = -(1.0/dt)*(self.uk[0]**2 - self.u0[0]**2)*v[0]*dx - inner(self.uk[0]**2*(grad(self.uk[1]) + self.E), grad(v[0]))*dx \
    - eps*inner(grad(self.uk[0]), grad(v[1]))*dx - self.uk[0]*(self.uk[2]-self.uk[1])*v[1]*dx \
    + lamda*inner(grad(self.uk[2]),grad(v[2]))*dx - self.uk[0]**2*v[2]*dx \
  self.J = derivative(self.F, self.uk, self.du)
  self.du = Function(W)
  self.problem = VariationalProblem(self.F, self.J, bc)
  self.problem.solve(self.uk)

  return self.uk

################################################################################
class StateT(object):
 def __init__(self):
  self.state = State()
  self.u = [Function(W)]*(timesteps + 1)

 def solve(self, uinit, E):
  self.u[0].vector()[:] = uinit.vector()
  self.E = E
  for i in range(0,timesteps):
   self.u[i+1] = Function(self.state.solve(self.u[i], self.E[i+1]))

  return self.u

################################################################################
class Functional(object):
 def __init__(self, E):
  self.statet = StateT()
  self.E = E
  self.uinit = Function(W)
  self.uinit = InitialConditions()
  self.uinit = interpolate(self.uinit, W)
  self.u = [Function(W)]*(timesteps + 1)
  self.u = self.statet.solve(self.uinit, self.E)

################################################################################
E = [Function(Q)]*(timesteps + 1)
for i in range(0,timesteps+1):
 E[i] = Constant((Exstart, Eystart))
 E[i] = interpolate(E[i],Q)
# Define Test Function for Forward Problem:
v = TestFunction(W)
# Define boundary condition for Forward Problem (for the first two equations):
bvalue = Constant(1.0)
bc = [DirichletBC(W.sub(0), bvalue, boundary), DirichletBC(W.sub(1), bvalue, boundary)]

functional = Functional(E)
# END OF CODE

The problem appears in the class state where the nonlinear variational problem is suppossed to be solved.
I have also tried different timegrids and meshsizes, but it doesn't help.

Have i done a coding mistake for the Neumann condition (i mean i simply skipped it in F)?
Again, any help is highly appreciated!

You'll need to make your code simpler to increase the likelihood of someone looking at it closely.

I wanted to provide a running piece of code, and this is the smallest i can get it without changing pieces that might influence the solution (as far as i see it).

But the solver is called in the class state, so this is most likely the place where the problem appears, so i just post this class again:

#START OF CODE
class State(object):
 def __init__(self):
  self.uk = Function(W)
  self.u0 = Function(W)
  self.E = Function(Q)
  self.problem = None

 def solve(self, u0, E):
  self.E = E
  self.u0 = u0
  self.uk = Function(W, self.u0)
  self.du = TrialFunction(W)
  self.F = -(1.0/dt)*(self.uk[0]**2 - self.u0[0]**2)*v[0]*dx - inner(self.uk[0]**2*(grad(self.uk[1]) + self.E), grad(v[0]))*dx \
    - eps*inner(grad(self.uk[0]), grad(v[1]))*dx - self.uk[0]*(self.uk[2]-self.uk[1])*v[1]*dx \
    + lamda*inner(grad(self.uk[2]),grad(v[2]))*dx - self.uk[0]**2*v[2]*dx \
  self.J = derivative(self.F, self.uk, self.du)
  self.du = Function(W)
  self.problem = VariationalProblem(self.F, self.J, bc)
  self.problem.solve(self.uk)

  return self.uk
#END OF CODE

Some tricks to simplify this code for debugging:

1) Get rid of all the self., they make everything less readable.

2) Replace problem coefficients with constants, 1 or 0 if possible to
simplify expressions.

3) Remove the class, make it a flat script.

Do this step by step in a fashion that lets you test as small changes as
possible.

Martin
Den 8. juni 2011 09.21 skrev "Sebastian Rau" <
<email address hidden>> følgende:
> Question #160453 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/160453
>
> Sebastian Rau gave more information on the question:
> I wanted to provide a running piece of code, and this is the smallest i
> can get it without changing pieces that might influence the solution (as
> far as i see it).
>
> But the solver is called in the class state, so this is most likely the
> place where the problem appears, so i just post this class again:
>
> #START OF CODE
> class State(object):
> def __init__(self):
> self.uk = Function(W)
> self.u0 = Function(W)
> self.E = Function(Q)
> self.problem = None
>
> def solve(self, u0, E):
> self.E = E
> self.u0 = u0
> self.uk = Function(W, self.u0)
> self.du = TrialFunction(W)
> self.F = -(1.0/dt)*(self.uk[0]**2 - self.u0[0]**2)*v[0]*dx - inner(self.uk
[0]**2*(grad(self.uk[1]) + self.E), grad(v[0]))*dx \
> - eps*inner(grad(self.uk[0]), grad(v[1]))*dx - self.uk[0]*(self.uk[2]-
self.uk[1])*v[1]*dx \
> + lamda*inner(grad(self.uk[2]),grad(v[2]))*dx - self.uk[0]**2*v[2]*dx \
> self.J = derivative(self.F, self.uk, self.du)
> self.du = Function(W)
> self.problem = VariationalProblem(self.F, self.J, bc)
> self.problem.solve(self.uk)
>
> return self.uk
> #END OF CODE
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Thanks for your answers, i will test this.

Can i also specify a normalization condition (something like integral over uk[2] - this is the variable with neumann boundaries - should be equal to 0) in order to ensure uniqueness of a solution? How would i code that?

I had a look at the demo file of the nonlinear Poisson Problem with Neumann boundaries, and added the normalization condition. It seems it works with that...