PaperlayerView exact revision on arXiv

Source-linked AI summary

A linear, second-order, unconditionally energy-stable scalar-auxiliary-variable scheme for the CHQ viscoelastic phase-separation model

Ritesh Bhola, Michael te Vrugt

arXiv:2610.30077v1math.NAcond-mat.softcs.NA

TL;DR

Viscoelastic phase separation couples phase evolution and elastic stress through a cross-diffusion flux, while the reference discretisation requires a nonlinear solve. This paper develops a linear SAV–BDF2 discretisation with extrapolated phase-dependent coefficients. It reports an exact discrete energy identity and lower solver cost than an independent nonlinear Newton–BDF2 solver.

  • Problem

    The phase and elastic-stress equations are coupled through cross-diffusion, and energy-stable schemes for uncoupled Cahn–Hilliard equations do not automatically extend to this system.

    Source: §1.1 · p. 1

  • Method

    The method combines SAV treatment of the double-well potential with BDF2 and extrapolates phase-dependent coefficients to keep each time step linear.

    Source: §1.1 · p. 2, §3.2 · p. 6

  • Results

    2.5–3.0× fewer linear iterations and 1.7–2.0× less wall time than an independent nonlinear Newton–BDF2 solver; the method also achieves second-order temporal convergence.

    Source: §8 · p. 17

  • Takeaways & Limitations

    The discrete energy identity also holds for the spectral-filtering variant used to produce the physics results, establishing that method’s stability.

    Source: §8 · p. 17

    The modified discrete energy differs from the true energy by O(∆t^2), so the true energy may show O(∆t^2) non-monotonicities even as the modified energy decreases monotonically.

    Source: §4 · p. 9

Abstract

from arXiv · show

The CHQ model describes viscoelastic phase separation in a polymer--solvent mixture: a Cahn--Hilliard equation for the polymer volume fraction $\phi$ is two-way coupled, through a phase-dependent cross-diffusion flux, to a damped relaxation equation for a scalar bulk elastic stress $q$. We propose a scalar-auxiliary-variable (SAV) time discretisation of this system, combining the SAV reformulation of the double-well potential with a second-order backward differentiation formula (BDF2) and a second-order extrapolation of all phase-dependent mobilities: a linear, computationally cheaper complement to the reference scheme's fully nonlinear Newton solve. The resulting method is (i) linear at every time step, with a coupling structure we prove has exactly one solution at every step, for every choice of time-step size; (ii) mass conserving, exactly, as an algebraic identity of the discrete equations themselves (the only drift ever seen numerically traces to the finite tolerance of the iterative linear solver, not a failure of the conservation property); and (iii) unconditionally stable for a modified discrete energy, via an exact one-step energy--dissipation identity. We show the energy identity extends to this filtered variant as well, with two explicit extra non-negative terms. Hence, the scheme actually producing the physics results is itself provably stable for the same energy. A Fourier--Galerkin implementation verifies that the energy identity holds at every time step to within floating-point round-off. It also confirms second-order temporal accuracy, including in the steep-interface regime, and agreement with an independently implemented nonlinear Newton--BDF2 solver. Compared with that independent nonlinear Newton--BDF2 solver, at equal grid resolution, the SAV method uses $2.5$--$3.0\times$ fewer linear-solver iterations and $1.7$--$2.0\times$ less wall-clock time.

1 Introduction

The paper studies a viscoelastic phase-separation model whose phase and elastic-stress equations are coupled through a phase-dependent cross-diffusion flux. It develops a linear SAV–BDF2 scheme designed to preserve the coupled model’s energy structure while addressing the nonlinear solve required by the reference discretisation.

1.1 The model

  • The model couples a Cahn–Hilliard-type phase equation to elastic-stress evolution through a cross-diffusion flux, making the extension of energy-stable discretisations nontrivial.The reference discretisation for this reduced model requires a nonlinear solve at every time step.Source: §1.1 · pp. 1–2

