Efficient primal--dual splitting methods for a Poisson-constrained JKO scheme for Poisson-Nernst-Planck models

Abstract: The Poisson--Nernst--Planck (PNP) equations strongly couple ionic transport and electrostatic interactions through the Poisson equation, posing substantial numerical challenges under small permittivity and complex potential boundary conditions. Underlying these equations is a natural Wasserstein gradient-flow structure, in which the Poisson equation serves as a local realization of the nonlocal electrostatic interaction energy. Exploiting this structure, we formulate each time step as a constrained convex minimization problem where the ionic continuity equations and the Poisson equation are incorporated as linear constraints, allowing the concentrations, fluxes, and electrostatic potential to be updated simultaneously. The variational structure of the scheme intrinsically guarantees the dissipation of the original free energy, mass conservation, and nonnegativity of ionic concentrations under general electrostatic boundary conditions. Moreover, the framework is structurally modular: extending from classical to modified PNP models with steric interactions and concentration-gradient corrections requires only modifying the energy functional, while all structure-preserving properties are automatically retained. To efficiently solve the resulting large-scale constrained problems, we develop preconditioned and transformed primal--dual algorithms equipped with tailored fast dual solvers, namely DCT-based direct and Schur-complement iterative methods, that exploit the coupled block structure of the PDE constraints. Numerical experiments on classical and modified PNP systems demonstrate the accuracy and structure-preserving properties of the scheme, and show that the proposed algorithms converge reliably in strongly coupled small-permittivity regimes without significant growth in computational cost.