# Poisson-Nernst-Planck equations with ill-conditioned sparse matrix

I am trying to solve Poisson-Nernst-Planck system of equations for ions diffusion problem using finite volume method. Nernst-Planck equation for mass transport and Poisson equation for electrostatic charge distribution.

My jacobian is sparse and ill-conditioned (condition number ~ 1E+14). For a system of 4 ions and three control volumes an schematic Jacobian looks like this:

I use Newton-Raphson method to solve it. The timesteps don't go larger than 1E-7 seconds and eventually it doesn't converge and crashes. In the mean time when I pause the run and check the results they are not correct.
One solution that I found was to assume $\epsilon/F = 0$ as $\epsilon/F \sim 1E-14$). In this case condition number becomes ~1E+12 although after ~ 2000 timesteps condition number jumps to 1E+30. But the code goes to larger timesteps (up to 1E-3 second) and the result seems right but not very precise. I need to speed up the code as this timestep is not practical.

1. Do you think using direct solver is an option?
2. Do you think higher quality initial guess for electric potential will help?
3. Is using solvers with quad precision helpful?
4. Is it possible that jacobian assembly is not right?

many thanks

• You might be able to diagnose your problem by simplifying it. For example, consider only ion diffusion, and drop all the electrostatic interactions. Then try using both diffusion and advection with a prescribed electric field, but no feedback of the ion concentration on the the electric field. Finally, add in the non-linearity. The PDE you're looking at is somewhat complicated, so any number of things could be going wrong. – Daniel Shapero Mar 28 '14 at 20:24
• Welcome to SciComp! This is a good question, but will get a better response if a bit more time is put into formatting. Consider using the built-in latex-like equation support, and perhaps naming some of the variables in the equations, as this will help people better understand the nature of your question. – meawoppl Mar 30 '14 at 1:56

For the past few months, I've been working on a similar problem. I only solve the steady-state potential equation for the moment, assuming homogeneous, steady-state ion concentrations.

Scaling your equations is going to be a rudimentary preconditioner. I found that I got good results using algebraic multigrid (AMG) for preconditioning instead, since my problem reduces to a coupled reaction-diffusion equation (i.e., a Poisson problem with nonlinear source term). You might try incomplete LU (ILU) to start, because it is a robust, slower preconditioner, and if that works, then try other preconditioners like AMG.

Using a direct solver might be helpful, but it really depends on the size of your problem. Even sparse LU decompositions will exhaust available memory due to fill-in for problems that are sufficiently large (that point occurs somewhere in the neighborhood of 100,000-1,000,000 variables). If your Jacobian matrix is rank-deficient (as was the case in my problem, due to a nonconductive filler material in some cells), then direct methods will also have problems.

A single higher quality guess for the nonlinear solver won't help. Your problem is transient, so you'll need a guess for every time step, if you're using an implicit method (and given that it's a diffusion problem, you probably should be, unless you want your time steps to be severely limited by the diffusion stability limit). Established libraries that implement numerical integration methods (i.e., time-steppers) do a good job of constructing initial guesses for each time step, and if the guesses are not great (they're typically derived from a problem-independent heuristic), the time steps will be reduced to preserve accuracy.

Quad precision might be helpful, if your problem defies preconditioning. There aren't many libraries I'm aware of that use quad precision. PETSc does, but I believe if you configure with quad precision, you're limited to PETSc-native solvers and preconditioners.

It's possible. You should test your assembled Jacobian against a finite difference approximation for debugging purposes. If you are using a finite difference approximation in your actual simulation, you will roughly halve the number of significant figures in your problem, which could severely erode your accuracy because your problem is ill-conditioned; that said, it doesn't seem like you're using that approximation in your simulation.

If you write your system as a single matrix, you may need to scale the two equations so that they have roughly the same order of magnitude. Rather than describing in great detail what to do and why, let me simply link to an example where we have done this: http://dealii.org/developer/doxygen/deal.II/step_32.html#Thescalingofdiscretizedequations

• If I understood it correctly, this method is not exactly normalization but more re-scaling. First equation dimension is [mol * m-3 * s-1], mass flux and second is [V * m-2], when it multiplies to F/epsilon. I multiplied the second equation to C * A * t-1 * V-1 [mol * m-1 * s-1 * V-1] = 1E-6. I got the similar results i.e small timestepping of 1E-7s. – Ben Feb 27 '14 at 2:18