nonlinear system with different user-defined source terms

Question

I'm trying to solve a modified nonlinear 1D Poisson equation. The source is nonlinear with f=exp(-u)-exp(u) (a Poisson-Boltzmann system). The modification is two different functions p0(x) and p1(x) overlaid on the nonlinear source terms, f=p0exp(-u)-p1exp(u). In general p0 is different from p1.

My sample dolfin script seems to run fine when p0=p1, the solution looks reasonable. But when I try to use two distinct perturbation functions (different magnitudes), dolfin does not converge. The error remains unchanged after iteration 1:

$ python test2.py 
Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 1.414e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 5: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 6: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 7: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 8: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 9: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
  Newton iteration 10: r (abs) = 4.672e-03 (tol = 1.000e-10) r (rel) = 3.303e-04 (tol = 1.000e-09)
Traceback (most recent call last):
  File "test2.py", line 60, in <module>
    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@fenicsproject.org
***
*** 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.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: unknown
*** 
*** DOLFIN version: 1.5.0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

Below is a minimal test code. The perturbation function is defined as an Expression subclass in UserDefinedExpression, instantiated in pList and used as p[0] and p[1] in source(u,p). It converges successfully as given, with secondAmplitude=firstAmplitude. But if I set secondAmplitude to a different value (e.g. with the line commented out), I get failure to converge.

I am new to dolfin, so it's quite possible I've done something stupid in the script.

===test.py===

from dolfin import *
import numpy

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

# This user-defined Expression adds a pertubation
# to the nonlinear source term
class UserDefinedExpression(Expression):
  def __init__(self, a):
    self.a = a

  def eval(self, values, x):
    values[0] = self.a*exp(-x[0])


# There are two source terms (exp(-u) and -exp(u)).
# I want to add different perturbation Expressions p[0] and p[1] to each of them
def source(u,p):
  return ( exp(-u)*exp(-p[0])
           -exp(u)*exp(-p[1]) )

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

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

a = inner(nabla_grad(u), nabla_grad(v))*dx

firstAmplitude=5
secondAmplitude=firstAmplitude    #  succeeds with this line when the two amplitudes are equal
#secondAmplitude=2*firstAmplitude   #  fails with this line to converge when the two amplitudes are different

pList=[]
pList.append(UserDefinedExpression(firstAmplitude))
pList.append(UserDefinedExpression(secondAmplitude))

f=source(u,pList)

L = f*v*dx

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

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

u=Function(V)

problem = NonlinearVariationalProblem(F, u, bc, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

plot(u, interactive=True)

Answer 1

Simple mistake towards the end. Try with these modifications

#u=Function(V)

# Use u_ for solution, because that is the Coefficient used by the forms F/J, 
# u_ is updated each iteration
problem = NonlinearVariationalProblem(F, u_, bc, J)
solver  = NonlinearVariationalSolver(problem)
solver.solve()

Comments

Perfect! Thanks Mikael. It solves my other question about linear self-consistency too.