Improving the Accuracy and Consistency of the Scalar Auxiliary Variable (SAV) Method with Relaxation

Many physical problems, such as interface dynamics[1, 20, 40], crystallization[9, 28], thin films[15, 33], polymers[14, 8], and liquid crystallization[12, 13] could be modeled by gradient flow systems which also agree with the second law of thermodynamics. If the total free energy is known, the gradient flow model could be obtained according to the mobility and the variation of free energy. Because of the nonlinear terms in the governing equation, neither the exact solution nor the numerical solution is easy to obtain. In general, consider the spatial-temporal domain $\Omega_{t}:=\Omega\times(0,T]$. The dissipative dynamics of the state variable $\phi$ are driven by

where $\mathcal{G}$ is a semi-positive definite operator known as the mobility operator, and $\mathcal{E}$ is a functional of $\phi$ known as the free energy. The triplet $(\phi,\mathcal{G},\mathcal{E})$ uniquely defines the dissipative system (gradient flow dynamics). For instance, given $\mathcal{G}=1$ and $\mathcal{E}=\displaystyle\int_{\Omega}\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}+\frac{1}{4}(\phi^{2}-1)^{2}d\mathbf{x}$, with $\displaystyle\varepsilon$ as a model parameter, we obtain the following Allen-Cahn equation

If we consider $\mathcal{G}=-\Delta$ and $\mathcal{E}=\displaystyle\int_{\Omega}\frac{\varepsilon^{2}}{2}|\nabla\phi|^{2}+\frac{1}{4}(\phi^{2}-1)^{2}d\mathbf{x}$, we obtain the Cahn-Hilliard equation

Given $\mathcal{G}=1$ and $\mathcal{E}=\displaystyle\int_{\Omega}D|\nabla\phi|^{2}d\mathbf{x}$, we have the heat equation

where $D$ is the diffusion coefficient.

All these models discussed above have an energy dissipation property. Mainly,

given the boundary terms are diminished to zero. Here we have used the inner producut notation $\displaystyle(f,g)=\int_{\Omega}fgd\mathbf{x}$, $\forall f,g\in L^{2}(\Omega)$. Numerical algorithms that solve such models shall also preserve the energy dissipation structure, i.e., follow the thermodynamic physical laws. Numerical schemes that preserve the energy dissipation structure are known as energy stable schemes. And if such structure-preserving doesn’t depend on the time step, the numerical schemes are known as unconditionally energy stable.

Many energy stable numerical schemes are proposed to approach the solutions of gradient flow models or the dissipative systems. The classical approaches are the fully implicit schemes [11]. Though some of them are unconditionally energy stable, solving such fully implicit schemes is not trivial. Nonlinear problems have to be solved in each time step. However, the existence and uniqueness of the solution usually have strong restrictions on the time steps which prevent those fully implicit numerical schemes from being widely used. One remedy is the convex splitting method [10], which splits the nonlinear terms of free energy into the subtraction of two convex functions. It is easy to check the convex-splitting schemes are unconditional energy stable and uniquely solvable. However, the general type of second-order convex splitting schemes is not available. So far, it is only possible to design second-order convex-splitting schemes case-by-case[34, 30, 3, 33]. Meanwhile, there are many other unconditional energy stable methods, such as stabilization method[28, 36], exponential time discretization method [35, 24, 6]. The stabilization method represents the nonlinear terms explicitly and adds some regularization terms to relax strict constraints for the time steps. Similarly, with the convex splitting method, it is usually limited to first-order accuracy. The exponential time discretization (ETD) method shows high-order accuracy by integrating the governing equation over a single time step and uses polynomial interpolations for the nonlinear terms. But the theoretical proofs for energy stability properties of high-order ETD schemes are still missing.

