Hi All,
I want to solve the steady-state PNP (Poisson-Nernst-Planck) equation using FEniCS. The domain is nano-scale 2D square domain with 2 ion species: K+ and Cl-.
I first solve the Poission equation to get electric potential ( concentrations of K+ and Cl- are fixed ), then put this electric potential into Nernst-Planck equation and solve for new concentrations of K+ and Cl-, respectively. Repeat above cycle iteratively until convergence is reached. Following is the pseudo code I used:
def solvePBproblem(cationconc,anionconc):
......... receive concentrations of K+ and Cl-,
solve for electric potential phi...........
return phi
def solveNPproblem(electric_potential,charge):
......... receive electric potential, charge,
solve for new concentration ...........
return new_concentration
while True: # do iterations until convergent
phi = solvePBproblem(cp1,cp2)
cp1 = solveNPproblem(phi,z1)
cp2 = solveNPproblem(phi,z2)
In my calculations, the convergence was reached only when I set the element size to be extremely small. If I slightly increased the element size, the calculation would abort after several iterations ( the electric potential suddenly became unexpectedly large ). I also tried to change the configurations of the solver from default to
solver.parameters["linear_solver"] = "minres" # also tried umfpack and gmres
solver.parameters["preconditioner"] = "amg" # also tried ilu as well
All above tests failed to reach convergence with the slightly increased element size. My question is : Is using the extremely small element size the only way to reach convergence or could I use some more advanced approaches with a fairly small element size to achieve convergence ?
Any suggestions would be highly appreciated !