Supplementary Variable Method for Developing Structure-Preserving Numerical Approximations to Thermodynamically Consistent Partial Differential Equations

Nonequilibrium phenomena require dynamical models derived from laws and principles of nonequilibrium thermodynamics to describe. The laws of thermodynamics, especially, the second law of thermodynamics or the equivalent generalized Onsager principle are fundamental laws/principles for developing such models for nonequilibrium phenomena not far from equilibria [12, 13, 16]. Assuming the thermodynamic variables describing nonequilibrium phenomena of a certain system are denoted as $\Phi\in{\bf R}^{n}$, where $n$ is the number of the variables, and the free energy of the system is described by a functional $F[\Phi]$ in isothermal cases (for nonisothermal cases, we use an entropy functional instead). The Onsager linear response theory yields the dynamical equation

where $\mathbf{M}$ and $\mathbf{R}$ are operators, $\frac{\delta F}{\delta\Phi}$ is the variation of $F$, $\mathbf{M}^{-1}\mathbf{R}$ is known as the friction operator and $\mathbf{R}^{-1}\mathbf{M}$ the mobility operator when they exist. This equation system provides relaxation dynamics for the nonequilibirum state to return to equilibria in dissipative systems or to oscillate in nondissipative systems. For simplicity, we consider the case of $\mathbf{R}=\mathbf{I}$ in this paper so that $\mathbf{M}$ is the mobility operator.

In general, the mobility operator can be decomposed into two parts:

where $\mathbf{M}_{a}$ is the antisymmetric part and $\mathbf{M}_{s}$ is the symmetric part. The time rate of change of the free energy is given by

under proper boundary conditions on $\Phi$. If $\mathbf{M}_{a}=0$ and $\mathbf{M}_{s}\geq 0$, the system is called a dissipative system since $\textrm{d}F/\textrm{d}t\leq 0$. While only $\mathbf{M}_{s}=0$, the system is a conservative system because $\textrm{d}F/\textrm{d}t=0$. Apparently, (1.3) is an equation deduced from (1.1). It can therefore be viewed as a consistent constraint for (1.1). When approximating equation (1.1) numerically, one would like to preserve energy dissipation property at the discrete level, i.e., to arrive at an approximate equation to (1.3) that is derivable from the approximate equation of (1.1). A numerical algorithm of such a property is called an energy stable algorithm or scheme.

In many thermodynamical models given in the form of (1.1), the total free energy is given by

where $(\bullet,\bullet)$ is the inner product in a $L^{2}$ space, $\mathcal{L}$ is a linear, self-adjoint, positive definite operator and $f$ or $(f,1)$ is bounded below for the physically accessible states of $\Phi$. This is also known as the gradient flow model in the literature. In this model, if we view $F$ as one of the thermodynamical variables, (1.1), (1.3), and (1.4) constitute an over-determined system of equations with $n+2$ equations and $n+1$ unknowns $(\Phi,F)$. This system of equations is consistent and solvable should gradient flow model (1.1) is solvable with free energy (1.4) and proper initial and boundary conditions. However, the structural consistency among the equations in the system can be easily broken under small perturbations to the system, leading to a system that is not solvable because of inconsistency among the equations. We label this scenario in the over-determined system of equations as structural unstable.

Traditionally, to retain consistency among the equations in the over-determined system in numerical approximations, one discretizes (1.1) and (1.4), and then try to put together a consistent discretized energy dissipation equation or inequality for (1.3) to arrive at energy stable algorithms. However, there is no guarantee that this approach would be successful. Then, there came about the convex splitting method for the gradient flow with a free energy given as a difference of two convex functions, in which one implements a mismatched temporal discretization on the two convex parts of the free energy, respectively, to achieve energy stability [6]; and the stabilizer approach [15], in which one modifies the gradient flow equation by adding higher order perturbations. All these approaches try to put together a consistent, discrete energy dissipation equation following the traditional approach mentioned above within the over-determined partial differential equation system.

