linear solution not consistent with nonlinear solution

Question

Suppose $u$ is the nonlinear weak solution for a nonlinear system
$a.u=f(u)$. Then $u$ should also be the solution for the linear system $a.u=g$, where $g$ is fixed using the nonlinear solution, $g=f(u)$.

I'm not getting this self-consistency in FEniCS. Taking $a$ as the standard Poisson operator and solving over the unit interval with Dirichlet boundary conditions ($u$=10 on the boundary), I get a sensible nonlinear solution (with $u>0$ everywhere) with $f=\exp(-u)-\exp(u)$.

But when I apply a linear solver after projecting $f$ using the nonlinear solution, the linear solution diverges wildly from the nonlinear solution, picking up a large negative value at the midpoint.

errornorm returns a residual of 415.7, where I'd want it to be 1e-6 or less.

Is this a bug in FEniCS?

minimal code:

from dolfin import *
import numpy

# Sub domain for Dirichlet boundary condition
class DirichletBoundary(SubDomain):
  def inside(self, x, on_boundary):
    return on_boundary

def source(u):
  return exp(-u)-exp(u)

# Create mesh and function space
mesh = UnitIntervalMesh(200)
V = FunctionSpace(mesh, "CG", 1)

# Define boundary condition
u0 = Constant(10)
bc = DirichletBC(V, u0, DirichletBoundary())


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

# standard Poisson operator
a = inner(nabla_grad(u), nabla_grad(v))*dx

# first get nonlinear solution
f_nonlinear=source(u)
L_nonlinear = f_nonlinear*v*dx

F=a-L_nonlinear
u_=Function(V)    
F=action(F,u_)
J=derivative(F,u_,u)

u_nonlinear=Function(V)
problem = NonlinearVariationalProblem(F, u_nonlinear, bc, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

# then fix source term using the nonlinear solution to get the linear solution
f_linear = project(source(u_nonlinear),V)
L_linear = f_linear*v*dx

u_linear=Function(V)
solve(a == L_linear, u_linear, bc)

print('error norm between nonlinear and linear solutions = ' + str(errornorm(u_nonlinear,u_linear)) )

plot(u_nonlinear, title='Nonlinear solution')
plot(u_linear, title='Linear solution', interactive=True)

Answer 1

Mikael Mortensen gave me the solution in the answer to a separate question, http://fenicsproject.org/qa/7732/nonlinear-system-with-different-user-defined-source-terms?show=7743#a7743

The nonlinear system needs to be solved for the same Function used to define it.

So where I have

u_=Function(V)    
F=action(F,u_)
J=derivative(F,u_,u)

defining the nonlinear components using u_, I also need to define the problem and solve it using the same u_ (and not using a new Function u_nonlinear):

problem = NonlinearVariationalProblem(F, u_, bc, J)

The way I had it, J was not being updated correctly.