pith. sign in

arxiv: 2508.19965 · v3 · submitted 2025-08-27 · 🧮 math.NA · cs.NA

High-order nonuniform time-stepping and MBP-preserving linear schemes for the time-fractional Allen-Cahn equation

Pith reviewed 2026-05-18 21:10 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords time-fractional Allen-Cahn equationmaximum-bound principleenergy stabilitylinear stabilized schemesnonuniform time-steppingL1 schemeL2-1σ scheme
0
0 comments X

The pith

New linear schemes for the time-fractional Allen-Cahn equation preserve both discrete energy stability and the maximum-bound principle under nonuniform time steps.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops high-order linear stabilized schemes for the time-fractional Allen-Cahn equation that maintain both discrete energy stability and the maximum-bound principle on nonuniform time grids. A prediction strategy produces a second-order solution that respects the bound, while a newly introduced nonnegative auxiliary functional supplies the stabilization term needed to control the explicit nonlinear potential. Two resulting schemes are shown to be unconditionally energy stable; the L1 version preserves the discrete maximum-bound principle for arbitrary steps, and the L2-1σ version does so under a mild restriction. An improved L2-1σ variant replaces the balanced stabilization with an unbalanced term that exploits the auxiliary functional's boundedness and monotonicity, extending reliable MBP preservation to larger steps. Numerical tests confirm that the methods retain accuracy, energy decay, and bound preservation.

Core claim

By pairing a new prediction strategy that yields a second-order MBP-preserving solution with an essential nonnegative auxiliary functional, the authors construct L1 and L2-1σ linear stabilized schemes that are unconditionally energy stable for the time-fractional Allen-Cahn equation; the L1 scheme preserves the discrete maximum-bound principle without time-step restrictions, while the L2-1σ scheme requires only a mild restriction, and an improved L2-1σ version further relaxes this restriction through an unbalanced stabilization term.

What carries the argument

The combination of a prediction strategy that delivers a second-order MBP-preserving solution and a nonnegative auxiliary functional that enables a stabilization term to dominate the predicted nonlinear potential.

If this is right

  • The L1 scheme guarantees discrete maximum-bound preservation for any positive time steps.
  • The L2-1σ scheme preserves the bound once the time steps satisfy a mild explicit restriction.
  • The improved unbalanced L2-1σ scheme extends reliable MBP preservation to substantially larger time steps.
  • Both schemes remain unconditionally energy stable on arbitrary nonuniform grids.
  • High-order accuracy is retained while the physical bounds and energy dissipation are preserved.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The auxiliary functional may transfer to other time-fractional phase-field models that obey a maximum principle.
  • Nonuniform grids paired with these schemes could enable efficient adaptive refinement near initial singularities typical of fractional equations.
  • The prediction-plus-auxiliary approach might reduce the need for iterative nonlinear solvers in long-time fractional simulations.

Load-bearing premise

The auxiliary functional remains nonnegative, bounded, and monotonic so that the stabilization term can always dominate the predicted nonlinear contribution and enforce the discrete bound.

What would settle it

A numerical run of the L1 scheme on a sequence of successively refined nonuniform grids that produces a discrete solution violating the maximum bound at some time step.

Figures

Figures reproduced from arXiv: 2508.19965 by Bingyin Zhang, Hongfei Fu.

