REVIEW 3 major objections 4 minor 39 references
A positivity-preserving hybrid DDG method for Poisson--Nernst--Planck systems
T0 review · 3 major / 4 minor · reviewed 2026-08-07 · deepseek-v4-flash
Pith's one-line read A modified numerical flux provably keeps ion concentrations positive in arbitrary high-order discontinuous Galerkin schemes for the Poisson–Nernst–Planck system.
desk verdict Solid 1D positivity theorem; multi-D proof has a Gauss-Lobatto exactness gap that undermines the claimed higher-dimensional result. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The central object is the modified DDG normal-derivative flux $\widehat{\partial_x p_i} = \hat{\partial}_x p_i + \frac{\tilde{\beta}_i}{2}[c_i]$, where $\tilde{\beta}_i = |\hat{\partial}_x p_i|/\{c_i\}$ when the interface average is positive. This local correction makes every boundary term in the Gauss-Lobatto quadrature representation of the updated cell average non-negative, so the average becomes a sum of positive quantities. The companion positivity-preserving reconstruction (3.19) rescales each polynomial toward its cell average so its minimum stays above a small threshold without changing the average.
What would settle it
A single run of scheme (3.15) with $c^n_{ih}>0$ that yields a non-positive cell average $\bar{c}^{n+1}_{ij}$ while the time step obeys (3.18) would refute Theorem 3.1.
Extended reading notes
Core claim
The central claim is Theorem 3.1: for the fully discrete DDG scheme (3.15), if the current concentration polynomial is pointwise positive and the time step satisfies (3.18), then the next cell average $\bar{c}^{n+1}_{ij}$ is positive. The proof rewrites the cell-average update as a Gauss-Lobatto quadrature and shows that the modified flux $\widehat{\partial_x p_i} = \widehat{\partial_x p_i} + \frac{\tilde{\beta}_i}{2}[c_i]$ with $\tilde{\beta}_i = |\widehat{\partial_x p_i}|/\{c_i\}$ makes all interface coefficients non-negative. The same argument is carried to tensor-product rectangular meshes in Theorem 4.1. To obtain pointwise positivity, the paper applies a positivity-preserving reconstruction limiter that keeps each polynomial above a small threshold $\delta$ without changing its cell average and without degrading accuracy when $\delta < h^{k+1}$. The result is an arbitrary high-order, provably positivity-preserving DDG method for the Poisson–Nernst–Planck system.
Load-bearing premise
The proof assumes the positivity-preserving limiter applied at each step keeps the concentration polynomial pointwise positive and the time-step condition valid, an interaction analyzed only in one dimension.
Editorial extensions
If this is right
- In one dimension, the modified flux plus forward Euler guarantees positive cell averages for any polynomial degree $k$, as long as the time step obeys (3.18).
- When combined with the reconstruction limiter, the scheme yields pointwise positive concentrations, keeping the logarithmic term $\log c_i$ well-defined at every step.
- On rectangular meshes in several dimensions, the same positivity propagation holds for tensor-product $Q^k$ elements under the multi-dimensional CFL condition (4.28).
- Numerical tests show that the modified flux does not degrade the expected $k+1$ order of accuracy, and that mass conservation and free-energy dissipation are observed in the hybrid scheme.
Reading between the lines
- Beyond the paper: the local, data-dependent choice of $\tilde{\beta}_i$ should transfer to other gradient-flow systems written as $\partial_t c = \nabla\cdot(c\nabla p)$ with $p=\delta E/\delta c$.
- Beyond the paper: the multi-dimensional theorem certifies only the modified-flux update, so pointwise positivity in 2D and 3D still relies on the limiter interaction the paper does not formalize.
- Beyond the paper: the numerical evidence for energy dissipation is not accompanied by a proof for the hybrid algorithm; proving dissipation for the modified flux would close the remaining structural gap.
- Beyond the paper: the $1/d$ factor in the multi-dimensional CFL bound makes high-order SSP-RK stages cheaper but time-step heavy, so an implicit or locally varying variant is a natural practical follow-up.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript proposes a hybrid direct discontinuous Galerkin (DDG) method for the Poisson-Nernst-Planck system that enforces positivity of ion concentrations. The main ingredients are: (i) a locally modified numerical flux, defined in (1.3) and (3.16)-(3.17), which makes the cell-average update nonnegative under a forward Euler discretization and a suitable CFL condition; (ii) a positivity-preserving reconstruction limiter, (3.19), applied at each step to obtain pointwise positivity; and (iii) an extension of the scheme to rectangular meshes in multiple dimensions, with a claimed cell-average positivity theorem (Theorem 4.1). Numerical experiments in 1D and 2D demonstrate high-order accuracy, positive cell averages, mass conservation, and free-energy dissipation, although the latter two properties are explicitly stated to be only numerically observed, not proven.
Significance. If the multi-dimensional positivity result were rigorously established, the paper would deliver a genuinely arbitrary-high-order, positivity-preserving DG scheme for PNP systems, which is a valuable contribution given the difficulty of maintaining positivity in high-order discretizations. The 1D argument is elegant and internally consistent: the modified flux is designed so that the Gauss-Lobatto decomposition of the cell average has nonnegative coefficients exactly when the CFL condition (3.18) holds. The numerical evidence for accuracy and positivity is strong, including a comparison showing failure of the unmodified scheme from [24]. However, the proof of the multi-dimensional extension contains a serious quadrature gap, and the hybrid limiter is not analyzed in multiple dimensions, so the central multi-D claim is not yet supported.
major comments (3)
- [§4.2, proof of Theorem 4.1] The proof is not complete as written. After deriving the update for the cell average, the authors apply the M-point Gauss-Lobatto quadrature rule in all directions except x^l to an integrand that contains the modified flux (4.25)-(4.26). Because beta-tilde_i^l = |\partial_{x_l} p_i^n| / c_i^n, the integrand involves the absolute value of a polynomial in the transverse variables and is not a polynomial. Gauss-Lobatto quadrature is therefore not exact for that integrand, and a positive quadrature sum does not imply that the exact cell average is positive. The stated CFL condition (4.28) controls |\partial_{x_l} p_i^n| only at the quadrature nodes, but a polynomial can attain a much larger maximum away from those nodes; for example, the transverse profile of a degree-3 element can peak at x=0.3 while the 3-point Gauss-Lobatto nodes are 0, 0.5, and 1. To conclude positivity of the exact integral, one would need the integrand to be nonnegative pointwise, which requires a pointwise (supremum) CFL bound rather than a node-only bound. Thus Theorem 4.1 is not established as stated. The one-dimensional Theorem 3.1 is unaffected by this issue.
- [Section 4 / Algorithm 1] The paper claims that the hybrid algorithm provides pointwise positivity in higher dimensions, but the positivity-preserving reconstruction (3.19) and Algorithm 1 are formulated for one-dimensional P^k elements. Section 4 introduces Q^k tensor-product elements and proves, even modulo the quadrature gap, only cell-average positivity for the modified flux; it does not define the reconstruction limiter on Q^k elements or analyze its effect on the multi-D update. Consequently, the pointwise-positive claim for the multi-D hybrid algorithm is not supported, even if the quadrature issue in Theorem 4.1 is repaired.
- [§3.3, Algorithm 1 and Remark 3.1] The algorithm as stated does not guarantee that the CFL hypothesis (3.18) holds at every step. Step 3 modifies c_i^n by the limiter, which changes p_i^n through (3.15b) and therefore changes \partial_x p_i^n; the practical time-step rule (3.21) contains an empirical factor gamma and is not proven to imply (3.18) for the reconstructed data. To claim that the hybrid algorithm preserves cell-average positivity, the authors should either prove a bound for the limiter's effect on \partial_x p_i^n, or explicitly state that (3.18) is verified after reconstruction and the time step is reduced when necessary.
minor comments (4)
- [Equations (3.18) and (4.28)] The notation |\partial_x p|^{-1} should be accompanied by the convention that this quantity is infinite when |\partial_x p| = 0; this is done implicitly in the proof of Theorem 3.1 for boundary edges but should be stated uniformly for all interior edges and for the multi-D condition (4.28).
- [Remark 3.2 and Section 6] Remark 3.2 states that discrete mass conservation and free-energy dissipation are not formally proven. Mass conservation for the modified scheme actually follows by summing (3.15a) over cells because the interface flux is single-valued; adding this one-line argument would strengthen the paper and sharpen the remark.
- [Table 5.5] For the case (m, gamma) = (3, 1/4), the reported convergence orders are non-monotone and the error increases from N=30 to N=40; the authors should comment on whether this is an artifact of the time-step rule or an indication that the CFL condition is not being enforced as required.
- [Labels in Theorem 2.1 and Lemma 3.1] The displayed labels in Theorem 2.1 and Lemma 3.1 contain corrupted text ('tieacce...lowercase'); these should be cleaned in the final version.
Circularity Check
No circularity found: the positivity theorem is a constructive consequence of the deliberately modified flux, and the self-citations to prior work are not load-bearing.
full rationale
I find no significant circularity. The central result, Theorem 3.1, is a direct theorem about the fully discrete scheme (3.15), in which the modified flux parameter is deliberately defined by (3.17) as beta_i = |d_dx p_i^n| / {c_i^n}. The proof uses exactly this choice to rewrite the interface terms as (d_dx p_i^n + beta_i {c_i^n}) c_i and (-d_dx p_i^n + beta_i {c_i^n}) c_i, each of which is nonnegative by the inequality a + |a| >= 0. This is a constructive design choice rather than a fitted parameter later renamed as a prediction; no quantity asserted as an output is fed back as an input. The pointwise positivity of the hybrid algorithm is enforced by the limiter (3.19), while the theorem supplies the cell-average positivity needed to keep the limiter applicable; this is a standard induction on the scheme, not a tautology. The paper cites the authors' earlier framework [24] for the DDG formulation, Theorem 2.1, and Lemma 3.1, but these are independent mathematical statements with stated assumptions that do not include the positivity result, so they do not make the argument circular. I also flag two limitations that are correctness concerns rather than circularity: Remark 3.2 and the concluding section explicitly state that discrete energy dissipation is not formally proven, and the proof of Theorem 4.1 applies Gauss-Lobatto quadrature in the transverse directions to integrands involving beta_i^l c_i with beta_i^l = |d_xl p_i^n| / c_i^n, which is not a polynomial in the transverse variables; thus the multi-D cell-average positivity theorem is not established as written. These observations do not undermine the 1D Theorem 3.1 and are not circular steps.
Assumptions & free parameters
free parameters (3)
- gamma (CFL safety factor) =
1/4, 1/2, 1 in Table 5.5
- delta (limiter threshold) =
not reported; constrained delta < h^{k+1}
- DDG flux coefficients beta0 and beta1 =
e.g., (2, -), (4, 1/20), (3, -), (9, 1/12), (19, 1/12), (16, 1/6)
assumptions (4)
- standard math Gauss-Lobatto quadrature with M >= (k+3)/2 points is exact for the k-th degree polynomial cell average and has positive weights.
- domain assumption The mesh is regular, h_max <= C h_min, and the DDG flux scales as O(h^{-1}).
- domain assumption The accuracy-preserving reconstruction lemma (Lemma 3.1, quoted from [24]) holds and the limiter preserves cell averages.
- standard math SSP-Runge-Kutta time discretizations are positive linear combinations of forward Euler steps and therefore inherit the positivity bound.
Cite this review
Pith. "Pith review of A positivity-preserving hybrid DDG method for Poisson--Nernst--Planck systems." pith.science (2026). https://pith.science/paper/MPHXYMFL
@misc{pith2026250524058,
author = {Pith},
title = {Pith review of: A positivity-preserving hybrid DDG method for Poisson--Nernst--Planck systems},
year = {2026},
howpublished = {\url{https://pith.science/paper/MPHXYMFL}},
note = {Machine review of arXiv:2505.24058}
}
read the original abstract
In earlier work [H. Liu and Z. Wang, J. Comput. Phys., 328(2017)], an arbitrary high-order conservative and energy-dissipative direct discontinuous Galerkin (DDG) scheme was developed. Although this scheme enforced solution positivity using cell averages as reference values, it lacked a theoretical guarantee for the positivity of those cell averages. In this study, we develop a novel arbitrary high-order DDG method with rigorously proven positivity-preserving properties. Specifically, the positivity of the cell averages is ensured through a modified numerical flux in combination with forward Euler time discretization. To achieve point-wise positivity of ion concentrations, we introduce a hybrid algorithm that integrates a positivity-preserving limiter. The proposed method is further extended to higher-dimensional problems with rectangular meshes. Numerical results confirm the scheme's high-order accuracy, guaranteed positivity preservation, and consistent discrete energy dissipation.
Figures
Figures from the paper (8 more)
Reference graph
Works this paper leans on
- [24]
-
[25]
H. Liu, Z. Wang, P. Yin, and H. Yu. Positivity-preserving third order DG schemes for Poisson-Nernst-Planck equations. J. Comput. Phys., 452:110777, 2022
work page 2022
-
[3]
J. A. Carrillo, H. Liu, and H. Yu. Positivity-preserving and energy-dissipating discontinuous Galerkin methods for non- linear nonlocal Fokker-Planck equations.Commun. Appl. Ind. Math., 2025
work page 2025
-
[1]
J.-D. Benamou and Y. Brenier. A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math., 84(3):375–393, 2000
work page 2000
-
[2]
J. A. Carrillo, A. Chertock, and Y. H. Huang. A finite-volume method for nonlinear nonlocal equations with a gradient flow structure.Commun. Comput. Phys., 17(1):233–258, 2015
work page 2015
-
[4]
H. Chen, H. Liu, and X. Xu. The Onsager principle and structure preserving numerical schemes.J. Comput. Phys., 523:113679, 2025
work page 2025
-
[5]
Datta.Electronic Transport in Mesoscopic Systems
S. Datta.Electronic Transport in Mesoscopic Systems. Cambridge University Press, 1997
work page 1997
-
[6]
J. Ding, Z. Wang, and S. Zhou. Positivity preserving finite difference methods for Poisson-Nernst-Planck equations with steric interactions: application to slit-shaped nanopore conductance.J. Comput. Phys., 397:108864, 2019
work page 2019
Show all 39 references
-
[7]
J. Ding, Z. Wang, and S. Zhou. Structure-preserving and efficient numerical methods for ion transport.J. Comput. Phys., 418:109597, 2020
2020
-
[8]
Eisenberg and W
B. Eisenberg and W. Liu. Poisson-Nernst-Planck systems for ion channels with permanent charges.SIAM J. Math. Anal., 38:1932–1966, 2007
1932
-
[9]
Glasstone.An introduction to Electrochemstry
S. Glasstone.An introduction to Electrochemstry. Van Nostrand Company, Inc., Princeton, NJ, 1942
1942
-
[10]
Gottlieb, C
S. Gottlieb, C. Shu, and E. Tadmor. Strong stability-preserving high-order time discretization methods.SIAM Rev., 43(1):89–112, 2001
2001
-
[11]
D. He, K. Pan, and X. Yue. A positivity preserving and free energy dissipative difference scheme for the Poisson-Nernst- Planck system.J. Sci. Comput., 81:436–458, 2019
2019
-
[12]
Hille.Ion channels of excitable membranes, volume 507
B. Hille.Ion channels of excitable membranes, volume 507. Sinauer Sunderland, MA, 2001. 21 Figure 5.10.The contours evolution ofc 1 −1/30,c 2 −2/30 andψ
2001
-
[13]
Hu and X
J. Hu and X. Huang. A fully discrete positivity-preserving and energy-dissipative finite difference scheme for Poisson- Nernst-Planck equations.Numer. Math., 145:77–115, 2020
2020
-
[14]
Jerome.Analysis of Charge Transport: A Mathematical Study of Semiconductor Devices
J. Jerome.Analysis of Charge Transport: A Mathematical Study of Semiconductor Devices. Springer, Berlin, 1996
1996
-
[15]
Jordan, D
R. Jordan, D. Kinderlehrer, and F. Otto. The variational formulation of the Fokker–Planck equation.SIAM J. Math. Anal., 29(1):1–17, 1998
1998
-
[16]
Kinderlehrer, L
D. Kinderlehrer, L. Monsaingeon, and X. Xu. A wasserstein gradient flow approach to Poisson-Nernst-Planck equations. ESAIM Control Optim. Calc. Var., 23(1):137–164, 2017
2017
-
[17]
Li.Electrokinetics in Microfluidics
D. Li.Electrokinetics in Microfluidics. Academic Press, 2004
2004
-
[18]
C. Liu, C. Wang, S. M. Wise, X. Yue, and S. Zhou. A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst-Planck system.Math. Comp., 90:2071–2106, 2021
2021
-
[19]
C. Liu, C. Wang, S. M. Wise, X. Yue, and S. Zhou. A second order accurate, positivity preserving numerical method for the Poisson-Nernst-Planck system and its convergence analysis.J. Sci. Comput., 97(23), 2023
2023
-
[20]
H. Liu. Analysis of direct discontinuous Galerkin methods for multi-dimensional convection–diffusion equations.Numer. Math., 147(4):839–867, 2021
2021
-
[21]
Liu and W
H. Liu and W. Maimaitiyiming. Unconditional positivity-preserving and energy stable schemes for a reduced Poisson- Nernst-Planck system.Commun. Comput. Phys., 27:1505–1529, 2020
2020
-
[22]
Liu and W
H. Liu and W. Maimaitiyiming. Efficient, positive, and energy stable schemes for multi-D Poisson-Nernst-Planck systems. J. Sci. Comput., 87:92, 2021
2021
-
[23]
Liu and W
H. Liu and W. Maimaitiyiming. A dynamic mass transport method for Poisson-Nernst-Planck equations.J. Comput. Phys, 473:111699, 2023. 22
2023
-
[26]
Liu and J
H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) methods for diffusion problems.SIAM J. Numer. Anal, 47:675–698, 2009
2009
-
[27]
Liu and J
H. Liu and J. Yan. The direct discontinuous Galerkin (DDG) method for diffusion with interface corrections.Commun. Comput. Phys., 8(3):541–564, 2010
2010
-
[28]
Liu and H
H. Liu and H. Yu. The entropy satisfying discontinuous Galerkin method for Fokker–Planck equations.J. Sci. Comput., 62:803–830, 2015
2015
-
[29]
Markowich.The Stationary Semiconductor Device Equations
P. Markowich.The Stationary Semiconductor Device Equations. Springer-Verlag, New York, 1986
1986
-
[30]
P. A. Markowich, C. A. Ringhofer, and C. Schmeiser.Semiconductor Equations. Springer-Verlag Inc, New York, 1990
1990
-
[31]
Metti, J
M. Metti, J. Xu, and C. Liu. Energetically stable discretizations for charge transport and electrokinetic models.J. Comput. Phys., 306:1–18, 2016
2016
-
[32]
M. S. Mock.Analysis of Mathematical Models of Semiconductor Devices, volume 3. Boole Press, 1983
1983
-
[33]
Shen and J
J. Shen and J. Xu. Unconditionally positivity preserving and energy dissipative schemes for Poisson–Nernst–Planck equations.Numer. Math., 148:671–697, 2021
2021
-
[34]
Su and H
S. Su and H. Tang. A positivity-preserving and free energy dissipative hybrid scheme for the poisson-nernst-planck equations on polygonal and polyhedral meshes.Comput. Math. Appl., 108:33–48, 2022
2022
-
[35]
Z. Sun, J. A. Carrillo, and C.-W. Shu. A discontinuous galerkin method for nonlinear parabolic equations and gradient flow problems with interaction potentials.J. Comp. Phys., 352:76–104, 2018
2018
-
[36]
P. Yin, Y. Huang, and H. Liu. An iterative discontinuous Galerkin method for solving the nonlinear Poisson-Boltzmann equation.Commun. Comput. Phys., 16(2):491–515, 2014
2014
-
[37]
P. Yin, Y. Huang, and H. Liu. Error estimates for the iterative discontinuous Galerkin method to the nonlinear Poisson- Boltzmann equation.Commun. Comput. Phys., 23(1):168–197, 2018
2018
-
[38]
X. Zhang. On positivity-preserving high order discontinuous galerkin schemes for compressible Navier-Stokes equations. J. Comput. Phys., 328:301–343, 2017
2017
-
[39]
Zhang and C.-W
X. Zhang and C.-W. Shu. On maximum-principle-satisfying high order schemes for scalar conservation laws.J. Comput. Phys., 229:3091–3120, 2010. † Mathematics Department, Iowa State University, Ames, IA 50011, USA Email address:hliu@iastate.edu ‡ Department of Mathematics and St...
2010
Reviewed August 7, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.