Recently, the numerical method named invariant energy quadratization(IEQ) or energy quadratization(EQ) are proposed [38, 37, 39, 41, 42, 43, 18, 17]. It is a generalization of the method of Lagrange multipliers or auxiliary variables from [2, 19]. The IEQ approach permits us to construct linear, second-order, unconditional energy stable schemes, and furthermore arbitrarily high-order unconditionally energy stable schemes[17, 16]. With many advantages of the IEQ or EQ approach, it usually leads to a coupled system with time-dependent coefficients. As a remedy, the SAV approach [31, 32, 25, 26, 27, 21, 29] have been proposed by introducing scalar auxiliary variables instead of auxiliary function variables. The SAV method also can be applied to a large class of gradient flow systems, which keeps the advantages of the EQ approach but usually leads to decoupled systems with constant coefficients. These properties make the SAV method easier to implement, so it is highly efficient. Besides, when the researchers applied the SAV approach to many different systems, several modified schemes were developed. Multiple scalar auxiliary variable(MSAV) approach [7] was proposed to solve the phase-field vesicle membrane model where two auxiliary variables introduced to match two additional penalty terms enforcing the volume and surface area. If using the introduced scalar variable to control both the nonlinear and the explicit linear terms, one high efficient SAV approach was developed[22], which spend half of time compared with the original SAV approach while keeping all its other advantages. One stabilized-scalar auxiliary variable(S-SAV)[23] approach was proposed to solve the phase-field surfactant model, which is a decoupled scheme and allowed to be solved step by step. For phase-field surfactant model, the authors in [44] also presented certain subtle explicit-implicit treatments for stress and convective terms to construct the linear, decoupled, unconditional energy stable schemes based on the classical SAV approach.

However, there is still a big gap for the IEQ or SAV method, making them not as perfect as expected. Mainly, these two methods preserve a ”modified” energy law according to the auxiliary variables instead of the original variables. This inconsistency introduces errors during the computation. In the end, even though the IEQ or SAV schemes preserve the ”modified” energy law, they are not necessarily preserving the original energy law, i.e., the energy law for the original PDE models might be violated by the numerical solutions. This is known in the community, but so far, no good remedy is available yet. This motivates our research in this paper. With a novel relaxation step, we effectively penalize the inconstancy between numerical solutions for the auxiliary variables and their continuous definitions. Thus we name this new approach as the relaxed-SAV (RSAV) method. Through the RSAV method, we are able to design novel linear, second-order, unconditionally energy stable schemes, which keep the advantage of the baseline SAV method and preserve the original energy laws. It turns out that the relaxation approach effectively improves the accuracy and consistency of the SAV method noticeably.

The rest of this paper is organized as follows. In Section 2, we revisit the baseline SAV approach and multiple SAV approaches for the general gradient flow system. Section 3 proposes the relaxed-SAV and relaxed-MSAV schemes. The energy stability properties of the RSAV and RMSAV methods are proved rigorously. Then, in Section 4, several specific examples and numerical tests are provided to verify the accuracy and effectiveness of the proposed relaxed SAV numerical schemes. In the end, we give a brief conclusion.

Consider the general gradient flow model

where $\phi$ is the state variable, $\mathcal{E}$ is the free energy, and $\mathcal{G}$ is a semi-positive definite operator for dissipative systems (and a skew-symmetric operator for reversible systems or Hamiltonian systems). In the rest of this paper, we consider periodic boundary conditions for simplicity, though all our results can be applied to models with thermodynamically consistent boundary conditions.

For the general gradient flow model (2.1), it has the following energy law

When $\mathcal{G}$ is semi-positive definite, we have $\frac{d\mathcal{E}}{dt}\leq 0$, and when $\mathcal{G}$ is a skew-symmetric operator, we have $\frac{d\mathcal{E}}{dt}=0$. Here we use the notation, $(f,g)=\displaystyle\int_{\Omega}fgd\mathbf{x}$, $\forall f,g\in L^{2}(\Omega)$. The induced norm will be denoted as $\|f\|=\sqrt{(f,f)}$.

Notice the definition in (2.6) tells us $\frac{q(t)}{Q(\phi)}=1$. Hence, we can observe the two energies (2.3) and (2.8) are equivalent in the PDE level. Meanwhile, the two PDE models (2.4) and (2.7) are equivalent. However, after temporal discretization, the numerical results of $q(t)$ and $Q(\phi)$ are not equal anymore, which means the discrete energies of (2.3) and (2.8) are not necessarily equivalent anymore. The major issue is that $q^{n+1}$ is no longer equal to $Q(\phi^{n+1})$ numerically. Thus, an energy stable scheme that satisfies the modified energy law in (2.8) does not necessarily satisfy the original energy law in (2.3).

To fix the inconsistency issue for $q^{n+1}$ and $Q(\phi^{n+1})$ (that are supposed to be equal as introduced in (2.6)), we propose a relaxation technique to penalize the difference between $q^{n+1}$ and $Q(\phi^{n+1})$. As will be clear in the following sections, the RSAV method introduces negligible extra computational cost, but it inherits all the baseline SAV method’s good properties.