1.2 Energy structure

  • The continuous energy contains interfacial, bulk free-energy, and elastic contributions, and its dissipation is non-negative under a coercivity condition on the model coefficients.The dissipation law provides the structure the discrete scheme is designed to reproduce.Source: §1.2 · p. 2

1.3 Related work

  • Prior work includes a fully nonlinear energy-identity scheme for this model, while the paper reports no earlier SAV, IEQ, or energy-quadratisation discretisation for this cross-diffusion structure.The authors describe their work as complementary to Ref. [6], not an improvement on its analysis, and do not claim a rigorous error estimate.Source: §1.3 · p. 3, §1.4 · p. 3

1.4 Contribution

  • The scheme combines SAV and BDF2 with a common extrapolated state for the phase-dependent coefficients, preserving the coupled cross-diffusion dissipation structure without a time-step restriction.Freezing the coefficients independently or at the new-time solution would break that structure and leave a time-step-dependent remainder.Source: §1.4 · p. 3
  • The authors report an exact discrete energy identity, unique per-step solvability, and numerical verification of second-order temporal accuracy and agreement with an independent nonlinear solver.The identity also holds for the post-step spectral-filtering variant used in the physics study.Source: §1.4 · p. 3

2 The continuous energy law and its SAV reformulation

The continuous model has an energy dissipation law under a coefficient coercivity condition. The SAV reformulation introduces a scalar to represent the nonlinear free-energy contribution, but its exact continuous-time relation to the phase field is not guaranteed by discrete time stepping.

2 The continuous energy law and its SAV reformulation

  • The continuous energy decreases at a rate given by the dissipation, which is non-negative when the coefficients satisfy the coercivity condition.The scheme is designed to reproduce this decay rate at the discrete level.Source: §2 · p. 4

2.1 The SAV reformulation

  • SAV introduces a scalar tracking the shifted nonlinear free-energy contribution, replacing the double-well derivative with a linear term and an additional scalar evolution equation.When the scalar matches the square root of the shifted nonlinear energy, the reformulated system recovers the original one; the quadratic elastic contribution needs no auxiliary variable.Source: §2.1 · p. 5
  • In continuous time the SAV scalar retains its defining relation to the phase field, but the discrete update need not preserve that relation exactly.Consequently, the modified and true discrete energies can differ by a small, controlled amount.Source: §2.1 · p. 5

3 The SAV–BDF2 scheme

The SAV–BDF2 scheme combines second-order time stepping with extrapolated, frozen coefficients to retain a linear formulation while preserving the coupled structure.

3.1 Notation

  • BDF2 approximates the new-time derivative from two preceding levels, while its second-difference quantity tracks the energy remainder rather than approximating a derivative.The BDF2 derivative is accurate to O(Δt^2).Source: §3.1 · p. 6

3.2 Freezing the coefficients

  • The extrapolated phase state approximates ϕ(tn+1) to O(Δt^2); clamping it to [δc, 1 − δc] guards against singular coefficient forms.Using ϕn instead would be only first-order accurate, while the clamp was inactive in the studied convergence regime.Source: §3.2 · p. 6
  • Freezing every phase-dependent coefficient at a shared extrapolated state keeps the equations linear while preserving the coupled cross-diffusion structure.The extrapolation is second-order accurate; the coercivity margin for the benchmark coefficients is εb, independent of the time step.Source: §3.2 · pp. 6–7

3 The SAV–BDF2 scheme

  • The scheme replaces time derivatives with BDF2, freezes phase-dependent coefficients, and adjoins the SAV auxiliary equation for r.It solves for ϕn+1, μn+1, qn+1, and rn+1 from the preceding two levels.Source: §3 · pp. 5–6, §3.3 · p. 7

3.3 The discrete scheme

  • Each equation is linear in the unknowns because coefficients use known frozen values, with the remaining coupling passing through the scalar rn+1.The default first step uses backward Euler; its local O(Δt^2) error does not degrade the second-order global rate in the reported check.Source: §3.3 · p. 7

