PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Wei Wu, Jin Zeng, Zhen Zhang, Chaozhen Wei

arXiv:2608.30693v1math.NAmath.OCphysics.comp-ph

TL;DR

PNP discretization must resolve tightly coupled ionic transport and electrostatics while preserving key physical structures under challenging boundary conditions and small permittivity. The paper formulates a Poisson-constrained JKO scheme and solves it with primal–dual methods and tailored dual solvers. Experiments report robust, efficient behavior for classical and modified PNP systems, especially in strongly coupled small-permittivity regimes.

  • Problem

    PNP discretization must simultaneously preserve nonnegativity, mass conservation, and original-energy dissipation under general electrostatic boundary conditions and strong coupling.

  • Method

    The paper formulates each time step as a Poisson-constrained convex JKO minimization and applies preconditioned and transformed primal–dual algorithms with tailored fast dual solvers.

  • Results

    The algorithms show superior robustness and efficiency in experiments with various boundary conditions and fixed charges, especially for small dielectric permittivity.

  • Takeaways & Limitations

    The constrained JKO framework retains structure-preserving properties for classical and modified PNP models while supporting strongly coupled simulations without significant computational-cost growth.

  • Takeaways & Limitations

    The proposed scheme is first-order accurate in time; higher-order variational schemes remain future work.

Abstract

from arXiv · show

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.

1. Introduction

PNP simulation must handle nonlinear ionic–electrostatic coupling while preserving nonnegativity, mass, and original free-energy dissipation. The paper proposes a Poisson-constrained JKO formulation with primal–dual solvers for efficient, structure-preserving computation.

  • PNP numerics must control ionic–electrostatic coupling while preserving concentration nonnegativity, species mass, and original free-energy dissipation.
  • Simultaneously preserving these structures remains difficult with general electrostatic boundary conditions, multispecies coupling, affordable nonlinear solves, and concentration-gradient corrections.
  • The paper develops a unified Poisson-constrained JKO scheme for classical and modified PNP models.
  • Each time step becomes a series of convex minimization problems with concentrations and potential as independent variables coupled by linear Poisson constraints.
  • PrePD and VPTPD algorithms with tailored dual solvers converge reliably in strongly coupled small-permittivity regimes without significant computational-cost growth.
  • The variational formulation makes original-energy dissipation and concentration nonnegativity intrinsic, while positivity is enforced through the admissible Wasserstein transport action.
  • The JKO framework incorporates additional physical effects through energy functionals or coupling constraints without compromising structure-preserving properties.

2. Poisson-Nernst-Planck equations and its extensions

The paper presents classical and modified PNP models as coupled Wasserstein gradient flows with electrostatic interactions represented through a Poisson equation. Modified models add steric and concentration-gradient effects while retaining the same constrained framework.

  • 2.1. Poisson-Nernst-Planck equations: The models admit a coupled Wasserstein gradient-flow interpretation for the corresponding free energy, with concentration-dependent mobility and chemical-potential driving force.
  • 2.1. Poisson-Nernst-Planck equations: General electrostatic boundary conditions decompose the boundary into Dirichlet, Neumann, and Robin parts, with compatibility and gauge conditions required in relevant cases.
  • 2.1. Poisson-Nernst-Planck equations: The electrostatic free energy may require boundary correction terms under mixed conditions with nonzero prescribed boundary data.
  • 2.1. Poisson-Nernst-Planck equations: Under periodic or no-flux ionic boundary conditions, each species conserves mass, and nonnegative initial concentrations remain nonnegative.
  • 2.2. Modified Poisson–Nernst–Planck Models: Modified PNP models add steric interactions and concentration-gradient corrections to represent finite ion sizes, correlations, solvent occupancy, and strong concentration variations.
  • 2.2. Modified Poisson–Nernst–Planck Models: Because modification changes only the free-energy functional and chemical potentials, classical and modified models share the same constrained JKO framework.

3. Variational schemes for Poisson-Nernst-Planck models