Recently, the energy quadratization (EQ) approach coined by Xiaofeng Yang, Jia Zhao and Qi Wang following the work of Ganzales and Tierra on liquid crystal models, provides a systematic approach to derive energy-dissipation-rate preserving schemes for thermodynamically consistent models [11, 17, 19]. In this approach, one introduces new unknowns (known as auxiliary variables) to reformulate the energy into a quadratic form and in the meantime supplement the system with the same number of dynamical equations for the new auxiliary variables. This process is known as the EQ reformulation, which effectively embeds the gradient flow model into a higher dimensional phase space with a quadratic free energy while retaining the thermodynamically consistent structure and property (including the energy dissipation property.) The EQ reformulated gradient flow model together with the definition of the quadratic free energy and the energy dissipation equation constitute a reformulated over-determined, consistent and solvable gradient flow system. One designs linear, energy stable schemes from the EQ reformulated gradient flow model. Analogously, the SAV method follows the same idea [14].

We illustrate the idea using an example here. We assume $f$ only depends on $\Phi$, but not its spatial derivatives. In this case, one introduces a $q$ variable to transform the free energy into a quadratic form:

where $A$ is a constant large enough to make $q$ well-defined. Then, one takes the time derivative of the algebraic equation of $q$ to arrive at an equivalent dynamical equation of $q$. The EQ-reformulated system in a higher dimensional phase space recasts into

Denoting $\Psi=\left(\Phi,\;q\right)^{T}$, the above system is rewritten into [10]

where $\mathcal{N}=\mathcal{A}^{*}\mathbf{M}\mathcal{A},$ $\mathcal{A}^{*}$ is the adjoint operator of $\mathcal{A},$

$\mathcal{B}$ is a linear, self-adjoint, positive definite operator. The energy dissipation rate is given by

Recently, Gong, Zhao and Wang showed that one can devise arbitrarily high order energy stable schemes to solve the over-determined system consisting of (1.5), (1.7) and (1.8) [10, 7, 9, 8]. In contrast to the traditional method, the EQ approach guarantees the consistency and solvability in the discretized system due to the reformulated free energy is quadratic whereas the traditional approach may not. Most recently, a class of new methods rooted in the scalar auxiliary variable approach have been devised to ensure linearity as well as energy stability in the numerical schemes developed for gradient flow models [1, 18].

Here, we take look at the issue of developing energy stable numerical approximations from a new perspective. We begin with the over-determined, consistent system of equations consisting of (1.1), (1.3), and (1.4). Now that this system is structurally unstable, it can easily lose the consistency among the equations subject to any small perturbations. How can we make it structurally stable. One way to make it structurally stable is to extend the system by adding supplementary variables to make it well-determined, i.e., the number of unknowns equals to the number of equations. There are virtually unlimited number of ways this can be done. To be consistent with the original problem, we accomplish this using perturbations and require that the well-determined, perturbed system include the original system as a special case. We would like to take a new point of view towards next present a new paradigm to devise energy dissipation rate preserving schemes for the over-determined system consisting of (1.1), (1.3) and (1.4) to retain structural consistency and solvability in the over-determined system. We name it the supplementary variable method (SVM). We will show that this includes the newly proposed Lagrange multiplier method based on the SAV approach by Cheng, Liu and Shen as well as the generalized SAV method developed by Yang and Dong lately as special cases [4, 18]. This method originates from the idea on how to enforce structural stability in the over-determined equation system by augmenting it with supplementary variables. Specifically, we perturb the system consisting of (1.1), (1.3) and (1.4) by adding a perturbation:

where $\alpha(t)(g[\Phi],h(t),j(t))^{T}$ is the perturbation, $\alpha(t)$ is a function of time, $g[\Phi]$ is a prescribed functional of $\Phi$, $h$ and $j$ are prescribed functions of t. $\alpha(t)$ is called the supplementary variable. If $\alpha(t)=0$, this perturbed equation system reduces to the original, consistent and over-determined system. By adding $\alpha$, however, we lift the dynamical system into a higher dimensional phase space to remove its over-determinedness and to effectively provide it with structural stability. The gained structural stability warrants us to obtain energy-dissipation-rate preserving numerical algorithms for (1.1) under very general conditions so long as we keep the perturbation the same order as the truncation error when we discretize the system. We next detail the implementation of this method in one simple gradient flow case.

We explicitly rewrite system (1.1), (1.3) as follows

To derive energy-dissipation-rate preserving numerical approximations, we firstly approximate the energy-dissipation-rate equation of the system given in (2.1)-(2.3) as follows

where $\tau$ is the time step,

It follows from (2.4) that