4 The discrete energy identity

The scheme has an exact discrete energy–dissipation identity, so its modified energy decreases for every positive time-step size. The result also applies to the filtered scheme, while modified-energy stability does not imply monotone decay of the physical energy.

4 The discrete energy identity

  • The modified discrete energy decreases at every step for every Δt > 0, through an exact identity with non-negative dissipation and remainder terms.The identity is subject to the stated coefficient and spatial-discretisation assumptions.Source: §4 · p. 8
  • The decreasing modified energy bounds ∇ϕ, q, and r for all times and every Δt > 0; mass conservation and a bounded gradient also bound ϕ.The energy reaches ϕ through its gradient, so fixed mass is used to control its spatial average.Source: §4 · p. 9
  • The filtered scheme retains the energy identity with two additional non-negative penalties, making the physics-run scheme provably stable.The filter is an orthogonal, gradient-commuting projection; both added terms remained strictly positive with margin on real trajectories.Source: §4 · p. 9
  • Monotone decay applies to the modified energy, not necessarily the physical energy E(ϕn, qn), which may show O(Δt^2) non-monotonicities.The modified and physical energies differ by O(Δt^2).Source: §4 · p. 9

5 Structural properties

The scheme’s linear system is uniquely solvable at every step for every positive time-step size under the stated assumptions, and its equations conserve mass exactly.

5 Structural properties

  • Under the energy-identity assumptions and b⋆ ≥ b1 > 0, the linear system is uniquely solvable for every Δt > 0.The proof uses the energy identity on the homogeneous system to show its only solution is trivial.Source: §5 · p. 10
  • The discrete equations conserve the spatial average of ϕ exactly at every step, including the backward-Euler startup step.The proof tests the phase equation with the constant function and uses the resulting recursion for consecutive mass values.Source: §5 · p. 10

6 Accuracy and implementation

This section justifies the scheme’s second-order accuracy and specifies the implementation used in the paper’s computations, including its spatial discretisation and linear-solver choices.

Accuracy and implementation

  • The scheme is globally O(∆t^2): BDF2, extrapolated coefficients, and the SAV substitution each contribute O(∆t^2) consistency, with the clamp inactive in the studied convergence regime.The argument gives local O(∆t^3) error accumulated over O(1/∆t) steps; a rigorous relative-energy error estimate is left for future work.Source: §6 · pp. 10–11
  • Each step extrapolates coefficients, solves the linear system for (ϕ^n+1, µ^n+1, q^n+1, r^n+1), and may filter ϕ and q in long-time runs; the first step uses backward Euler.The first step skips extrapolation, and the optional filter applies the 2/3-rule.Source: §6 · p. 11
  • Periodic Fourier–Galerkin collocation supplies the discrete integration-by-parts and commutation properties required by the stability proof; zeroing Nyquist modes in the Laplacian is essential to preserve its identity.Using the naive −|k|^2 symbol makes the discrete gradient and Laplacian inconsistent at Nyquist modes, causing (12) to fail if those modes are present.Source: §6 · p. 11
  • The linear system is solved by a dense direct method on small grids or matrix-free GMRES with a Fourier-diagonal preconditioner on larger grids; an independent Newton solver cross-checks SAV effects.The direct and iterative methods agree to 10^-12 where both are affordable; the Newton reference uses the same spatial discretisation and treats f′(ϕ^n+1) fully implicitly.Source: §6 · p. 11
  • The reference configuration uses Ω = (0, 1)^2 and d0 = C0 = 1, and satisfies the model’s structural conditions.These parameters define the reference configuration used in the implementation discussion.Source: Table 1 · p. 12
  • The true-energy curves for the two finer ∆t values track closely until about t ≈ 3, while ∆t = 10^-2 diverges earlier and settles higher.The supplied figure caption reports the curve comparison but gives no numerical energy values.Source: Figure 1 · p. 12

7 Numerical results