The paper formulates PNP time steps as Wasserstein/JKO variational problems with transport and Poisson constraints, then discretizes them into convex programs. The construction supports general electrostatic boundary conditions and preserves key physical structure, although convergence to continuous PDE solutions remains open in some complex settings.

  • Variational formulation: PNP dynamics are formulated through a Wasserstein gradient-flow and JKO framework for the free energy.The construction uses transport variables for ionic concentrations and treats the Poisson equation as an additional coupling constraint.
  • Variational formulation: The semi-discrete JKO step combines Benamou–Brenier transport with the Poisson constraint while updating endpoint concentrations and potential together.The auxiliary transport variable is discretized with one step, consistent with the first-order accuracy of the outer JKO step.
  • Fully discrete scheme: The fully discrete formulation is a convex objective with linear PDE constraints for concentrations, fluxes, and electrostatic potential.The discretization uses cell-centered grids, centered-difference continuity equations, and boundary-aware discrete Laplace operators.
  • Scope and limitations: For multi-species PNP with anisotropic diffusion and complex boundary conditions, convergence of JKO solutions to the continuous PDE remains open.Well-posedness and weak convergence are available only under suitable assumptions cited in the paper.
  • Fully discrete scheme: The scheme accommodates classical and modified PNP models by changing the energy functional, including steric and concentration-gradient terms.Natural boundary conditions are imposed for concentration gradients, while mixed Dirichlet–Neumann conditions are retained in the electrostatic energy.
  • Structure preservation: The discrete variational scheme preserves original energy dissipation, mass conservation of p and n, and nonnegativity of both concentrations.These properties are established at the discrete level for the fully discrete JKO construction.

the discrete Poisson constraint, the optimality of uk+1 yields

The discrete optimality argument yields the principal structure-preserving properties of the scheme: energy dissipation, conservation of ionic masses, and nonnegative concentrations.

  • Energy dissipation: The discrete energy inequality reproduces the original energy dissipation law.The action terms are nonnegative, so the next-step discrete free energy does not exceed the previous-step value.
  • Mass conservation: Summing the discrete continuity equation with no-flux conditions proves conservation of p mass and, analogously, n mass.The proof applies the boundary conditions to cancel the discrete flux contributions.
  • Positivity: The admissible set enforces p_i,j ≥ 0 and n_i,j ≥ 0, so the JKO minimizer remains nonnegative.Positivity follows directly from the definition of the transport action and the admissible set.

4. Primal–dual splitting methods

The paper develops preconditioned and transformed primal–dual methods for the linearly constrained convex JKO problem. Their efficiency depends on conditioning the coupled saddle-point system and solving the global dual subproblem with structure-exploiting methods.

  • Problem formulation: The fully discrete JKO scheme is cast as a linearly constrained convex problem Au = b and relaxed into a proximal optimization formulation.The relaxation parameter δ permits a small constraint violation while enabling efficient proximal algorithms.
  • PrePD: Direct PD3O-type discretizations suffer severe step-size restrictions and slow convergence for dynamic JKO minimization.This motivates preconditioned primal–dual methods.
  • PrePD: PrePD uses block-diagonal primal and dual preconditioners, with Tv = AAT, to improve saddle-point conditioning and accelerate convergence.The resulting primal update is local for the transport action, while the dual update contains the main global linear solve.
  • VPTPD: VPTPD combines a Schur-complement block-triangular transformation with variable-dependent preconditioning.Its design aims to make the transformed saddle-point system nearly upper triangular and retain pointwise separability of the primal proximal step.
  • VPTPD: The primal proximal operator separates across species and grid cells when the preconditioner is diagonal, enabling inexpensive parallel computation.Local action subproblems reduce to scalar cubic equations, solved by closed-form formulas or Newton iterations.
  • Fast dual solvers: PNP dual updates are the computational bottleneck because the dual operator couples the Poisson constraint with both ionic continuity equations.The paper therefore develops fast solvers exploiting the coupled block structure under different potential boundary conditions.

5. Fast solvers for the coupled dual subproblem

