Initial guess for newton in coupled non-linear pde system

Question

Hello,

I want to solve a non-linear coupled PDE with a dispersion term (Klausmeier-Model for vegetation patterns). The newton solver does not converge and I want to provide an intial guess for the two fields "w" and "n" in each step (just as in http://fenicsproject.org/qa/3657/providing-initial-guess-for-newton_solver).

My questions:

1. How to linearize coupled PDEs if "du" is the trailfunction and not "u"? "du" only appears as in a = derivative(L, u, du). It is already linear in du! Should I instead linearize L?

2. Is there another way to provide a good initial guess (e.g. short grid-searching?)

For completeness, I attached the relevant pieces of code below. It is essentially an adapated version of Cahn-Hillard.demo

Class for interfacing with the Newton solver

class Klausmeierquation(NonlinearProblem):
def __init__(self, a, L):
    NonlinearProblem.__init__(self)
    self.L = L
    self.a = a
    self.reset_sparsity = True
def F(self, b, x):
    assemble(self.L, tensor=b)
def J(self, A, x):
    assemble(self.a, tensor=A, reset_sparsity=self.reset_sparsity)
    self.reset_sparsity = False

du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

Define functions

u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

Split mixed functions (water and biomass)

dw, dn = split(du)
w,  n  = split(u)
w0, n0 = split(u0)

Create intial conditions and interpolate

u_init = InitialConditions()
u.interpolate(u_init)
u0.interpolate(u_init)

Create derivative

udx = project(u.dx(0), ME) # u.dx(1) for y-derivative
wdx,ndx = split(udx)

mu_(n+theta), theta = 0.5

w_mid = (1.0-theta)*w0 + theta*w
n_mid = (1.0-theta)*n0 + theta*n

Weak statement of the equations (Parameters: A = rain, M = mortality, S = water flow)

L0 = w*v*dx - w0*v*dx - A*v*dt*dx + w0*v*dt*dx + \\
w0*n_mid*n_mid*v*dt*dx - S*wdx*v*dt*dx

L1 = n*q*dt*dx - n0*q*dt*dx - w_mid*n0*n0*q*dt*dx \\
+ M*n0*q*dt*dx + dot(grad(n0), grad(q))*dt*dx

L = L0 + L1

Compute directional derivative about u in the direction of du (Jacobian)

a = derivative(L, u, du)

Create nonlinear problem and Newton solver

problem = Klausmeierquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

Time-Step

while (t < T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    file << (u.split()[1], t)