Figure 1
Figure 1. Figure 1: The initial highly-oscillated phase field [PITH_FULL_IMAGE:figures/full_fig_p014_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Time evolutions of the maximum norm and energy of si [PITH_FULL_IMAGE:figures/full_fig_p014_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Time evolutions of the maximum norm and energy of si [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: The dynamic snapshots of the numerical solution [PITH_FULL_IMAGE:figures/full_fig_p015_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Time evolutions of the maximum norm (left), energy [PITH_FULL_IMAGE:figures/full_fig_p016_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: The dynamic snapshots of the numerical solution [PITH_FULL_IMAGE:figures/full_fig_p017_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: Time evolutions of the maximum norm (left), energy [PITH_FULL_IMAGE:figures/full_fig_p017_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: The dynamic snapshots of the numerical solution [PITH_FULL_IMAGE:figures/full_fig_p018_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Time evolutions of the energy (left) and time steps [PITH_FULL_IMAGE:figures/full_fig_p018_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Plots of the iso-surfaces (value 0) of the numerical solution at different time instants obtained by the L1-sESAV scheme for the tFAC model with α = 0.5 0 20 40 60 80 100 time 0 0.2 0.4 0.6 0.8 1 1.2 maximum norm 0 50 100 0.9999999999996 0.9999999999997 0.9999999999998 0.9999999999999 1 (a) maximum-norm of φ 0 20 40 60 80 100 time 0 0.005 0.01 0.015 0.02 0.025 0.03 energy (b) energy 0 20 40 60 80 100 time… view at source ↗
Figure 11
Figure 11. Figure 11: Time evolutions of the maximum norm (left), energ [PITH_FULL_IMAGE:figures/full_fig_p019_11.png] view at source ↗
Figure 12
Figure 12. Figure 12: Time evolutions of the maximum norm and energy of s [PITH_FULL_IMAGE:figures/full_fig_p027_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: Time evolutions of the maximum norm and energy of s [PITH_FULL_IMAGE:figures/full_fig_p027_13.png] view at source ↗
Figure 14
Figure 14. Figure 14: The dynamic snapshots of the numerical solution [PITH_FULL_IMAGE:figures/full_fig_p028_14.png] view at source ↗
Figure 15
Figure 15. Figure 15: Time evolutions of the maximum norm (left), energ [PITH_FULL_IMAGE:figures/full_fig_p029_15.png] view at source ↗
Figure 16
Figure 16. Figure 16: The dynamic snapshots of the numerical solution [PITH_FULL_IMAGE:figures/full_fig_p029_16.png] view at source ↗
Figure 17
Figure 17. Figure 17: Time evolutions of the maximum norm (left), energ [PITH_FULL_IMAGE:figures/full_fig_p030_17.png] view at source ↗
read the original abstract

In this paper, we present a class of nonuniform time-stepping, high-order linear stabilized schemes that can preserve both the discrete energy stability and maximum-bound principle (MBP) for the time-fractional Allen-Cahn equation. To this end, we develop a new prediction strategy to obtain a second-order and MBP-preserving predicted solution, which is then used to handle the nonlinear potential explicitly. Additionally, we introduce an essential nonnegative auxiliary functional that enables the design of an appropriate stabilization term to dominate the predicted nonlinear potential, and thus to preserve the discrete MBP. Combining the newly developed prediction strategy and auxiliary functional, we propose two unconditionally energy-stable linear stabilized schemes, L1 and L2-$1_\sigma$ schemes. We show that the L1 scheme unconditionally preserves the discrete MBP, whereas the L2-$1_\sigma$ scheme requires a mild time-step restriction. Furthermore, we develop an improved L2-$1_\sigma$ scheme with enhanced MBP preservation for large time steps, achieved through a novel unbalanced stabilization term that leverages the boundedness and monotonicity of the auxiliary functional. Representative numerical examples validate the accuracy, effectiveness, and physics-preserving of the proposed methods.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 2 minor

Summary. The paper develops high-order linear stabilized schemes with nonuniform time-stepping for the time-fractional Allen-Cahn equation. It introduces a prediction strategy yielding a second-order MBP-preserving predictor and a nonnegative auxiliary functional used to construct stabilization terms that dominate the explicit nonlinear potential. Two schemes are analyzed: an L1 scheme claimed to be unconditionally energy-stable and MBP-preserving, and an L2-1σ scheme that is unconditionally energy-stable but requires a mild time-step restriction for MBP preservation, together with an improved L2-1σ variant using unbalanced stabilization for large steps. Numerical examples are provided to illustrate accuracy and preservation properties.

Significance. If the energy-stability and MBP-preservation proofs hold for arbitrary nonuniform meshes, the work would supply practical high-order methods for long-time simulation of fractional phase-field models while respecting key structural properties. The auxiliary-functional approach and predictor construction offer a template that may extend to other time-fractional nonlinear PDEs.

major comments (2)
  1. [Section 3.3] The proof that the auxiliary functional remains bounded and monotonic independently of the maximum step-size ratio is not supplied in sufficient detail. Without a uniform bound, the algebraic domination argument used to conclude |u^{n+1}| ≤ 1 for the L1 scheme (the inequality immediately after the definition of the stabilization term) may fail on meshes with large consecutive step ratios, undermining the unconditional MBP claim.
  2. [Theorem 4.2] Theorem 4.2 asserts unconditional energy stability for both schemes, yet the proof relies on the same auxiliary functional whose monotonicity under nonuniform steps is not verified independently of the step-ratio bound; this makes the energy-stability result conditional on the same unproven estimate required for MBP preservation.
minor comments (2)
  1. [Section 3.4] The definition of the L2-1σ kernel and the precise form of the unbalanced stabilization term in the improved scheme should be restated explicitly rather than referenced only by equation number.
  2. [Section 5] In the numerical section, the convergence tables do not report the observed maximum step-size ratios used in the nonuniform meshes; adding this information would help readers assess the practical range of the claimed robustness.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for the thorough review and valuable comments on our manuscript arXiv:2508.19965. The points raised highlight areas where the exposition of the auxiliary functional's properties can be strengthened. We address each major comment below and will incorporate the necessary clarifications and expanded proofs in the revised version.

read point-by-point responses
  1. Referee: [Section 3.3] The proof that the auxiliary functional remains bounded and monotonic independently of the maximum step-size ratio is not supplied in sufficient detail. Without a uniform bound, the algebraic domination argument used to conclude |u^{n+1}| ≤ 1 for the L1 scheme (the inequality immediately after the definition of the stabilization term) may fail on meshes with large consecutive step ratios, undermining the unconditional MBP claim.

    Authors: We agree that the current exposition in Section 3.3 would benefit from a more detailed and self-contained proof of the boundedness and monotonicity of the auxiliary functional for arbitrary step-size ratios. In the revision, we will expand this section to include a complete derivation: starting from the definition of the auxiliary functional, we will explicitly show its nonnegativity, derive a uniform upper bound independent of the maximum ratio r = max(τ_{k+1}/τ_k), and prove monotonicity using the discrete fractional integral properties and the second-order MBP-preserving predictor. This will directly support the algebraic domination argument immediately following the stabilization term definition, thereby rigorously justifying the unconditional MBP preservation for the L1 scheme on general nonuniform meshes. revision: yes

  2. Referee: [Theorem 4.2] Theorem 4.2 asserts unconditional energy stability for both schemes, yet the proof relies on the same auxiliary functional whose monotonicity under nonuniform steps is not verified independently of the step-ratio bound; this makes the energy-stability result conditional on the same unproven estimate required for MBP preservation.

    Authors: We acknowledge the logical dependence noted by the referee. The energy stability proof in Theorem 4.2 indeed invokes the monotonicity and boundedness properties of the auxiliary functional. By supplying the expanded, self-contained proof in the revised Section 3.3 (as described above), which establishes these properties independently of any step-ratio restriction, the energy stability argument will hold unconditionally for arbitrary nonuniform time steps. We will also restructure the proof of Theorem 4.2 to explicitly reference the newly detailed estimates, ensuring the unconditional claims for both the L1 and L2-1σ schemes are fully substantiated. revision: yes

Circularity Check

0 steps flagged

No significant circularity detected in derivation chain

full rationale

The paper introduces a new prediction strategy for a second-order MBP-preserving predictor and a nonnegative auxiliary functional whose boundedness and monotonicity are used to construct a stabilization term that dominates the explicit nonlinear potential. These elements are presented as novel contributions enabling the proofs of unconditional energy stability for both schemes and unconditional MBP preservation for the L1 scheme. No quoted step reduces the target preservation claims to a fitted parameter, self-definition, or prior self-citation by construction; the algebraic inequalities for |u^{n+1}| ≤ 1 follow from the stated properties of the auxiliary functional on the nonuniform mesh rather than from renaming or tautological re-use of the result itself. The derivation therefore remains self-contained against external benchmarks.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 1 invented entities

The work rests on standard properties of fractional derivatives and discrete operators, plus the new auxiliary functional introduced to enable stabilization.

axioms (1)
  • standard math Standard properties of the Caputo fractional derivative and summation-by-parts formulas for discrete fractional operators hold.
    Invoked to derive the schemes and prove stability and preservation.
invented entities (1)
  • Nonnegative auxiliary functional no independent evidence
    purpose: To construct a stabilization term that dominates the predicted nonlinear potential and enforces the discrete MBP.
    Newly introduced in the paper; no independent evidence outside the construction is provided.

pith-pipeline@v0.9.0 · 5747 in / 1421 out tokens · 41620 ms · 2026-05-18T21:10:22.884725+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

What do these tags mean?
matches
The paper's claim is directly supported by a theorem in the formal canon.
supports
The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends
The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses
The paper appears to rely on the theorem as machinery.
contradicts
The paper's claim conflicts with a theorem or certificate in the canon.
unclear
Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

46 extracted references · 46 canonical work pages

  1. [1]

    Ainthworth and Z

    M. Ainthworth and Z. Mao. Analysis and approximation of a fractional Cahn– Hilliard equation. SIAM J. Numer. Anal. , 55:1689–1718, 2017

  2. [2]

    Akagi, G

    G. Akagi, G. Schimperna, and A. Segatti. Fractional Cahn –Hilliard, Allen–Cahn and porous medium equations. J. Differ. Equ. , 261:2935–2985, 2016

  3. [3]

    Alikhanov

    A. Alikhanov. A new difference scheme for the time fractio nal diffusion equation. J. Comput. Phys. , 280:424–438, 2015

  4. [4]

    Allen and J

    S. Allen and J. Cahn. A microscopic theory for antiphase b oundary motion and its application to antiphase domain coarsening. Acta Metall. , 27:1085–1095, 1979

  5. [5]

    Anderson, G

    D. Anderson, G. McFadden, and A. Wheeler. Diffuse-interf ace methods in fluid mechanics. Annu. Rev. Fluid Mech. , 30:139–165, 1998

  6. [6]

    Cahn and J

    J. Cahn and J. Hilliard. Free energy of a nonuniform syste m I: Interfacial free energy. J. Chem. Phys. , 28:258–267, 1958

  7. [7]

    L. Chen, J. Zhao, W. Cao, H. Wang, and J. Zhang. An accurate and efficient algorithm for the time-fractional molecular beam epitaxy m odel with slope selection. Comput. Phys. Commun. , 245:106842, 2019. 32

  8. [8]

    Q. Du, L. Ju, X. Li, and Z. Qiao. Maximum principle preserv ing exponential time differencing schemes for the nonlocal Allen–Cahn equation. SIAM J. Numer. Anal. , 57:875–898, 2019

  9. [9]

    Q. Du, L. Ju, X. Li, and Z. Qiao. Maximum bound principles f or a class of semilinear parabolic equations and exponential time-differencing sch emes. SIAM Rev., 63:317– 359, 2021

  10. [10]

    Q. Du, C. Liu, and X. Wang. A phase field approach in the num erical study of the elastic bending energy for vesicle membranes. J. Comput. Phys. , 198:450–468, 2004

  11. [11]

    Q. Du, J. Yang, and Z. Zhou. Time-fractional Allen–Cahn equations: Analysis and numerical methods. J. Sci. Comput. , 85:42, 2020

  12. [12]

    Hawkins-Daarud, K

    A. Hawkins-Daarud, K. Zee, and J. Oden. Numerical simul ation of a thermodynam- ically consistent four-species tumor growth model. Int. J. Numer. Methods Biomed. Eng., 8:3–24, 2012

  13. [13]

    D. Hou, L. Ju, and Z. Qiao. A linear second-order maximum bound principle- preserving BDF scheme for the Allen–Cahn equation with a gen eral mobility. Math. Comp., 92:2515–2542, 2023

  14. [14]

    Hou and C

    D. Hou and C. Xu. Highly efficient and energy dissipative s chemes for the time fractional Allen–Cahn equation. SIAM J. Sci. Comput. , 43:A3305–A3327, 2021

  15. [15]

    T. Hou, T. Tang, and J. Yang. Numerical analysis of fully discretized Crank-Nicolson scheme for fractional-in-space Allen–Cahn equations. J. Sci. Comput. , 72:1–18, 2017

  16. [16]

    Huang, Z

    Q. Huang, Z. Qiao, and H. Yang. Maximum bound principle a nd non-negativity preserving ETD schemes for a phase field model of prostate can cer growth with treatment. Comput. Methods Appl. Mech. Engrg. , 426:116981, 2024

  17. [17]

    M. Inc, A. Yusuf, A. Aliyu, and D. Baleanu. Time-fractio nal Cahn–Allen and time- fractional Klein–Gordon equations: lie symmetry analysis , explicit solutions and convergence analysis. Physica A Stat. Mech. Appl. , 493:94–106, 2018

  18. [18]

    B. Ji, H. Liao, Y. Gong, and L. Zhang. Adaptive second-or der Crank–Nicolson time-stepping schemes for time-fractional molecular beam epitaxial growth models. SIAM J. Sci. Comput. , 42:B738–B760, 2020

  19. [19]

    B. Ji, H. Liao, and L. Zhang. Simple maximum principle pr eserving time-stepping methods for time-fractional Allen–Cahn equation. Adv. Comput. Math. , 46:37, 2020

  20. [20]

    B. Jin, R. Lazarov, and Z. Zhou. An analysis of the L1, sch eme for the subdiffusion equation with nonsmooth data. IMA J. Numer. Anal. , 36:197–221, 2016

  21. [21]

    L. Ju, X. Li, and Z. Qiao. Generalized SA V-exponential i ntegrator schemes for Allen–Cahn type gradient flows. SIAM J. Numer. Anal. , 60:1905–1931, 2022

  22. [22]

    L. Ju, X. Li, and Z. Qiao. Stabilized exponential-SA V sc hemes preserving energy dissipation law and maximum bound principle for the Allen–C ahn type equations. J. Sci. Comput. , 92:66, 2022. 33

  23. [23]

    Z. Li, H. Wang, and D. Yang. A space-time fractional phas e-field model with tunable sharpness and decay behavior and its efficient numerical simu lation. J. Comput. Phys., 347:20–38, 2017

  24. [24]

    H. Liao, D. Li, and J. Zhang. Sharp error estimate of the n onuniform L1 formula for linear reaction-subdiffusion equations. SIAM J. Numer. Anal. , 56:1112–1133, 2018

  25. [25]

    H. Liao, W. Mclean, and J. Zhang. A second-order scheme w ith nonuniform time steps for a linear reaction-subdiffusion problem. Commun. Comput. Phys. , 30:567– 601, 2021

  26. [26]

    H. Liao, T. Tang, and T. Zhou. A second-order and nonunif orm time-stepping maximum-principle preserving scheme for time-fractional Allen–Cahn equations. J. Comput. Phys. , 414:109473, 2020

  27. [27]

    H. Liao, T. Tang, and T. Zhou. An energy stable and maximu m bound preserving scheme with variable time steps for time fractional Allen–C ahn equation. SIAM J. Sci. Comput. , 43:A3503–A3526, 2021

  28. [28]

    H. Liao, T. Tang, and T. Zhou. Positive definiteness of re al quadratic forms resulting from the variable-step L1-type approximations of convolution operators. Sci. China Math., 67:237–252, 2024

  29. [29]

    H. Liao, X. Zhu, and H. Sun. Asymptotically compatible e nergy and dissipation law of the nonuniform L2-1σ scheme for time fractional Allen–Cahn model. J. Sci. Comput., 99:46, 2024

  30. [30]

    Lin and C

    Y. Lin and C. Xu. Finite difference/spectral approximat ions for the time-fractional diffusion equation. J. Comput. Phys. , 225:1533–1552, 2007

  31. [31]

    H. Liu, A. Cheng, H. Wang, and J. Zhao. Time-fractional A llen–Cahn and Cahn– Hilliard phase-field models and their numerical investigat ion. Comput. Math. Appl. , 76:1876–1892, 2019

  32. [32]

    Mustapha, B

    K. Mustapha, B. Abdallah, and K. Furati. A discontinuou s Petrov-Galerkin method for time-fractional diffusion equations. SIAM J. Numer. Anal. , 52:2512–2529, 2014

  33. [33]

    Qi and X

    R. Qi and X. Zhao. A unified design of energy stable scheme s with variable steps for fractional gradient flows and nonlinear integro-differe ntial equations. SIAM J. Sci. Comput. , 46:A130–A155, 2024

  34. [34]

    T. Qian, X. Wang, and P. Sheng. Molecular scale contact l ine hydrodynamics of immiscible flows. Phys. Rev. E , 68:016306, 2003

  35. [35]

    Z. Qiao, Z. Zhang, and T. Tang. An adaptive time-steppin g strategy for the molec- ular beam epitaxy models. SIAM J. Sci. Comput. , 33:1395–1414, 2011

  36. [36]

    Quan and B

    C. Quan and B. Wang. Energy stable L2 schemes for time-fr actional phase-field equations. J. Comput. Phys. , 458:111085, 2022

  37. [37]

    Shen and J

    J. Shen and J. Xu. Convergence and error analysis for the scalar auxiliary variable (SA V) schemes to gradient flows. SIAM J. Numer. Anal. , 56:2895–2912, 2018. 34

  38. [38]

    J. Shen, J. Xu, and J. Yang. The scalar auxiliary variabl e (SA V) approach for gradient flows. J. Comput. Phys. , 353:407–416, 2018

  39. [39]

    Stynes, E

    M. Stynes, E. O’Riordan, and J. Gracia. Error analysis o f a finite difference method on graded meshes for a time-fractional diffusion equation. SIAM J. Numer. Anal. , 55:1057–1079, 2017

  40. [40]

    Sun and X

    Z. Sun and X. Wu. A fully discrete difference scheme for a d iffusion-wave system. Appl. Numer. Math. , 56:193–209, 2006

  41. [41]

    Tang and J

    T. Tang and J. Yang. Implicit-explicit scheme for the Al len–Cahn equation preserves the maximum principle. J. Comput. Math. , 34:471–481, 2016

  42. [42]

    T. Tang, H. Yu, and T. Zhou. On energy dissipation theory and numerical stability for time-fractional phase-field equations. SIAM J. Sci. Comput. , 41:A3757–A3778, 2019

  43. [43]

    S. Wise, J. Lowengrub, H. Frieboes, and V. Cristini. Thr ee-dimensional multispecies nonlinear tumor growth I: Model and numerical method. J. Theor. Biol. , 253:524– 543, 2008

  44. [44]

    Xu and T

    C. Xu and T. Tang. Stability analysis of large time-step ping methods for epitaxial growth models. SIAM J. Numer. Anal. , 44:1759–1779, 2006

  45. [45]

    Zhang, C

    G. Zhang, C. Huang, A. Alikhanov, and B. Yin. A high-orde r discrete energy decay and maximum-principle preserving scheme for time fraction al Allen–Cahn equation. J. Sci. Comput. , 96:39, 2023

  46. [46]

    Zhang and Z

    Z. Zhang and Z. Qiao. An adaptive time-stepping strateg y for the Cahn–Hilliard equation. Commun. Comput. Phys. , 11:1261–1278, 2012. 35