The coupled dual system is solved either by block Gauss–Seidel or Schur-complement PCG, with transform-based and structured linear-algebra routines exploiting the PDE blocks. Schur-PCG is more robust for strong coupling, while both methods can achieve O(N log N) iterations without assembling the full matrix.

  • Coupled dual system: The constraint operator couples continuity-equation divergence blocks, density restrictions, and a Poisson Laplace block under specified boundary conditions.The block structure separates transport and electrostatic components while retaining their coupling through the linear constraints.
  • Dual solvers: BGS decouples the variables through stationary block updates, whereas Schur-PCG eliminates ionic dual variables and applies PCG to a reduced electrostatic system.BGS has lower cost per inner iteration, while Schur-PCG provides Krylov acceleration on the reduced system.
  • Solver comparison: Schur-PCG generally becomes more robust and efficient than BGS in strongly coupled or ill-conditioned regimes, although BGS is faster for moderate permittivity.The Schur reduction removes explicit block coupling, while BGS sweeps increase significantly as ϵ decreases.
  • Computational cost: O(N log N) cost per BGS sweep or Schur-PCG iteration avoids storing the assembled full block matrix.The fast actions use elliptic-operator inversions accelerated by FFT-based, sparse Cholesky, PCG, or multigrid solvers depending on the framework.
  • Transform-based realizations: 1D and multidimensional tensor-product operators are diagonalized by discrete sine or cosine transforms, enabling transform-space division with special handling of Neumann null modes.Homogeneous Neumann conditions require care because the corresponding operator has a null mode.
  • Transform-based realizations: DCT- and DST-based inversions efficiently apply the transport and Poisson block inverses for Dirichlet boundary conditions without assembling the coupled matrix.For PrePD, the transport block uses DCT-based fast algorithms and the Poisson-related block uses DST-based algorithms.

6. Numerical results

The numerical study validates the Poisson-constrained JKO scheme, compares primal–dual methods and dual solvers, and examines classical, modified, and small-permittivity PNP behavior.

  • Study scope: The experiments assess scheme accuracy, structure preservation, comparisons with existing numerical and optimization methods, and primal–dual algorithm performance.The study also investigates ionic interaction phenomena through numerical tests.
  • Study scope: Modified PNP experiments illustrate how concentration-gradient energy and spatial interactions affect the computed solutions.The modified-model tests extend the validation beyond classical PNP.

6.1. Validation tests

Validation tests show first-order temporal accuracy and preservation of energy dissipation, mass conservation, and concentration positivity. Small-permittivity simulations reproduce diffuse-charge asymptotics and retain these properties while remaining computationally reliable.

  • Benchmark experiments for accuracy: First-order temporal accuracy is observed for p, n, and ϕ in the one-dimensional Poisson-constrained JKO scheme.The numerical solution at t = 0.1 is compared with a fine-reference computation.
  • Benchmark experiments for accuracy: The discrete total energy decays monotonically, relative mass errors remain controlled, and ionic concentration extrema confirm positivity preservation.These observations support energy dissipation, mass conservation, and nonnegativity in the one-dimensional test.
  • Modified PNP tests: Concentration-gradient energy penalizes sharp spatial variations and produces smoother modified-PNP concentration profiles.The effect is shown for different values of the gradient-correction strength σ.
  • Diffuse-charge dynamics: In weak-voltage thin-double-layer tests, enrichment and depletion layers form near electrodes while the bulk remains nearly electroneutral and field-free.Most of the voltage drop becomes confined to diffuse layers.
  • Diffuse-charge dynamics: The numerical cathodic-charge curves match the analytical asymptotic solution for ϵD ∈ {0.1, 0.01, 0.001}, with smaller ϵD better capturing long-time equilibrium.The scheme also preserves energy dissipation, positivity, and mass conservation at ϵ = 2 × 10^-6.
  • Diffuse-charge dynamics: In the strongly nonlinear regime, the potential agrees with the Gouy–Chapman composite profile, while salt depletion and localized interfacial charge reproduce the predicted bulk response.The interior remains approximately electroneutral as excess salt and opposite diffuse charge localize near the electrodes.

6.2. Comparison tests