This is a 2th order approximation to the energy-dissipation-rate equation with a truncation error in the order of $O(\tau^{2})$.

Secondly, we modify (2.1) by a time-dependent supplementary variable $\alpha(t)$ together with a user supplied functional $g[\Phi]$:

where $g[\Phi]$ is a given functional that may depend on $\Phi$ and its derivatives. There is a grate deal of flexibility in determining how to modify the gradient flow model with a supplementary variable. It is clearly an open problem for this approach. We simply present two special implementations in this study to illustrate the idea.

• 1.

  One is

  $$g[\Phi]=\mathbf{M}f^{\prime}(\Phi),$$

  (2.10)

  which leads to

  $$\partial_{t}\Phi=-\mathbf{M}\Big{(}\mathcal{L}\Phi+(1-\alpha)f^{\prime}(\Phi)\Big{)}.$$

  (2.11)

  Here, the chemical potential is modified into

  $$\mu=\mathcal{L}\Phi+(1-\alpha)f^{\prime}(\Phi).$$

  (2.12)

• 2.

  The other is

  $$g[\Phi]=-\mathbf{M}\Big{(}\mathcal{L}\Phi+f^{\prime}(\Phi)\Big{)},$$

  (2.13)

  which implies

  $$\partial_{t}\Phi=-(1+\alpha)\mathbf{M}\big{(}\mathcal{L}\Phi+f^{\prime}(\Phi)\big{)}.$$

  (2.14)

  Here, the mobility is modified into $(1+\alpha)\mathbf{M}$.

We note that a system consisting of (2.3), (2.11) and (2.12) is equivalent to the reformulated model proposed in [4] using what they call the Lagrange multiplier method. The supplementary variable $\alpha$ can also be viewed as a perturbation variable since when $\alpha(t)=0$, the modified/perturbed system reduces to the original one. By supplementing the over-determined system with a new variable without introducing an equation, we effectively remove the over-determinedness by making the modified system a well-determined in a higher dimensional phase space! In addition, the modified model reduces to the original one at $\alpha=0$. We reiterate that there is a great deal of flexibility to place the supplementary variable in the gradient flow model, which differentiates it from the Lagrange multiplier method and the generalized SAV method.

Thirdly, we apply the implicit-explicit Crank-Nicolson scheme in time to (2.9) to arrive at the following semi-discrete system

where $\widetilde{\Phi}^{n+\frac{1}{2}},~{}\widetilde{F}^{n+1}$ have been obtained from (2.8).

Next we discuss how to solve system (2.15) efficiently . Let

It follows from the first equation of (2.15) that

Then we plug $\Phi^{n+1}$ into the second equation in (2.15) to obtain

which is an algebraic equation for $\beta.$ In general, it can have multiple solutions, but one of them must be close to 0 and it approaches to zero as $\tau\to 0$. So, we solve for this solution using an iterative method such as the Newton iteration with 0 as the initial condition, it generally converges to a solution close to 0 when $\tau$ is not too large. After obtaining $\beta$, we update $\Phi^{n+1}$ using (2.20). Following the work of Refs. [2, 3], the existence of solution $\beta$ is guaranteed under the conditions of the following theorem.

We present a numerical example to demonstrate the practicability, accuracy, as well as energy stability of the proposed schemes. For simplicity, we use periodic boundary conditions for the numerical example below and name the resulting scheme SVM-I and SVM-II when (2.10) and (2.13) are adopted, respectively.

The numerical results demonstrate the proposed schemes based on the supplementary variable approach can predict fast coarsening dynamics of the Cahn-Hilliard model accurately and outperform SAV-CN in solution accuracy and FICN in efficiency. Here we simply present two convenient implementations of supplementary variables for developing energy-dissipation-rate preserving schemes. There can be many other ways guided by this paradigm to achieve energy stability, better computational efficiency and accuracy. We expect to see more property preserving numerical algorithms developed for thermodynamically consistent models guided by this paradigm in the future.

Research is partially supported by the Foundation of Jiangsu Key Laboratory for Numerical Simulation of Large Scale Complex Systems (202001, 202002), the Natural Science Foundation of Jiangsu Province (award BK20180413), the National Natural Science Foundation of China (award 11801269 and NSAF-U1930402) and National Science Foundation of US (award DMS-1815921 and OIA-1655740) and a GEAR award from SC EPSCoR/IDeA Program.