I am attempting to solve the electrolyte domain of li-ion battery model. The concentration problem is somewhat like a nonlinear diffusion/heat equation but it gets messy because it has a div(grad(ln(c))) term. So I would be looking at something like this:

-div(k1*grad(c)) + div(k2*grad(ln(c))) = 0

Which we can rewrite as:

-div((k1 - k2/c)*grad(c)) = 0

If k2 is relatively small enough, the newton solver should converge with little problem. The problem I am having is with the linear solve. If I discretize the above using CG1, no solver works so long as k2 >> 0. With DG1, I can only use direct solvers like LU to solve the resulting linear system of equations. I can't seem to find a good iterative solver/pc to solve this - HYPRE gives me a DIVERGED_NANORINF, and jacobi/bjacobi do not converge at all.

I've attached the code for you guys to play around with. I don't know what is a good PC for grad(ln(c)), my guess is that I'm going to need to provide an assembled preconditioner which approximates the inverse of the k2/c*grad(c) term but it doesn't seem clear to me how I can provide an explicit PC to the nonlinear solver without resorting to hacking at petsc4py.

Any thoughts appreciated

Justin