Numerical tests verify energy decay and mass conservation, convergence, and agreement with an independent solver. The coupled-versus-uncoupled test also shows slower phase separation with viscoelastic coupling.

7.1 Headline demonstration: energy, mass, ∆t-robustness

  • At ∆t = 10−2, E(10) = 0.1199, roughly 2× the finer runs’ value; the finer runs give 0.0704 versus 0.0618, a 12% gap.The coarse step diverges from the finer runs earlier, despite the energy identity holding throughout; the remaining finer-step gap is unresolved.Source: §7.1 · p. 13

7.2 The energy identity, convergence, and independent solvers

  • Across random frozen states and ∆t from 10−4 to 10^2, the N = 13 linear-system condition number stays below 2 × 10^6.The smallest singular value remains σmin ≳4 × 10−4 in those tested cases.Source: §7.2 · p. 13
  • Dense-solve energy-identity residuals stay at ∼10−14–10−17 per step for both resolved smooth and under-resolved noisy initial data.In the noisy case, the clamp is active at no more than 9% of grid points; mass stays at round-off in these dense-solve runs.Source: §7.2 · p. 14, Figure 3 · p. 14
  • Temporal tests give pt ≈2.00 for backward-Euler startup, 1.68 to 1.95 for Richardson startup as ∆t shrinks, and physical-case self-convergence of 1.98 and 2.00.The smooth-data spatial test’s error falls by a factor of 107 as N increases from 16 to 96.Source: §7.2 · p. 14, Figure 4 · p. 15
  • Across four successive halvings, the SAV–Newton max-norm ϕ difference shrinks at orders 2.04, 2.02, 2.00, and 2.01.The difference between their true energies also shrinks at a matching O(∆t^2) rate.Source: §7.2 · p. 15, Figure 5 · p. 16
  • 2.5–3.0× fewer linear iterations and 1.7–2.0× less wall time for SAV than Newton–BDF2 at N = 64–128, with substantially less memory.The comparison uses γ = 0.20, where both schemes compute the same solution.Source: §7.2 · p. 15
  • Sweeping εq to zero preserves the energy identity and monotonicity, moving the solution by only 4 × 10−6; varying C0 over five orders changes ϕ, q by at most 8 × 10−11.The tests check that stress smoothing and the SAV constant do not materially alter the computed solution.Source: §7.2 · p. 16

7.3 Coupled-vs-uncoupled demonstration

  • At T = 10, the coupled model releases 59% of its initial energy versus 83% for the A ≡0 control and remains on a downward energy slope.Mass is conserved to ∼5 × 10−12 in both runs, and the coupled model needs a median 333 GMRES iterations/step versus 16 for the control.Source: §7.3 · p. 16

8 Conclusion

The paper derives and verifies a linear, second-order, unconditionally energy-stable SAV–BDF2 scheme for the CHQ model. Its numerical results include improved solver costs and viscoelastic suppression of phase separation.

8 Conclusion

  • An exact discrete energy–dissipation identity establishes unconditional stability, and the per-step linear system is uniquely solvable for every ∆t > 0.The identity also holds for the post-step spectral-filtering variant used in the physics study.Source: §8 · pp. 16–17
  • The SAV scheme uses 2.5–3.0× fewer linear iterations and 1.7–2.0× less wall time than an independent nonlinear Newton–BDF2 solver.It also achieves second-order temporal convergence, including for the thin interface γ = 10−3.Source: §8 · p. 17
  • On the reference configuration, the coupled system forms a persistent, still-evolving gel-like network, while the uncoupled control forms droplets that coarsen to a fully phase-separated state.The morphology panels show the two models already differ qualitatively within the first time unit.Source: §8 · p. 17, Figure 6 · p. 17
  • The energy–dissipation identity holds numerically to round-off at every time step, including for εq = 0 and across five orders of magnitude in C0.The method also achieves spectral accuracy in space for the smooth test problems considered.Source: §8 · p. 17

A Proofs and derivations