Comparison tests show that the JKO and primal–dual approaches remain stable and efficient in challenging small-permittivity settings. Schur-PCG and inexact dual solves improve robustness or reduce cost, while VPTPD accelerates outer convergence relative to PrePD.

  • Comparison with existing numerical methods: For ϵ = 0.0025, PJM fails to maintain stability under the tested discretization, whereas JKO preserves stability and energy dissipation.PJM can handle this regime only with a smaller time step and finer mesh at similar CPU cost.
  • Comparison with existing numerical methods: The JKO scheme is more robust in small-permittivity regimes without substantially increasing computational time despite requiring minimization at each time step.The comparison is made against PJM in the two-dimensional Neumann-boundary test.
  • Comparison with existing numerical methods: PrePD requires fewer iterations and less total CPU time than AEPG for moderate ϵ, while AEPG fails to converge within the iteration limit at ϵ = 0.0025.PrePD still converges and preserves monotone discrete-energy decay in the small-permittivity test.
  • Fixed-charge tests: As ϵ decreases, concentration and potential profiles become sharper and more localized near fixed charges, consistent with thinner electrostatic screening layers.Mobile-charge redistribution partially compensates fixed charge and localizes the potential variation.
  • Performance of BGS and Schur-PCG dual solvers: Schur-PCG becomes more robust and efficient than BGS for small ϵ, although BGS is faster at moderate ϵ.The outer primal–dual iteration count is the same for both exact dual solvers, but their inner solver costs differ.
  • Inexact dual solves: Using one inexact inner dual iteration preserves the outer primal–dual iteration count even for small ϵ and substantially reduces total CPU time.Exact and inexact solvers eventually follow the same convergence trajectory for one JKO step.
  • Primal–dual method comparison: VPTPD and VPTPD(λ) reduce primal–dual iterations and total CPU time relative to PrePD, with adaptive stepsizes providing additional efficiency.The reduction in outer iterations outweighs the increase in inner BGS sweeps, especially in the tested small-ϵ comparison.

6.3. Extended experiments for modified PNP models

Extended modified-PNP experiments show how concentration-gradient regularization and electrostatic effects reshape ionic distributions over time and across parameter settings. The three-dimensional simulation exhibits smoothing, anisotropic redistribution, complementary ion profiles, and boundary-dominated late-time states.

  • 2D parameter study: Concentration-gradient energy produces smoother, more spatially coherent ionic profiles in the modified PNP model.As the coefficient increases, local variations are suppressed and profiles become less sensitive to individual fixed-charge interfaces.
  • 2D parameter study: As the concentration-gradient coefficient increases, the electrostatic potential changes more moderately because it is nonlocally coupled to smoothed charge density.
  • 3D evolution: At early times, the three-dimensional ionic concentrations occupy two separated spherical regions with high-concentration zones at opposite positions.
  • 3D evolution: Diffusion and concentration-gradient regularization smooth sharp interfaces, while localized fixed charge drives anisotropic ionic redistribution.
  • 3D evolution: Opposite electrostatic drift directions produce increasingly complementary ion profiles as initially localized structures expand, deform, and spread through the domain.
  • 3D evolution: By t = 0.015, ionic distributions become boundary-dominated while the potential approaches a smooth quasi-steady configuration governed by fixed charge and mixed boundary conditions.

7. Conclusion

The paper proposes a unified Poisson-constrained JKO framework and efficient primal–dual solvers for strongly coupled PNP models. The scheme preserves key variational structures and performs robustly across experiments, especially at small dielectric permittivity, while remaining first-order in time.

  • Scheme and algorithms: The unified scheme extends the JKO framework for Wasserstein gradient flows to strongly coupled multi-variable PNP models.
  • Scheme and algorithms: Two efficient primal–dual splitting algorithms, PrePD and VPTPD, use fast dual solvers for the resulting constrained optimization problems.
  • Structure-preserving properties: The Poisson equation is incorporated as an additional linear constraint, enabling preservation of energy dissipation, ionic positivity, and mass conservation.
  • Numerical performance: The algorithms show superior robustness and efficiency in experiments with varied boundary conditions and fixed charges, especially for small dielectric permittivity.
  • Scope and outlook: The proposed scheme is first-order accurate in time, with higher-order variational structure identified as future work.
Loading 2608.30693v1…