The appendix provides derivations and proofs for the SAV reformulation and the scheme’s continuous and discrete stability properties, including unique solvability of the linear scheme for every time step and arbitrary prior data. It also extends the modified-energy stability identity to the filtered scheme, with two additional nonnegative terms.

Derivation conventions

  • Periodic-domain integration by parts removes boundary terms in the test-function identities, while extrapolated coefficients are fixed independently of the unknowns and keep each time step linear.The BDF2 quantity Φ combines two time levels to produce a telescoping identity.Source: A Proofs and derivations · pp. 20–21

Exact equivalence of the SAV reformulation

  • The SAV reformulation is exactly equivalent to the original system because r remains equal to the shifted potential energy, making rG(ϕ) = f′(ϕ).The SAV evolution equation and the shifted potential-energy expression have the same time derivative and initial value; conversely, the original solution defines r by that expression.Source: A Proofs and derivations · p. 21

Proof of Proposition 2.1

  • The continuous energy satisfies dE/dt = −D, with non-negative dissipation revealed by completing the square in the coupled flux terms.The coupling sign makes the cross term appear twice when the phase and stress equations are tested, allowing the square completion.Source: A Proofs and derivations · p. 23

Proof of the G-stability identity (13)

  • The algebraic BDF2 G-stability identity rewrites tested time differences as telescoping Φ differences plus a non-negative remainder.It applies to any three elements of the same space and is identified as the Nevanlinna–Odeh multiplier for BDF2, used for G-stability.Source: A Proofs and derivations · pp. 21, 24

Proof of Theorem 4.1

  • Testing the discrete equations, completing the square, and applying G-stability yields the one-step energy identity (12).The shared coupling term appears twice in the phase and stress equations; the spatial terms become non-negative dissipation, while the BDF2 terms telescope.Source: A Proofs and derivations · pp. 24, 26

Energy stability proof

  • The unfiltered energy proof telescopes the BDF2 terms and yields the stability identity because every dissipation term is nonnegative.All coefficients on the right-hand side are frozen independently of ϕn+1 and qn+1, so the identity applies exactly.Source: A Proofs and derivations · pp. 26–27

Filtered-scheme stability extension

  • The filtered scheme satisfies the modified-energy identity with two additional nonnegative dissipative terms, so the stored energy does not increase for any Δt > 0.Filtering lowers the energy by 5γ/4 ∥(I − F)∇ϕn+1∥2 + 5/4∥(I − F)qn+1∥2.Source: A Proofs and derivations · p. 28

Proof of Theorem 5.1

  • The scheme has a unique solution for every Δt > 0 and arbitrary prior data, provided b⋆ ≥ b1 > 0.The proof shows the homogeneous system has only the trivial solution; the positive lower bound on b⋆ is used to conclude ∇ν = 0.Source: A Proofs and derivations · p. 29

B The manufactured solution used in Section 7.2

The manufactured-solution check uses explicit smooth fields for ϕ and q, then adds forcing terms so those fields solve the modified system exactly. The resulting run measures temporal truncation-error order, not stability.

B The manufactured solution used in Section 7.2

  • The manufactured fields are ϕex = 0.5 + 0.15 cos(2πx) cos(2πy) e−3t and qex = 0.10 sin(2πx) sin(2πy) e−2t on Ω = (0, 1)2.The test uses γ = 0.30, chosen so ϕ = 0.5 is not spinodally unstable for the low wavenumber.Source: B The manufactured solution used in… · p. 29
  • Forcing terms equal to the residuals of the two evolution equations make the manufactured fields an exact solution of the forced system.The forcing is evaluated at tn+1 and added to (S1) and (S3), while (S4) is left unchanged.Source: B The manufactured solution used in… · p. 30
  • Because the manufactured fields do not solve the unforced system, this test checks truncation-error order rather than the energy identity or stability.The exact solution at each ∆t allows direct measurement of ∥ϕn − ϕex(tn)∥∞.Source: B The manufactured solution used in… · pp. 29–30
Loading 2610.30077v1…