A Roadmap for Discretely Energy-Stable Schemes for Dissipative Systems Based on a Generalized Auxiliary Variable with Guaranteed Positivity

TL;DR

This paper introduces a unified framework for creating discretely energy-stable schemes for dissipative systems using a generalized auxiliary variable that guarantees positivity and stability regardless of time step size.

Contribution

The paper proposes the gPAV method, a novel approach that ensures positivity and energy stability in numerical schemes for dissipative systems, applicable to a wide class of problems.

Findings

Guaranteed positivity of the auxiliary variable at discrete level.

Energy stability of the proposed schemes for various dissipative systems.

Effective performance and robustness demonstrated through numerical experiments.

Abstract

We present a framework for devising discretely energy-stable schemes for general dissipative systems based on a generalized auxiliary variable. The auxiliary variable, a scalar number, can be defined in terms of the energy functional by a general class of functions, not limited to the square root function adopted in previous approaches. The current method has another remarkable property: the computed values for the generalized auxiliary variable are guaranteed to be positive on the discrete level, regardless of the time step sizes or the external forces. This property of guaranteed positivity is not available in previous approaches. A unified procedure for treating the dissipative governing equations and the generalized auxiliary variable on the discrete level has been presented. The discrete energy stability of the proposed numerical scheme and the positivity of the computed auxiliary variable have been proved for general dissipative systems. The current method, termed gPAV (generalized Positive Auxiliary Variable), requires only the solution of linear algebraic equations within a time step. With appropriate choice of the operator in the algorithm, the resultant linear algebraic systems upon discretization involve only constant and time-independent coefficient matrices, which only need to be computed once and can be pre-computed. Several specific dissipative systems are studied in relative detail using the gPAV framework. Ample numerical experiments are presented to demonstrate the performance of the method, and the robustness of the scheme at large time step sizes.

Tables

Table 1. Table 1: Simulation parameter values for convergence tests of Cahn-Hilliard equation.

parameter	value	parameter	value
$C_0$	$1$	$\lambda$	$0.01$
$m_0$	$0.01$	$\eta$	$0.1$
$t_0$	$0.1$	$t_f$	$0.2$ (spatial tests) or $1.1$ (temporal tests)
Element order	(varied)	Elements	$2$
$\Delta t$	(varied)	$\Delta t_{\min }$	$1e-4$
S	$1$ (variable mobility solver)	$S$	$\sqrt{\frac{4{\gamma }_0\lambda }{m_0\Delta t}}$ or $\sqrt{\frac{4{\gamma }_0\lambda }{m_0\Delta t_{\min }}}$ (constant mobility solver)
$ℱ\left(R\right)$	$R$	${\varphi }_0$	${\varphi }_{in}$
$m_c\left({\varphi }_0\right)$	$m\left({\varphi }_0\right)$

Equations

$$2f(x)f^{\prime }(x)=1.$$

$$\frac{\partial \mathbf{u}}{\partial t}=\mathbf{F}(\mathbf{u})+\mathbf{f}(\mathbf{x},t)$$

$$\mathbf{B}(\mathbf{u})={\mathbf{f}}_b,\text{on}\Gamma$$

$$\mathbf{u}(\mathbf{x},t=0)={\mathbf{u}}_{in}(\mathbf{x})$$

$$E_{tot}(t)=E_{tot}[\mathbf{u}]={\int }_{\Omega }e(\mathbf{u})d\Omega ,$$

$$\frac{d{E}_{tot}}{dt}={\int }_{\Omega }\frac{\partial e}{\partial \mathbf{u}}\cdot \frac{\partial \mathbf{u}}{\partial t}d\Omega ={\int }_{\Omega }\frac{\partial e}{\partial \mathbf{u}}\cdot \left[\mathbf{F}(\mathbf{u})+\mathbf{f}\right]d\Omega ,$$

$${\int }_{\Omega }\frac{\partial e}{\partial \mathbf{u}}\cdot [\mathbf{F}(\mathbf{u})+\mathbf{f}]d\Omega =-{\int }_{\Omega }V(\mathbf{u})d\Omega +{\int }_{\Omega }{V}_s(\mathbf{f},\mathbf{u})d\Omega +{\int }_{\Gamma }{B}_s({\mathbf{f}}_b,\mathbf{u})d\Gamma ,$$

$$V_s(\mathbf{f},\mathbf{u})=0,\text{if}\mathbf{f}=0.$$

$$\frac{d{E}_{tot}}{dt}=-{\int }_{\Omega }V(\mathbf{u})d\Omega +{\int }_{\Omega }{V}_s(\mathbf{f},\mathbf{u})d\Omega +{\int }_{\Gamma }{B}_s({\mathbf{f}}_b,\mathbf{u})d\Gamma .$$

$$B_s({\mathbf{f}}_b,\mathbf{u})=0\text{if}{\mathbf{f}}_b=0,\text{on}\Gamma .$$

$$V(\mathbf{u})⩾0.$$

$$E(t)=E[\mathbf{u}]={\int }_{\Omega }e(\mathbf{u})d\Omega +C_0,$$

$$\{\begin{array}{rl}& \mathcal{F}(\chi )>0,\text{for}\chi >0;\\ & \mathcal{G}(\chi )>0,\text{for}\chi >0.\end{array}$$

$$R(t)=\mathcal{G}(E),$$

$$E(t)=\mathcal{F}(R),$$

$${\mathcal{F}}^{\prime }(R)\frac{dR}{dt}={\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \frac{\partial \mathbf{u}}{\partial t}d\Omega$$

$$\mathcal{F}(\chi )={\chi }^m,\mathcal{G}(\chi )={\chi }^{1/m},m\in {\mathcal{Z}}^{+}=\{1,2,3,...\};$$

$$\mathscr{F}(\chi)=\frac{e_{0}}{2}\ln\Big{(}\frac{\kappa_{0}+\chi}{\kappa_{0}-\chi}\Big{)},\quad\mathscr{G}(\chi)=\kappa_{0}\tanh\left(\frac{\chi}{e_{0}}\right),$$

$$\frac{\partial\bm{u}}{\partial t}=\bm{F}_{L}(\bm{u})+\frac{\mathscr{F}(R)}{E}\Big{(}\bm{F}(\bm{u})-\bm{F}_{L}(\bm{u})\Big{)}+\bm{f},$$

$$\begin{array}{rl}{\mathcal{F}}^{\prime }(R)\frac{dR}{dt}=& {\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \frac{\partial \mathbf{u}}{\partial t}d\Omega +\left[\frac{\mathcal{F}(R)}{E}-1\right]{\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \left[{\mathbf{F}}_L(\mathbf{u})+\mathbf{f}\right]d\Omega \\ & +\frac{\mathcal{F}(R)}{E}\left({\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \left[\mathbf{F}(\mathbf{u})-{\mathbf{F}}_L(\mathbf{u})\right]d\Omega -{\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \left[\mathbf{F}(\mathbf{u})-{\mathbf{F}}_L(\mathbf{u})\right]d\Omega \right)\\ & +\left[1-\frac{\mathcal{F}(R)}{E}\right]\left|{\int }_{\Omega }{V}_s(\mathbf{f},\mathbf{u})d\Omega +{\int }_{\Gamma }{B}_s({\mathbf{f}}_b,\mathbf{u})d\Gamma \right|\\ =& {\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \frac{\partial \mathbf{u}}{\partial t}d\Omega -{\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \left({\mathbf{F}}_L(\mathbf{u})+\frac{\mathcal{F}(R)}{E}\left[\mathbf{F}(\mathbf{u})-{\mathbf{F}}_L(\mathbf{u})\right]+\mathbf{f}\right)d\Omega \\ & +\frac{\mathcal{F}(R)}{E}{\int }_{\Omega }\frac{\partial e}{\partial \mathbf{u}}\cdot \left[\mathbf{F}(\mathbf{u})+\mathbf{f}\right]d\Omega +\left[1-\frac{\mathcal{F}(R)}{E}\right]\left|{\int }_{\Omega }{V}_s(\mathbf{f},\mathbf{u})d\Omega +{\int }_{\Gamma }{B}_s({\mathbf{f}}_b,\mathbf{u})d\Gamma \right|\end{array}$$

$$\begin{array}{rl}{\mathcal{F}}^{\prime }(R)\frac{dR}{dt}=& {\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \frac{\partial \mathbf{u}}{\partial t}d\Omega -{\int }_{\Omega }{e}^{\prime }(\mathbf{u})\cdot \left({\mathbf{F}}_L(\mathbf{u})+\frac{\mathcal{F}(R)}{E}\left[\mathbf{F}(\mathbf{u})-{\mathbf{F}}_L(\mathbf{u})\right]+\mathbf{f}\right)d\Omega \\ & +\frac{\mathcal{F}(R)}{E}\left[-{\int }_{\Omega }V(\mathbf{u})d\Omega +{\int }_{\Omega }{V}_s(\mathbf{f},\mathbf{u})d\Omega +{\int }_{\Gamma }{B}_s({\mathbf{f}}_b,\mathbf{u})d\Gamma \right]\\ & +\left[1-\frac{\mathcal{F}(R)}{E}\right]\left|{\int }_{\Omega }{V}_s(\mathbf{f},\mathbf{u})d\Omega +{\int }_{\Gamma }{B}_s({\mathbf{f}}_b,\mathbf{u})d\Gamma \right|.\end{array}$$

$$R(0)=\mathcal{G}(E(0)),\text{where}E(0)={\int }_{\Omega }e({\mathbf{u}}_{in})d\Omega +C_0.$$

$${\chi }^{n+\frac{3}{2}}=\frac{3}{2}{\chi }^{n+1}-\frac{1}{2}{\chi }^n,{\chi }^{n+\frac{1}{2}}=\frac{3}{2}{\chi }^n-\frac{1}{2}{\chi }^{n-1},$$

$$\displaystyle\frac{\partial\chi}{\partial t}\Big{|}^{n+1}=\frac{\chi^{n+\frac{3}{2}}-\chi^{n+\frac{1}{2}}}{{\Delta}t}=\frac{1}{\Delta{t}}\Big{(}\frac{3}{2}\chi^{n+1}-2\chi^{n}+\frac{1}{2}\chi^{n-1}\Big{)},$$

$${\bar{\chi }}^{n+1}=2{\chi }^n-{\chi }^{n-1},$$

$$D_{\mathscr{F}}(\chi)\big{|}^{n+1}=\frac{\mathscr{F}(\chi^{n+\frac{3}{2}})-\mathscr{F}(\chi^{n+\frac{1}{2}})-\mathscr{F}^{\prime}(\chi^{n+1})\cdot(\chi^{n+\frac{3}{2}}-\chi^{n+\frac{1}{2}})}{\|\chi^{n+\frac{3}{2}}-\chi^{n+\frac{1}{2}}\|^{2}}(\chi^{n+\frac{3}{2}}-\chi^{n+\frac{1}{2}})+\mathscr{F}^{\prime}(\chi^{n+1}),$$

$$D_{\mathscr{F}}(\chi)\Big{|}^{n+1}\cdot\left(\frac{3}{2}\chi^{n+1}-2\chi^{n}+\frac{1}{2}\chi^{n-1}\right)=D_{\mathscr{F}}(\chi)\Big{|}^{n+1}\cdot\left(\chi^{n+\frac{3}{2}}-\chi^{n+\frac{1}{2}}\right)=\mathscr{F}(\chi^{n+\frac{3}{2}})-\mathscr{F}(\chi^{n+\frac{1}{2}}).$$

$$D_{\mathscr{F}}(\chi)\big{|}^{n+1}=\frac{\mathscr{F}(\chi^{n+\frac{3}{2}})-\mathscr{F}(\chi^{n+\frac{1}{2}})}{\chi^{n+\frac{3}{2}}-\chi^{n+\frac{1}{2}}}=\frac{\mathscr{F}(\chi^{n+\frac{3}{2}})-\mathscr{F}(\chi^{n+\frac{1}{2}})}{\frac{3}{2}\chi^{n+1}-2\chi^{n}+\frac{1}{2}\chi^{n-1}},$$

$$\displaystyle\frac{\partial\bm{u}}{\partial t}\Big{|}^{n+1}=\bm{F}_{L}(\bm{u}^{n+1})+\xi\Big{[}\bm{F}(\bar{\bm{u}}^{n+1})-\bm{F}_{L}(\bar{\bm{u}}^{n+1})\Big{]}+\bm{f}^{n+1},$$

$$\xi =\frac{\mathcal{F}(R^{n+3/2})}{E[{\tilde{\mathbf{u}}}^{n+3/2}]},$$

$$E[{\tilde{\mathbf{u}}}^{n+3/2}]={\int }_{\Omega }e({\tilde{\mathbf{u}}}^{n+3/2})d\Omega +C_0,$$

$$\mathbf{B}({\mathbf{u}}^{n+1})={\mathbf{f}}_b^{n+1},\text{on}\Gamma ,$$

$$\begin{split}D_{\mathscr{F}}(R)\big{|}^{n+1}&\left.\frac{dR}{dt}\right|^{n+1}=\int_{\Omega}e^{\prime}({\bm{u}}^{n+1})\cdot\left.\frac{\partial\bm{u}}{\partial t}\right|^{n+1}d\Omega\\ &-\int_{\Omega}e^{\prime}(\bm{u}^{n+1})\cdot\left(\bm{F}_{L}(\bm{u}^{n+1})+\xi\Big{[}\bm{F}(\bar{\bm{u}}^{n+1})-\bm{F}_{L}(\bar{\bm{u}}^{n+1})\Big{]}+\bm{f}^{n+1}\right)d\Omega\\ &+\xi\left[-\int_{\Omega}V(\tilde{\bm{u}}^{n+1})d\Omega+\int_{\Omega}V_{s}(\bm{f}^{n+1},\tilde{\bm{u}}^{n+1})d\Omega+\int_{\Gamma}B_{s}(\bm{f}_{b}^{n+1},\tilde{\bm{u}}^{n+1})d\Gamma\right]\\ &+(1-\xi)\left|\int_{\Omega}V_{s}(\bm{f}^{n+1},\tilde{\bm{u}}^{n+1})d\Omega+\int_{\Gamma}B_{s}(\bm{f}_{b}^{n+1},\tilde{\bm{u}}^{n+1})d\Gamma\right|.\end{split}$$

$$\frac{\mathcal{F}(R^{n+\theta })}{E[{\tilde{\mathbf{u}}}^{n+\theta }]}=1+\mathcal{O}(\Delta t{)}^2,$$

Full text

Keywords: *energy stability; unconditional stability; dissipative systems; conservative systems; auxiliary variables; positivity *

1 Introduction

Dissipative systems are of immense interest to science and engineering. Physical systems encountered in the real world are dissipative, thanks to the second law of thermodynamics. In dissipative systems there exists a storage function that is bounded from below Willems1972 . We will refer to this function as the energy in the current work. Dissipative systems are distinguished from general dynamical systems by the dissipation inequality, which basically states that the increase in storage of the system over a time interval cannot exceed the supply to the system during that interval Willems1972 ; Willems2007 . The governing partial differential equations (PDE) describing dissipative systems are typically nonlinear, and they satisfy a balance equation for the energy (or entropy) as an embodiment of the dissipation inequality GrootM1984 ; Ottinger2005 ; AndersonMW1998 ; LowengrubT1998 ; AbelsGG2012 ; Dong2018 .

A highly desirable property for numerical algorithms for dissipative systems is the preservation of the energy dissipation (or conservation) on the discrete level. This not only preserves one important aspect of the underlying structure of the continuous system HairerLW2006 , but more practically also provides a control on the numerical stability in actual computer simulations. The history for such strategies is long and they can be traced to at least the work of CourantFL1928 on discrete energy conservation for finite difference approximations in the 1920s. While energy-stable schemes for specific domains of science and engineering have been under intensive studies and these efforts have borne invaluable fruits, the schemes and methods developed usually have only limited applicability across domains. The energy-stable schemes for one area are hardly transferable to a different field, and they can hardly shed light on the development of such types of schemes in new unexplored domains. Unified techniques that can be broadly applied to treat different PDEs from different domains for devising energy-stable schemes are generally lacking. The metaphor used in Iserles2008 (page 139) to compare the motley collection of PDEs to a hugh unhappy family (each unhappy in its own way; Tolstoy, “Anna Karenina”) seems fitting in describing this situation (see also Celledonietal2012 ).

Occasionally, certain methods appear and seem to be broadly applicable to a wide class of problems spanning different areas. The average vector field (AVF) method Celledonietal2012 ; QuispelM2008 and the discrete variational derivative method (DVDM) FurihataM2011 , both of which can be traced to the idea of discrete gradients Gonzalez1996 ; McLachlanQR1999 , are two such examples. For gradient systems that can be expressed into the form $\frac{\partial\bm{u}}{\partial t}=\bm{L}\cdot\frac{\delta H}{\delta\bm{u}},$ where $\bm{L}$ is an anti-symmetric or negative semi-definite matrix, $\bm{u}$ is the field variable, $H(\bm{u})$ is the energy functional and $\frac{\delta H}{\delta\bm{u}}$ denotes the variational derivative, the AVF and DVDM methods can preserve the energy conservation (resp. energy dissipation) discretely. We refer the reader to e.g. Furihata1999 ; DahlbyO2011 ; MiyatakeM2014 ; CaiLW2018 ; EidnesOR2018 (among others) for related and variants of these methods. A potential drawback of these methods is their computational cost. Because these are fully implicit schemes and the governing PDEs are in general nonlinear, these methods will entail the solution of nonlinear algebraic equations on the discrete level. Consequently, some nonlinear algebraic solver (e.g. Newton type methods) will be required for computing the field functions, and the associated computational cost can be substantial.

In the current work we present a framework for devising energy-stable schemes for general dissipative systems that can potentially be useful and applicable to different domains. Our method does not require the governing PDEs to be in any particular form, as long as they are dissipative (or conserving). When devising the energy-stable numerical schemes, we are particularly mindful of the computational cost involved therein. The resultant energy-stable schemes from our method involve only the solution of linear algebraic equations when computing the field functions within a time step, and no nonlinear algebraic solver is needed. Furthermore, with appropriate choice of the operator in the scheme, the resultant linear algebraic systems upon discretization can involve only constant and time-independent coefficient matrices, which only need to be computed once and can be pre-computed during pre-processing. Thanks to these properties, the presented method and the resultant energy-stable schemes are computationally very competitive and attractive. In terms of the computational cost the presented method enjoys a notable advantage when compared with the aforementioned methods.

The key to achieving the above useful properties for general dissipative systems in the presented method lies in the introduction of a generalized auxiliary variable. The generalized auxiliary variable introduced here is inspired by the scalar auxiliary variable (SAV) approach proposed by ShenXY2018 , and to a lesser extent, by the invariant energy quadratization (IEQ) method Yang2016 , both of which are devised for gradient flows; see also e.g. ShenX2018 ; GongZYW2018 ; ChengS2018 ; Zhaoetal2018 ; KouSW2018 ; LiZW2019 ; YangLD2019 ; Yang2019 (among others) for extensions and applications of these techniques. In SAV a scalar-valued auxiliary variable is defined, as the square root of the shifted potential energy integral. In IEQ an auxiliary field variable is defined, as the square root of the shifted potential energy density function. With these auxiliary variables, energy-stable schemes can be devised for gradient flows and their discrete energy stability can be proven in the SAV and IEQ methods. In both SAV and IEQ, the use of the square root function is critical to the proof of the discrete energy stability of the resultant numerical schemes, due to the interesting property that the square root is the only function form that satisfies the relation

[TABLE]

In the current work we will show that the square root function is not essential to devising energy-stable schemes. In the generalized auxiliary variable method developed here, the auxiliary variable (a scalar number) can be defined by a rather general class of functions (conditions specifically given in Section 2.1) in terms of the energy functional, which is why the method is termed “generalized”, and the resultant numerical schemes can be proven to be discretely energy stable.

The method presented here is applicable to general dissipative systems, which is another key difference from previous auxiliary-variable approaches. The ability to deal with general dissipative systems hinges on how the governing PDEs are treated based on the generalized auxiliary variable and how the generalized auxiliary variable is numerically treated on the discrete level. A unified procedure for treating discretely the dissipative governing equations and the generalized auxiliary variable has been presented. These numerical treatments have drawn inspirations from the recent developments in LinYD2019 ; YangD2018 for incompressible Navier-Stokes equations and for the incompressible two-phase flows, which are not gradient-type systems.

The generalized auxiliary variable method proposed herein has another remarkable property: The computed values for the auxiliary variable are guaranteed to be positive on the discrete level. Such a property is not available in the SAV (or IEQ) method. In both SAV and IEQ, as well as in the current method, the auxiliary variable is computed discretely by solving an associated dynamic equation, which is derived based on the definition of the auxiliary variable in terms of the square root function in SAV and IEQ or a general function in the current method. The auxiliary variable physically should be positive according to its definition. However, this positivity property is in general not guaranteed in the computed values for the auxiliary variable, because they are obtained by numerically solving a differential equation. Indeed, in numerical experiments we have observed negative values for the computed auxiliary variable using the previous methods, especially at large time step sizes. With the current method, on the other hand, we can prove that the computed values for the generalized auxiliary variable are guaranteed to be positive, regardless of the time step sizes or the external forces. The guaranteed positivity of the auxiliary variable in the current method is intimately related to and is critical to the proof of discrete energy stability of the proposed numerical schemes.

Because of these crucial properties, we will refer to the framework proposed herein as “gPAV”, which stands for the generalized Positive Auxiliary Variable method.

In this paper we consider general dissipative systems and outline the gPAV procedure for devising discretely energy-stable schemes. The discrete energy stability of the proposed numerical scheme and the positivity property of the computed auxiliary variable will be proven for general dissipative systems. As already mentioned, the gPAV method requires only the solution of linear algebraic equations within a time step, and with appropriate choice of the operator in the algorithm, the resultant linear algebraic systems involve only constant and time-independent coefficient matrices that can be pre-computed. We demonstrate the gPAV procedure by looking into three specific dissipative systems: a chemo-repulsion model Gonzalez2019 , the Cahn-Hilliard equation CahnH1958 with constant and variable mobility, and the nonlinear Klein-Gordon equation Strauss1978 . Ample numerical experiments are provided for each system to demonstrate the performance of the algorithm and the effects of the parameters.

The current work contains several new aspects: (i) the framework for developing discretely energy-stable schemes for general dissipative systems; (ii) the generalized auxiliary variable introduced herein; and (iii) the guaranteed positivity of the computed auxiliary variable on the discrete level. Some other aspects, such as the generalization of the numerical algorithm as discussed in Remarks 2.5 and 2.6, are also potentially useful to other researchers and the community.

The remainder of this paper is structured as follows. In Section 2 we introduce a generalized auxiliary variable and present the gPAV framework for devising discretely energy-stable schemes for general dissipative systems. The discrete energy stability of the presented algorithm and the positivity of the computed auxiliary variable will be proven. The solution algorithm for implementing the proposed energy-stable scheme will be presented. An alternative formulation for the energy-stable scheme will also be discussed in this section. Then in the three subsequent sections (Sections 3–5) we apply the gPAV framework to three specific dissipative systems (a chemo-repulsion model, Cahn-Hilliard equation with constant and variable mobility, and Klein-Gordon equation). Ample numerical experiments are provided to demonstrate the performance of the method for each system, and numerical results with large time step sizes are presented to show the robustness of the proposed scheme. Section 6 concludes the discussions with some closing remarks. In Appendix A we provide a method for approximating the variables for the first time step, which guarantees the positivity of the computed auxiliary variable to start off. This startup procedure is important for the proof of discrete energy stability of the presented numerical scheme.

2 The gPAV Framework for Energy-Stable Schemes for Dissipative Systems

Consider a domain $\Omega$ in two or three dimensions and a dissipative system on this domain, whose dynamics is described by,

[TABLE]

where $\bm{x}$ and $t$ denote the spatial coordinate and time, $\bm{u}(\bm{x},t)$ denotes the state variables of the system and can be a scalar- or vector-valued field function, and $\bm{f}(\bm{x},t)$ is an external source term (hereafter referred to as the external force). $\bm{F}(\bm{u})$ is an operator that gives rise to the dissipative dynamics of the system and can be nonlinear in general. Equation (2.1) is supplemented by the boundary condition

[TABLE]

where $\Gamma$ denotes the domain boundary, $\bm{f}_{b}$ is an external source term on the boundary, which will be referred to as the external boundary force hereafter, and $\bm{B}$ is assumed to be a linear operator for the sake of simplicity. The initial condition is

[TABLE]

where $\bm{u}_{in}(\bm{x})$ is the initial distribution of the state variable.

Because the system is dissipative, there exists a storage function that is bounded from below Willems1972 , which hereafter will be referred to as the energy,

[TABLE]

where $e(\bm{u})$ is the energy density function. The evolution of the energy is described by

[TABLE]

where we have used equation (2.1). With integration by part, the right-hand-side (RHS) of equation (2.5) can be transformed into

[TABLE]

where $V_{s}(\bm{f},\bm{u})=\frac{\partial e}{\partial\bm{u}}\cdot\mathbf{f}$ denotes the volume terms involving the external force $\bm{f}$, which satisfies the property

[TABLE]

The rest of the volume terms are denoted by $-V(\bm{u})$, not involving $\bm{f}$. $B_{s}(\bm{f}_{b},\bm{u})$ denotes the boundary terms, which may involve the boundary source term ($\bm{f}_{b}$) through the boundary conditions.

Substituting equation (2.6) into equation (2.5), we arrive at the following energy balance equation for the system,

[TABLE]

We assume that the boundary conditions (2.2) satisfy the following property,

[TABLE]

The dissipative nature of the system ensures that $\frac{dE_{tot}}{dt}\leqslant 0$ in the absence of the external forces (i.e. $\bm{f}=0$ and $\bm{f}_{b}=0$). Because the domain $\Omega$ can be arbitrary, it follows that $V(\bm{u})$ must be non-negative, i.e.

[TABLE]

2.1 Reformulated Equivalent System

To facilitate energy-stable numerical approximations of the system (2.1), we define a shifted energy of the following form

[TABLE]

where $C_{0}$ is a chosen energy constant such that $E(t)>0$ for $0\leqslant t\leqslant T$, and $T$ is the time interval on which the computation is to be carried out. Note that for a physical system the energy is bounded from below, and thus $C_{0}$ can always be found.

Let $\mathscr{F}$ denote a one-to-one increasing differentiable function, with its inverse $\mathscr{F}^{-1}=\mathscr{G}$, satisfying the property

[TABLE]

We define a scalar variable $R(t)$ by

[TABLE]

where $E(t)$ is the shifted energy given by (2.11). $R(t)$ then satisfies the following evolution equation,

[TABLE]

which is obtained by taking the time derivative of equation (2.13b) and using equation (2.11).

Remark 1**.**

The choice for $\mathscr{F}$ and $\mathscr{G}$ is rather general. Some examples are,

[TABLE]

or

[TABLE]

where $\kappa_{0}$ and $e_{0}$ are positive constants. It is important to notice that a function like $\mathscr{F}(\chi)=\chi^{2m+1}$ (with an integer $m\geqslant 0$) or $\mathscr{F}(\chi)=\ln(1+\chi)$ does not automatically guarantee that $\mathscr{F}(\chi)>0$ with arbitrary $\chi$. However, if one can ensure that the argument satisfies $\chi>0$, the property $\mathscr{F}(\chi)>0$ can be guaranteed with such choices of functions when defining $R(t)$. This point is critical in the subsequent development of the numerical algorithm.

Noting that $\frac{\mathscr{F}(R)}{E}=1,$ we rewrite equation (2.1) into an equivalent form

[TABLE]

where $\bm{F}_{L}(\bm{u})$ is a chosen linear operator about $\bm{u}$. $\bm{F}_{L}(\bm{u})$ should be of the same spatial order as $\bm{F}(\bm{u})$. For improved accuracy $\bm{F}_{L}(\bm{u})$ should be an approximation of $\bm{F}(\bm{u})$ in some way, such as the linear component of $\bm{F}(\bm{u})$ or a linearized approximation of $\bm{F}(\bm{u})$. For improved numerical efficiency $\bm{F}_{L}(\bm{u})$ should be easy to compute and implement.

Remark 2.1**.**

$\bm{F}(\bm{u})$* often consists of linear components and nonlinear components for many systems, and oftentimes one can choose the linear components as the $\bm{F}_{L}$ operator. One can also add/subtract certain linear operators, and treat one part freely and the other part together with $\frac{\mathscr{F}(R)}{E}$ as in equation (2.17). By choosing an $\bm{F}_{L}$ operator that involves only time-independent (or constant) coefficients, the resultant method will become computationally very efficient, because the coefficient matrices for the linear algebraic systems upon discretization will be time-independent and therefore can be pre-computed when solving the field variables. This point will become clearer from later discussions.*

We reformulate equation (2.14) as follows,

[TABLE]

where it can be noted that a number of zero terms have been incorporated. In the above equation \big{|}(\cdot)\big{|} denotes the absolute value of $(\cdot)$. In light of (2.6), we transform equation (2.18) into the final reformulated equivalent form

[TABLE]

The reformulated system consists of equations (2.17) and (2.19), the boundary conditions (2.2), the initial condition (2.3) for $\bm{u}$, and the following initial condition for $R(t)$,

[TABLE]

In the reformulated system, the dynamic variables are $\bm{u}$ and $R(t)$, which are coupled in the equations (2.17) and (2.19). $E(t)$ is given by equation (2.11). Note that in this system $R(t)$ is determined by solving the coupled system of equations, not by using the equation (2.13a).

{comment}

Remark 2.2**.**

In the developments of this section, we will assume periodic boundary conditions for $\bm{u}$ (or essentially no boundaries for the domain). This is because with no knowledge about the specific structure (e.g. spatial order) of the operator $\bm{F}(u)$, it is difficult to discuss the boundary conditions in general terms. The goal of this section is to demonstrate the general strategy for the development of energy-stable schemes for these systems. When looking into the applications of the general framework to specific systems in subsequent sections, some commonly-used boundary conditions will be considered. Numerical treatment of those boundary conditions will be illustrated in these later sections.

2.2 An Energy-Stable Scheme

We next present an energy-stable scheme for the reformulated system consisting of (2.17) and (2.19), together with the boundary condition (2.2) and the initial conditions (2.3) and (2.20).

Let $n\geqslant 0$ denote the time step index, and $(\cdot)^{n}$ represent the variable $(\cdot)$ at time step $n$, corresponding to the time $t=n\Delta t$, where $\Delta t$ is the time step size. If a real-valued parameter $\theta$ is involved, $(\cdot)^{n+\theta}$ represents the variable $(\cdot)$ at time step ($n+\theta$), corresponding to the time $(n+\theta)\Delta t$.

Let $\chi$ denote a generic scalar or vector-valued variable. We consider the following second-order approximations:

[TABLE]

where (2.21b) is the second-order backward differentiation formula (BDF) and $\bar{\chi}^{n+1}$ is an explicit approximation of $\chi^{n+1}$. We also consider the following second-order approximation of \left.\frac{d\mathscr{F}(\chi)}{d\chi}\right|^{n+1}=\mathscr{F}^{\prime}(\chi)\Big{|}^{n+1} based on the discrete directional derivative Gonzalez1996 ,

[TABLE]

which satisfies the property

[TABLE]

Note that in these equations $\chi^{n+3/2}$ and $\chi^{n+1/2}$ are given by (2.21a). If $\chi$ represents a scalar-valued variable, one can also approximate \mathscr{F}^{\prime}(\chi)\Big{|}^{n+1} by

[TABLE]

which satisfies the same property (2.23).

We propose the following scheme to approximate the reformulated system:

[TABLE]

In the above equations, $\left.\frac{\partial{\bm{u}}}{\partial t}\right|^{n+1}$ and $\left.\frac{dR}{dt}\right|^{n+1}$ are defined by (2.21b), $\left.D_{\mathscr{F}}(R)\right|^{n+1}$ is defined by (2.22) (or (2.24)), $\bar{\bm{u}}^{n+1}$ is defined by (2.21c), and $R^{n+3/2}$ is defined by (2.21a). $\tilde{\bm{u}}^{n+1}$ and $\tilde{\bm{u}}^{n+3/2}$ are second-order approximations of $\bm{u}^{n+1}$ and $\bm{u}^{n+3/2}$, respectively, to be specifically defined later in (2.43).

Remark 2.3**.**

It is critical to note that in the scheme (2.25a)–(2.25e), $\frac{\mathscr{F}(R)}{E[\bm{u}]}$ is approximated at step ($n+\frac{3}{2}$) while the other variables are approximated at step ($n+1$). This feature, together with the approximation (2.22), allows $R^{n+1}$ to be computed from a linear algebraic equation (no nonlinear algebraic solver), and endows the scheme with the property that the computed $R^{n+1}$ and $\mathscr{F}(R^{n+1})$ (resp. $R^{n+3/2}$ and $\mathscr{F}(R^{n+3/2})$, for all $n\geqslant 0$) are guaranteed to be positive. These points will become clear from later discussions. It should be noted that the approximation $\frac{\mathscr{F}(R^{n+3/2})}{E[\tilde{\bm{u}}^{n+3/2}]}$ at step ($n+3/2$) is a second-order approximation of $\frac{\mathscr{F}(R)}{E}=1$. In fact, the approximation involving any real parameter $\theta$,

[TABLE]

is a second-order approximation of $\frac{\mathscr{F}(R)}{E}=1$, as long as $R^{n+\theta}$ and $\tilde{\bm{u}}^{n+\theta}$ are second-order approximations of $R(t)$ and $\bm{u}(t)$ at time $(n+\theta)\Delta t$. Therefore, the approximation in (2.25b) does not affect the second-order accuracy of the scheme.

The scheme given by (2.25a)–(2.25e) has the following property.

Theorem 2.1**.**

In the absence of the external force and external boundary force (i.e. $\bm{f}=\bm{0}$ and $\bm{f}_{b}=0$), the following relation holds with the scheme (2.25):

[TABLE]

if the approximation of $R(t)$ at time step $\frac{1}{2}$ is positive, i.e. Y_{0}=R^{n+1/2}\Big{|}_{n=0}>0.

Proof.

By equations (2.21b) and (2.22), we have

[TABLE]

Taking the $L^{2}$ inner product between equation (2.25a) and $e^{\prime}({\bm{u}}^{n+1})$, and adding the resultant equation to equation (2.25e) and noting equation (2.28), we arrive at

[TABLE]

where we have used equation (2.25b), and $S_{0}$ is defined by

[TABLE]

Then it follows that, if $\bm{f}=0$ and $\bm{f}_{b}=0$,

[TABLE]

where we have used the relations (2.7) and (2.9).

Note that $E[\tilde{\bm{u}}^{n+3/2}]>0$ and $V(\tilde{\bm{u}}^{n+1})\geqslant 0$, in light of (2.11) and (2.10). If $Y_{0}=R^{n+1/2}|_{n=0}>0$, then $\mathscr{F}(Y_{0})>0$ based on the property (2.12). By induction, we can conclude from equation (2.31) that $\mathscr{F}(R^{n+3/2})>0$ for all $n\geqslant 0$. The inequality in (2.27) then holds. We therefore conclude that, if $R^{n+1/2}|_{n=0}>0$,

[TABLE]

Thus, the scheme is unconditionally energy stable with respect to the modified energy $\mathcal{F}(R)$, if the approximation of $R(t)$ at time step $\frac{1}{2}$ is positive. ∎

There are many ways to approximate $R(t)$ to ensure that it is positive at time step $\frac{1}{2}$ and that the overall scheme is second-order accurate in time. One such method is given in the Appendix A. Therefore we have the following result:

Theorem 2.2**.**

With $\bm{u}^{1}$ and $R^{1}$ approximated using the method from Appendix A, in the absence of external forces ($\bm{f}=0$ and $\bm{f}_{b}=0$), the scheme represented by (2.25a)–(2.25e) is unconditionally energy-stable in the sense of the relation (2.32).

Remark 2.4**.**

If the functional form of $\mathscr{F}(\chi)$ is such that $\mathscr{F}(\chi)\geqslant 0$ for all $\chi\in(-\infty,\infty)$, e.g. $\mathscr{F}(\chi)=\chi^{2m}$ (with an integer $m\geqslant 1$), then the scheme given by (2.25a)–(2.25e) is unconditionally energy stable regardless of the approximation of $R(t)$ at the time step $\frac{1}{2}$.

Remark 2.5**.**

The scheme (2.25) is devised by enforcing the system of equations consisting of (2.17), (2.19) and (2.2) at time step ($n+1$), approximating $\frac{\mathscr{F}(R)}{E}$ at time step ($n+\frac{3}{2}$), and employing the approximations (2.21a)–(2.22). Inspired by the recent work YangLD2019 , we can generalize this scheme by enforcing the system of equations at time step ($n+\theta$), where $\theta$ is a real-valued parameter, to arrive at a family of energy-stable schemes.

In brief, let us consider the following second-order approximations at time step ($n+\theta$) with $\theta\geqslant\frac{1}{2}$: ($\chi$ denoting a generic variable, and $\beta\geqslant 0$ denoting a real parameter below)

[TABLE]

and the following approximation of $\left.\frac{d\mathscr{F}(\chi)}{d\chi}\right|^{n+\theta}=\left.\mathscr{F}^{\prime}(\chi)\right|^{n+\theta}$ based on discrete directional derivative,

[TABLE]

These approximations satisfy the following properties:

[TABLE]

Note that the parameter $\beta\geqslant 0$ in (2.33b) can often be used to control the numerical dissipation of the approximations, which will be useful for approximating energy-conserving systems. An example will be given with the Klein-Gordon equation in a later section. The scheme given in (2.25) corresponds to $\theta=1$ and $\beta=\frac{1}{4}$.

By approximating the terms in equations (2.17), (2.19) and (2.2) at time step ($n+\theta$), except for the term $\frac{\mathscr{F}(R)}{E}$, which will be approximated at time step ($n+\theta+\frac{1}{2}$), and employing the approximations (2.33a)–(2.34), one can prove that the resultant family of schemes (with $\theta$ and $\beta$ as parameters) is unconditionally energy-stable. The details will not be provided here.

{comment}

Note that in Yang2019 , a family of energy stable scheme for Cahn-Hilliard type equations has been proposed and the key idea is to approximate the variable $\chi$ at time step $(n+\theta)$ as follows

[TABLE]

However, these approximations are not sufficient to obtain an energy stable scheme at time step $(n+\theta)$ for the current method. In equations (2.21b) and (2.22), the approximation at the time step $(n+\frac{3}{2})$ and $(n+\frac{1}{2})$ (adjacent to the time step $(n+1)$) is used such that \mathscr{F}^{\prime}[R]\Big{|}^{n+1}\frac{dR}{dt}\Big{|}^{n+1} leads to the difference of $\mathscr{F}[R]$ at the adjacent time steps, see equation (2.28). In view of this, we adopt the corresponding second-order approximations at time step $(n+\theta)$ as follows:

where $\theta\in[0.5,1.5]$ and $\beta$ in (2.33b) is a constant to adjust the artificial dissipation of the approximation. Therefore, by using the above approximations at time step $(n+\theta)$ for equations (2.25a)-(2.25e) and approximate equation (2.25b) at time step $(n+\theta+\frac{1}{2}),$ it is direct to show that the resultant scheme is unconditionally energy stable.

2.3 Solution Algorithm

Let us now consider how to implement the algorithm represented by equations (2.25a)-(2.25e). We first introduce some notations ($\chi$ again denoting a generic variable):

[TABLE]

Then the approximation in (2.21b) can be written as

[TABLE]

Inserting notation (2.38) into equation (2.25a), we have

[TABLE]

Note that $\bar{\bm{u}}^{n+1}$ and $\hat{\bm{u}}$ are both explicitly known, and $\xi$ is an unknown depending on $\bm{u}^{n+1}$. Taking advantage of the fact that $\xi$ is a scalar number instead of a field function and the linearity of the operator $\bm{B}$ in the boundary condition (2.2), we introduce two field functions $(\bm{u}_{1}^{n+1},\bm{u}_{2}^{n+1})$ as solutions to the following two linear systems:

[TABLE]

[TABLE]

Since the operator $\bm{F}_{L}$ is chosen to be a linear operator and relatively easy to compute, $\bm{u}_{1}^{n+1}$ and $\bm{u}_{2}^{n+1}$ can be solved efficiently from these equations. Then we have the following result.

Theorem 2.3**.**

Given scalar value $\xi,$ the following function solves the system consisting of equations (2.25a) and (2.25d):

[TABLE]

where $\bm{u}_{1}^{n+1}$ and $\bm{u}_{2}^{n+1}$ are given by the equations (2.40a)-(2.41b).

The scalar value $\xi$ still needs to be determined. Define

[TABLE]

which are second-order approximations of $\bm{u}^{n+1}$ and $\bm{u}^{n+3/2}$. These field variables can be explicitly computed after ${\bm{u}}_{1}^{n+1}$ and $\bm{u}_{2}^{n+1}$ are obtained. By equation (2.25b), we have

[TABLE]

Note that equation (2.25e) can be transformed into equation (2.29). Inserting equation (2.44) into equation (2.29) leads to the solution for $\xi$,

[TABLE]

where $\tilde{\bm{u}}^{n+1}$ and $\tilde{\bm{u}}^{n+3/2}$ are given by (2.43), $S_{0}$ is given by equation (2.30), and $E[\tilde{\bm{u}}^{n+3/2}]$ is computed by equation (2.25c).

In light of equations (2.44) and (2.21a), we can then compute $R^{n+1}$ by

[TABLE]

The following result holds.

Theorem 2.4**.**

The scalar value $\xi$ computed by equation (2.45) and the variable $R^{n+1}$ ($n\geqslant 0$) computed by equation (2.46) are always positive, if the approximation of $R(t)$ at time step $\frac{1}{2}$ is positive, i.e. $Y_{0}=R^{n+1/2}|_{n=0}>0$.

Proof.

If $Y_{0}=R^{n+1/2}|_{n=0}>0$, then $\mathscr{F}(Y_{0})>0$ based on (2.12). Since $E(\bm{u})$ is a positive function, $V(\bm{u})\geqslant 0$ and $|S_{0}|-S_{0}\geqslant 0$, we conclude by induction $\xi$ computed from (2.45) is always positive.

Note that $R^{0}=R(0)>0$ according to equation (2.20). In light of the property (2.12), we conclude that $R^{n+3/2}$ and $R^{n+1}$ computed from equation (2.46) are both positive. ∎

Using the method from the Appendix A can ensure the positiveness of the approximation of $R(t)$ at the time step $\frac{1}{2}$. We have the following result.

Theorem 2.5**.**

With $\bm{u}^{1}$ and $R^{1}$ computed based on the method from Appendix A, the $\xi$ given by (2.45) and $R^{n+1}$ and $R^{n+3/2}$ given by (2.46) satisfy the property

[TABLE]

for all $n\geqslant 0$, regardless of the external forces $\bm{f}$ and $\bm{f}_{b}$ and the time step size $\Delta t$.

Combining the above discussions, we arrive at the solution procedure for solving the system consisting of equations (2.25a)-(2.25e). Given $(\bm{u}^{n},R^{n})$, we compute $(\bm{u}^{n+1},R^{n+1})$ through the following steps:

Solve equations (2.40a)–(2.40b) for $\bm{u}_{1}^{n+1}$;

Solve equations (2.41a)–(2.41b) for $\bm{u}_{2}^{n+1}$. 2. 2.

Compute $\tilde{\bm{u}}^{n+1}$ and $\tilde{\bm{u}}^{n+{3}/{2}}$ based on equation (2.43);

Compute $E[\tilde{\bm{u}}^{n+\frac{3}{2}}]$, $\int_{\Omega}V(\tilde{\bm{u}}^{n+1})$ and $S_{0}$ based on equations (2.11), (2.6) and (2.30). 3. 3.

Compute $\xi$ based on equation (2.45). 4. 4.

Compute $\bm{u}^{n+1}$ based on equation (2.42). Compute $R^{n+1}$ based on equation (2.46).

It can be noted that the numerical scheme and the solution algorithm developed in this section has several attractive properties: (i) Only linear systems need to be solved for the field variables $\bm{u}$ within a time step. Moreover, with appropriate choice for the $\bm{F}_{L}$ operator, the system can involve only constant and time-independent coefficient matrices, which can be pre-computed. Therefore, the solution for $\bm{u}$ will be computationally very efficient. (ii) The auxiliary variables $R$ and $\xi$ can be computed by a well-defined explicit formula, and no nonlinear algebraic solver is involved. Their computed values are guaranteed to be positive. (iii) The auxiliary variable $R$ can be defined by a rather general class of functions ($\mathscr{F}$ and $\mathscr{G}$) using the method developed here. (iv) The scheme is unconditionally energy-stable for general dissipative systems.

2.4 An Alternative Formulation and Energy-Stable Scheme

The numerical formulation presented in the previous subsections is not the only way to devise energy-stable schemes for dissipative systems. In this subsection we outline an alternative formulation and associated energy-stable scheme. The process is analogous to the developments in the sections 2.1–2.3. So many details will be omitted in the following discussions.

The main idea with the alternative formulation is to realize that $\frac{R(t)}{\mathscr{G}(E)}=1$ with the auxiliary variable $R(t)$ defined in (2.13a). Therefore, one can potentially employ $\frac{R}{\mathscr{G}(E)}$, instead of $\frac{\mathscr{F}(R)}{E}$, in the numerical formulations. With appropriate reformulation and treatments of different terms, it turns out that a discretely energy-stable scheme can be obtained with similar attractive properties, such as the guaranteed positiveness of the computed values for the variable $R(t)$.

Note that $R(t)$ is defined by (2.13a), where $\mathscr{G}$ is a one-to-one increasing differentiable function with $\mathscr{G}(\chi)>0$ and $\mathscr{G}^{\prime}(\chi)>0$ for $\chi>0$. $R(t)$ satisfies the following dynamic equation

[TABLE]

where $E(t)$ is defined by (2.11).

We reformulate equation (2.1) into

[TABLE]

where the notations follow those defined in previous subsections. Analogously, by incorporating appropriate zero terms we can transform (2.48) into

[TABLE]

The reformulated system now consists of equations (2.49) and (2.50), the boundary condition (2.2), and the initial conditions (2.3) and (2.20).

We discretize the reformulated system as follows:

[TABLE]

In these equations $\bar{\bm{u}}^{n+1}$ is defined by (2.21c), $R^{n+3/2}$ and $R^{n+1/2}$ are defined by (2.21a), and $\tilde{\bm{u}}^{n+1}$ and $\tilde{\bm{u}}^{n+3/2}$ are second-order approximations of $\bm{u}^{n+1}$ and $\bm{u}^{n+3/2}$ respectively to be specified later.

Taking the $L^{2}$ inner product between $\mathscr{G}^{\prime}(E[\tilde{u}^{n+3/2}])e^{\prime}(\bm{u}^{n+1})$ and equation (2.51a), and summing up the resultant equation and equation (2.51e), we get

[TABLE]

where $S_{0}$ is given by the equation (2.30). In the absence of external forces ($\bm{f}=0$ and $\bm{f}_{b}=0$), $S_{0}=0$ and equation (2.52) leads to

[TABLE]

where we have used (2.51b). Note that $E[\tilde{\bm{u}}^{n+3/2}]>0$, $V(\tilde{\bm{u}}^{n+1})\geqslant 0$, and that $\mathscr{G}(\chi)>0$ and $\mathscr{G}^{\prime}(\chi)>0$ for $\chi>0$. By induction we can conclude from (2.53) that $R^{n+3/2}\geqslant 0$ (for all $n\geqslant 0$) if the approximation of $R(t)$ at time step $\frac{1}{2}$ is non-negative. Equation (2.52) then leads to the following result.

Theorem 2.6**.**

In the absence of external forces ($\bm{f}=0$ and $\bm{f}_{b}=0$), if the approximation of $R(t)$ at time step $\frac{1}{2}$ is non-negative, the scheme given by (2.51a)–(2.51e) is unconditionally energy-stable in the sense that

[TABLE]

In the Appendix A, we have presented a method for computing the first time step, which can ensure that the approximation of $R(t)$ at step $\frac{1}{2}$ is positive. This leads to the following result.

Theorem 2.7**.**

In the absence of external forces ($\bm{f}=0$ and $\bm{f}_{b}=0$), when the first time step is approximated using the method from Appendix A, the numerical scheme given by (2.51a)–(2.51e) is unconditionally energy-stable in the sense of equation (2.54).

The scheme represented by (2.51a)–(2.51e) can be implemented in a similar way to that of Section 2.3, with the following steps:

• •

Compute $\bm{u}_{1}^{n+1}$ and $\bm{u}_{2}^{n+1}$ by solving equations (2.40a)–(2.41b).

• •

Define $\tilde{\bm{u}}^{n+1}$ and $\tilde{\bm{u}}^{n+3/2}$ again by equations (2.43). These variables can be computed.

• •

Compute $\xi$ based on equation (2.52), specifically by

[TABLE]

where $S_{0}$ is given by (2.30).

• •

Compute ${\bm{u}}^{n+1}$ by equation (2.42). Compute $R^{n+1}$ by

[TABLE]

where we have used equations (2.51b) and (2.21a).

Noting the positiveness of energy $E(t)$ and the other functions involved in equations (2.55) and (2.56), we have the following result.

Theorem 2.8**.**

If the first time step is approximated using the method from Appendix A, regardless of the external forces $\bm{f}$ and $\bm{f}_{b}$ and the time step size $\Delta t$, the computed values for $\xi$ and $R^{n+1}$ with the scheme (2.51a)–(2.51e) satisfy the property,

[TABLE]

for all time steps.

Remark 2.6**.**

In the current paper we have used the total energy (shifted) $E_{tot}(t)$ (see equation (2.11)) to define the auxiliary variable $R(t)$. One can also define an auxiliary variable based on a part of the total energy. Suppose the total energy of the system can be written as

[TABLE]

where each of the energy components $E_{1}[\bm{u}]$ and $E_{2}[\bm{u}]$ is bounded from below. One can define an auxiliary variable $R(t)$ based on e.g. $E_{2}(t)$ (shifted appropriately),

[TABLE]

where the chosen energy constant $C_{0}$ is to ensure that $E_{s}(t)>0$. By appropriate reformulation of the system one can devise energy-stable schemes in an analogous way. We refer the reader to YangD2018 for such an energy-stable scheme for incompressible two-phase flows with different densities and viscosities for the two fluids, which corresponds to a specific mapping function $\mathscr{F}(R)=R^{2}$. A drawback with this lies in that one needs to solve a nonlinear algebraic equation (or a quadratic equation), albeit about a scalar number, when computing the auxiliary variable, and that the property for guaranteed positiveness of the computed auxiliary-variable values will be lost.

In the subsequent sections, we consider three dissipative (or conserving) systems (a chemotaxis model, Cahn-Hilliard equation, and Klein-Gordon equation) as specific applications and demonstrations of the gPAV method developed in this section.

3 A Chemo-Repulsion Model

3.1 Model and Numerical Scheme

Consider the following repulsive-productive chemotaxis model with a quadratic production term (see e.g. Gonzalez2019 ) in a domain $\Omega$ (with boundary $\Gamma$):

[TABLE]

where $p(u)=u^{2}$ is the quadratic production term, $u(\bm{x},t)\geq 0$ is the cell density, and $v(\bm{x},t)\geq 0$ is the chemical concentration. $f_{1}$, $f_{2}$, $d_{a}$ and $d_{b}$ denote the volume and boundary source terms, respectively. $u_{in}$ and $v_{in}$ are the initial distributions of the field variables. This system is dissipative in the absence of the source terms, with the total energy given by (see Gonzalez2019 )

[TABLE]

By taking the $L^{2}$ inner products between (3.1a) and $u,$ and between (3.1b) and $-\dfrac{1}{2}\nabla^{2}v$, summing them up and performing integration by part and imposing boundary conditions in (3.1c), we can obtain the following energy balance equation:

[TABLE]

Following the gPAV procedure from section 2, we define a shifted energy according to equation (2.11)

[TABLE]

where $C_{0}$ is a chosen energy constant such that $E(t)>0$. Define a scalar auxiliary variable $R(t)$ according to equation (2.13a). Thus, equation (2.14) becomes

[TABLE]

Following equations (2.17)-(2.19), we reformulate equations (3.1a)-(3.1b) into the following equivalent form:

[TABLE]

By incorporating the following zero terms into the right hand side of equation (3.5),

[TABLE]

we can transform this equation into

[TABLE]

where we have used the fact $\frac{\mathscr{F}(R)}{E}=1$ and the boundary conditions (3.1c).

The reformulated equivalent system consist of equations (3.6a)-(3.7) and (3.1c)-(3.1d). The energy-stable scheme for this system is as follows:

[TABLE]

and

[TABLE]

In these equations, $\left.\frac{\partial u}{\partial t}\right|^{n+1}$, $\left.\frac{\partial v}{\partial t}\right|^{n+1}$ and $\left.\frac{dR}{dt}\right|^{n+1}$ are defined by equation (2.21b). $\bar{u}^{n+1}$ and $\bar{v}^{n+1}$ are defined by (2.21c). $\tilde{u}^{n+1}$ and $\tilde{v}^{n+1}$ are second-order approximations of $u^{n+1}$ and $v^{n+1}$ to be specified later in (3.21). $\tilde{u}^{n+3/2}$ and $\tilde{v}^{n+3/2}$ are second-order approximations of $u^{n+3/2}$ and $v^{n+3/2}$ to be specified later in (3.22). $S_{0}$ in equation (3.9) is given by

[TABLE]

where

[TABLE]

These equations are supplemented by the following initial conditions

[TABLE]

Theorem 3.1**.**

In the absence of the external force $f_{1}=f_{2}=0,$ and with homogeneous boundary conditions $d_{a}=d_{b}=0,$ the scheme consisting of (3.8a)-(3.9) is unconditionally energy stable in the sense that:

[TABLE]

if the approximation of $R(t)$ at the time step $\frac{1}{2}$ is non-negative.

This theorem can be proved in a way analogous to Theorem 2.1. We can apply the method from Appendix A to this chemo-repulsion model for the first time step, and this ensures that $R^{n+1/2}|_{n=0}>0$.

3.2 Solution Algorithm and Implementation

Using the notation (2.38), we rewrite equations (3.8a)-(3.8b) into

[TABLE]

Barring the unknown scalar $\xi,$ (3.14) and (3.15) are two decoupled Helmholtz-type equations about $u^{n+1}$ and $v^{n+1},$ respectively.

Note that $\xi$ is a scalar number instead of a field function, we define two sets of variables $(u_{i}^{n+1},v_{i}^{n+1})$ $(i=1,2)$ as the solutions to the following equations:

[TABLE]

Then we have the following result: Given the scalar number $\xi,$ the following field functions solve the system consisting of equations (3.14)-(3.15):

[TABLE]

where $(u_{i}^{n+1},v_{i}^{n+1})$ $i=1,2$ is given by equations (3.16)-(3.19), respectively.

Once $(u_{i}^{n+1},v_{i}^{n+1})$ $i=1,2$ are known, we determine $\tilde{u}^{n+1}$, $\tilde{v}^{n+1}$, $\tilde{u}^{n+3/2}$ and $\tilde{v}^{n+3/2}$ according to (2.43), specifically by

[TABLE]

In light of equations (3.1b), (3.21) and (3.11) , we compute $\nabla^{2}\tilde{v}^{n+1}$ in equation (3.9) by

[TABLE]

where $\left.\frac{\partial v}{\partial t}\right|^{*,n+1}$ is given by (3.11).

Combining equations (3.8a)–(3.8b) and (3.9), and using the property (2.23), we have

[TABLE]

This gives rise to

[TABLE]

in which $S_{0}$ is given by (3.10), $\nabla^{2}\tilde{v}$ is to be computed by (3.23), and $E[\tilde{u}^{n+3/2},\tilde{v}^{n+3/2}]$ is given by (3.8d). With $\xi$ known, $R^{n+1}$ and $(u^{n+1},v^{n+1})$ can be evaluated directly by (2.46) and (3.20), respectively.

We employ $C^{0}$-continuous high-order spectral elements for spatial discretizations in our implementation. Note that equations (3.16)–(3.19) involve Helmholtz type equations with Neumann type boundary conditions. The weak formulations of these equations are: Find $u_{i}^{n+1}$ and $v_{i}^{n+1}$ $\in H^{1}(\Omega)$ for $i=1,2,$ such that

[TABLE]

for $\forall\varphi\in H^{1}(\Omega)$, where

[TABLE]

These weak forms can be discretized using $C^{0}$ spectral elements in the standard way KarniadakisS2005 .

3.3 Numerical Results

3.3.1 Convergence Rate

We first employ a manufactured analytical solution to the chemo-repulsion model to demonstrate the spatial and temporal convergence rates of the proposed algorithm.

Consider the computational domain $\Omega=[0,1]^{2}$ and the following contrived solution to the system (3.1) on this domain

[TABLE]

The external forces $f_{1}(\bm{x},t),$ $f_{2}(\bm{x},t)$ and boundary forces $d_{a}(\bm{x},t),$ $d_{b}(\bm{x},t)$ therein are chosen such that the expressions in (3.27) satisfy (3.1) .

The domain is discretized with four equal-sized quadrilateral elements. The initial cell density $u_{in}$ and initial chemical concentration $v_{in}$ are given according to the analytic expressions in (3.27) by setting $t=0.$ We simulate this problem from $t=0$ to $t=t_{f}.$ Then we compare the numerical solutions of $u$ and $v$ at $t=t_{f}$ with the analytic solutions in (3.27) and various norms of the errors are computed. The element order and time step sizes are varied systematically in order to investigate their effects on the numerical errors. We employ the function $\mathscr{F}(R)=R$ for defining the auxiliary variable $R(t)$ and the energy constant $C_{0}=1$ in the following convergence tests.

We first study the spatial convergence rate. A fixed $t_{f}=0.1$ and $\Delta{t}=0.001$ is employed and the element order is varied systematically between 2 and 20. We record the errors at $t=t_{f}$ between the numerical solution and the contrived solution (3.27) in both $L^{\infty}$ and $L^{2}$ norms with respect to the element orders. Figure 3.1(a) shows these numerical errors as a function of the element order. We observe an exponential decrease of the numerical errors with increasing element order, and a level-off of the error curves beyond element order 10 and 8, respectively for $u$ and $v$, due to the saturation of temporal errors.

The study of the temporal convergence rate is summarized by the results in Figure 3.1(b). Here we fix the integration time $t_{f}=1.0$ and the element order at a large value 18, and vary $\Delta{t}$ systematically between $0.2$ and $1.953125\times 10^{-4}.$ This figure demonstrates the $L^{\infty}$ and $L^{2}$ errors of $u$ and $v$ as a function of $\Delta{t}$. It is evident that the proposed scheme has a second-order convergence rate in time.

3.3.2 Study of Unconditional Stability and Effect of Algorithmic Parameters

We next consider the test problem used in Gonzalez2019 , and show the efficiency and unconditional stability of the method proposed here. Consider the domain $\Omega=[0,2]^{2}$ and the initial distributions for the cell density $u$ and chemical concentration $v$ in this domain given by

[TABLE]

The external forces and boundary forces in (3.1) are set to $f_{1}=f_{2}=d_{a}=d_{b}=0.$ The computational domain is discretized with 400 equal-sized quadrilateral elements, and the element order is fixed to be 10.

Figures 3.2 and 3.3 demonstrate the dynamics of the system. These results are obtained with $\Delta{t}=10^{-5}$, $\mathscr{F}(R)=R$ and $C_{0}=1$ in the numerical algorithm. Figure 3.2 shows the evolution of the cell density $u(\bm{x},t)$ with a temporal sequence of snapshots of the distribution visualized by the contour plots. The $z$ coordinate corresponds to $u$ in these plots. The system exhibits a very rapid dynamics. The initial cell density has a Gaussian type distribution, taking a minimal value 0.0001 at the domain center $\bm{x}_{0}=(1,1)$ and gradually approaching the maximal value 10.0001 near the domain boundary. In a very short time $t=10^{-2},$ the maximal density increases to around 16, attained near the boundary of a circular region with radius 0.6 and center at $\bm{x}_{0}$; see Figure 3.2(b). Then the maximal density gradually moves from the circular boundary to the domain boundary between $t=2\times 10^{-2}$ and $t=7.5\times 10^{-2}$; see Figure 3.2(c)-(f). The high density near the domain boundary then appears to diffuse to the low density region near the center $\bm{x}_{0}$, and the system finally reaches an equilibrium state between $t=0.1$ and $t=0.5$ with a constant density level; see Figure 3.2(g)-(i). Figure 3.3 illustrates the evolution of the chemical concentration $v(\bm{x},t)$. Figure 3.3(a) shows the distribution of the initial chemical concentration. It has also a Gaussian type distribution, with a maximal value 100.0001 at the origin $\bm{x}_{0}$ and decreasing to 0.0001 gradually near the domain boundary. The concentration diffuses rapidly between $t=0$ to $t=5\times 10^{-2}$ (Figures 3.3(a)-(e)), and the maximal concentration decreases to around 10 at the origin. From $t=7.5\times 10^{-2}$ to $t=0.2,$ the contrast in the concentration levels in the domain becomes even smaller (Figure 3.3(f)-(h)), and the concentration reaches its equilibrium with a constant level around 36.6 (Figure 3.3(i)).

Figure 3.4 shows time histories of three quantities: $E(t)$, $\mathscr{F}(R)$, and $\xi=\frac{\mathscr{F}(R)}{E(t)}$, corresponding to three time step sizes $\Delta t=10^{-5}$, $10^{-4}$ and $10^{-3}$. Note that $E(t)$ is computed based on equation (3.4), $\mathscr{F}(R)$ is computed based on the $R(t)$ obtained from the algorithm, and $\xi$ is computed based on equation (3.25). These results are obtained with $\mathscr{F}(R)=R$ and $C_{0}=1$ in the algorithm. It is observed from Figure 3.4(a) that both $E(t)$ and $\mathscr{F}(R)$ decrease over time and gradually level off at certain levels over time. A comparison of the $E(t)$ histories obtained using different $\Delta{t}$ indicates that they are quite close, with only some slight difference on the interval between $t=0.002$ and $t=0.15$. Note that $\mathscr{F}(R)$ is an approximation of $E(t)$ in the current method, and the evolution equation for $R(t)$ stems from this relation; see equations (2.13a)–(2.14). Therefore, the difference between $E(t)$ and $\mathscr{F}(R)$, and also the quantity $\xi=\frac{\mathscr{F}(R)}{E(t)}$, can serve as an indicator of the accuracy of the simulations. If the difference between $E(t)$ and $\mathscr{F}(R)$ is small, or the deviation of $\xi$ from the unit value is small, then the simulation tends to be more accurate. On the other hand, when the difference between $E(t)$ and $\mathscr{F}(R)$ is pronounced, or the deviation between $\xi$ and the unit value is significant, it implies that $\mathscr{F}(R)$ is no longer an accurate approximation of $E(t)$ and the simulation will contain large numerical errors. Here it can be observed that $E(t)$ and $\mathscr{F}(R)$ computed with $\Delta{t}=10^{-5}$ essentially overlap with each other, indicating $\mathscr{F}(R)$ approximates well the quantity $E(t).$ However, the time histories for $E(t)$ and $\mathscr{F}(R)$ obtained with $\Delta{t}=10^{-4}$ and $10^{-3}$ exhibit noticeable discrepancies. This suggests that in these cases $\mathscr{F}(R)$ is no longer an accurate approximation of $E(t).$ We also observe from Figure 3.4(b) that $\xi$ computed by $\Delta{t}=10^{-5}$ is essentially 1, while with larger values $\Delta{t}=10^{-4}$ and $\Delta{t}=10^{-3}$ the computed $\xi$ attains values significantly smaller than 1. These results indicate that with the larger time step sizes $\Delta t=10^{-4}$ and $10^{-3}$ the simulation results contain pronounced errors and they are not accurate any more. Because this problem exhibits very rapid dynamics (see Figures 3.2 and 3.3), to capture such dynamics accurately the requirement on $\Delta t$ is very stringent.

Thanks to its energy-stable nature, our algorithm can produce stable simulation results even with very large $\Delta{t}$ values. This is demonstrated by Figure 3.5 with several large time step sizes, ranging from ${\Delta}t=0.01$ to $\Delta t=10$, with $\mathscr{F}(R)=R$ and $C_{0}=1$ in the algorithm. We show the time histories of the total energy $E_{tot}(t)$ (see equation (3.2)) and the ratio $\xi=\frac{\mathscr{F}(R)}{E}$ for a much longer simulation (up to $t=1000$). The long time histories demonstrate that the computations with these large $\Delta{t}$ values are indeed stable using the current algorithm. On the other hand, because these $\Delta t$ values are very large, we cannot expect that the results will be accurate. This is evident from the values of $\xi$ in Figure 3.5(b). These time histories for $\xi$ tend to level off at very small but positive values, with large deviations from the unit value. It is noted that the simulations are nonetheless stable, regardless of $\Delta t$.

When defining the modified energy $E(t)$ (see equation (3.4)) we have incorporated an energy constant $C_{0}$. The goal of $C_{0}$ is to ensure that $E(t)>0$ for all time, even in certain extreme cases such as when $E_{tot}=0$, so that $\frac{1}{E(t)}$ (as in $\frac{\mathscr{F}(R)}{E(t)}$) is always well-defined. We observe that the choice of the $C_{0}$ value seems to have some influence on the numerical results. This effect is illustrated by Figure 3.6. Here we employ $\mathscr{F}(R)=R$ and $\Delta{t}=10^{-5}$ and $10^{-4}$, and depict the time histories of $E_{tot}(t)$ and $\xi$ obtained with several $C_{0}$ values ($C_{0}=1$, $10^{3}$, $10^{6}$ and $10^{10}$). With the smaller $\Delta t=10^{-5}$, the obtained $E_{tot}$ histories corresponding to different $C_{0}$ values overlap with one another. The computed $\xi$ values are essentially $1$, with a discrepancy on the order of magnitude of $10^{-6}$. This discrepancy between the computed $\xi$ and the unit value is associated with the smaller $C_{0}=1$ and $10^{3}$. With the larger $C_{0}=10^{6}$ and $10^{10}$, no difference can be observed at this scale. This suggests that with a small $\Delta t$ (so that the simulation result is generally accurate) a larger $C_{0}$ value tends to give rise to more accurate $\xi$ in terms of its discrepancy from the unit value. Figures 3.6(c) and (d) are the corresponding result obtained with a larger $\Delta t=10^{-4}$, in which case the simulation result is no longer accurate. In this case it is observed that with the larger $C_{0}=10^{6}$ and $10^{10},$ the energy $E_{tot}$ history curves exhibit a bump, apparently artificial; see Figure 3.6(c). In contrast, with the smaller $C_{0}=1$ and $10^{3}$, such a bump is not quite obvious from the energy history curves. In addition, with the larger $C_{0}=10^{6}$ and $10^{10}$, the computed $\xi$ attains a very small value (close to 0), while $\xi$ attains a value around $0.2$ with the smaller $C_{0}=1$ and $10^{3}$. This indicates that, with larger $\Delta t$ (when simulation loses accuracy), the simulation results obtained with a smaller $C_{0}$ may be better than those obtained with a larger $C_{0}$, even though all the results become inaccurate. The results of this group of tests suggest the following. With small $\Delta t$ values, a larger $C_{0}$ tends to give rise to more accurate results in the sense that the computed $\xi$ tends to be closer to the unit value. However, a $C_{0}$ that is very large seems to have an adverse effect when $\Delta t$ becomes large, because it can lead to computed $\xi$ values that deviate from the unit value more severely. The majority of simulations in this section are performed using $C_{0}=1$.

The method developed in the current work can employ a general function $\mathscr{F}(R)$ (with inverse $\mathscr{G}$) to define the auxiliary variable $R(t)$, as long as $\mathscr{F}$ is a one-to-one increasing differentiable function satisfying (2.12). We observe that the choice for the specific mapping $\mathscr{F}$ seems to have very little or no influence on the simulation results using the current method. This point is demonstrated by Figure 3.7. Here we have considered several functions, $\mathscr{F}(R)=R^{m}$ ($m=1,2,3,4,6$) and $\mathscr{F}(R)=\frac{e_{0}}{2}\ln(\frac{\kappa_{0}+R}{\kappa_{0}-R})$ with $e_{0}=8040$ and $\kappa_{0}=10^{3}$. Figure 3.7 shows the time histories of $E_{tot}(t)$ and $\xi$ obtained using these mappings, together with a fixed $C_{0}=1$ and two time step sizes $\Delta{t}=10^{-4}$ and $10^{-5}$. It can be observed that the time history curves for both $E_{tot}(t)$ and $\xi$ corresponding to different $\mathscr{F}$ functions overlap with one another, suggesting no or very little difference in the simulation results. In particular, Figure 3.7(d) shows the $\xi$ history curves corresponding to different $\mathscr{F}$ obtained with the smaller $\Delta t$, with the vertical axis $\xi$ magnified around the unit value. It can be observed that the difference between various curves is on the order of magnitude $10^{-6}$. Since little difference in the numerical results is observed with different mapping functions $\mathscr{F}(R)$ using the current method, the majority of numerical tests reported in this and subsequent sections will be carried out using the simplest mapping $\mathscr{F}(R)=R$.

In Section 2.4 we have discussed another unconditionally energy-stable scheme (referred to as “alternative method”), which is based on an alternative formulation with $\xi=\frac{R}{\mathscr{G}(E)}$. The dynamic equation for the auxiliary variable $R(t)$ is accordingly replaced by equation (2.48). Figure 3.8 is a comparison of the time histories for $E_{tot}(t)$ and $\xi$ obtained using these two methods. The results in Figure 3.8(a) and (b) are obtained with a mapping function $\mathscr{F}(R)=R^{2}$ (or equivalently $\mathscr{G}(E)=\sqrt{E}$), and those in (c) and (d) correspond to $\mathscr{F}(R)=R^{3}$ (or $\mathscr{G}(E)=\sqrt[3]{E}$). We observe that there seems to be little difference in the computed total energy $E_{tot}(t)$. But some difference can be noted with the $\xi$ histories. The computed $\xi$ values using the current method (with $\frac{\mathscr{F}(R)}{E}$) seem to be consistently larger than those using the alternative method (with $\frac{R}{\mathscr{G}(E)}$). While all these values deviate from the unit value substantially because of the time step size $\Delta t=10^{-4}$, the deviation with the current method appears noticeably smaller than that with the alternative method. This seems to suggest that, while the simulation results using these methods are not very much different, the formulation using $\frac{\mathscr{F}(R)}{E}$ may be somewhat better than the alternative formulation using $\frac{R}{\mathscr{G}(E)}$.

4 Cahn-Hilliard Equation with Constant and Variable Mobility

We apply the gPAV method to simulate the Cahn-Hilliard equation CahnH1958 in this section. This equation has widespread applications in the phase-field modeling of materials science, two-phase and multiphase flows (see e.g. LowengrubT1998 ; Chen2002 ; LiuS2003 ; YueFLS2004 ; KimL2005 ; DingSS2007 ; DongS2012 ; Dong2012 ; Dong2014 ; LiuSY2015 ; WuX2017 ; XuLWB2019 , among others). Consider the Cahn-Hilliard equation on a domain $\Omega$ (with boundary $\Gamma$):

[TABLE]

supplemented by the initial condition

[TABLE]

In these equations, $\phi(\bm{x},t)\in[-1,1]$ is the phase field function, $f(\bm{x},t)$, $d_{a}(\bm{x},t)$ and $d_{b}(\bm{x},t)$ are prescribed source terms for the purpose of convergence testing only, and will be set to $f(\bm{x},t)=d_{a}(\bm{x},t)=d_{b}(\bm{x},t)=0$ in actual simulations. $E_{tot}$ is the free energy functional,

[TABLE]

in which $\eta$ is the characteristic interfacial thickness scale, and $\lambda$ is referred to as the mixing energy density coefficient and is related to other physical parameters. For example, for two-phase flow problems $\lambda$ is given by $\lambda=\frac{3}{2\sqrt{2}}\sigma\eta,$ where $\sigma$ is the surface tension. $\mu$ is referred to as the chemical potential, and the nonlinear term $h(\phi)$ is given by $h(\phi)=H^{\prime}(\phi)$. $H(\phi)$ is referred to as the potential free energy density function, which can take many different forms. In this paper we only consider the double-well form as given in (4.3). $m\geqslant 0$ is the mobility, and in this work we consider two cases: (i) $m=m_{0}$, and (ii) $m=m(\phi)=\max(m_{0}(1-\phi^{2}),0)$, with $m_{0}$ being a given positive constant.

We take the $L^{2}$ inner product between (4.1a) and $\mu$, perform integration by part and impose the boundary condition (4.1d). This leads to the energy balance equation,

[TABLE]

Based on equations (2.11) and (4.4), we define the shifted total energy by

[TABLE]

where $C_{0}$ is chosen to ensure $E(t)>0$. Let us define $\mathscr{F}$ and $\mathscr{G}$ and $R(t)$ based on equations (2.13a)–(2.13b). Following equation (2.14) and using (4.5), we have

[TABLE]

where the boundary condition (4.1d) has been used.

4.1 Constant Mobility

Assume that $m(\phi)=m_{0}>0$ is a constant. We reformulate equations (4.1a)–(4.1c) as follows,

[TABLE]

where $S$ is chosen constant satisfying a condition to be specified later. Note that a zero term $S(\phi-\phi)$ is added in these equations. We reformulate equation (4.6) as follows,

[TABLE]

where $\mu$ is given by (4.1b), and the following zero terms have been incorporated into the RHS,

[TABLE]

The energy-stable scheme for the equations (4.7a)–(4.7b), (4.1d) and (4.8) is as follows:

[TABLE]

and

[TABLE]

These are supplemented by the initial conditions

[TABLE]

In the above equations, $\left.\frac{\partial\phi}{\partial t}\right|^{n+1}$ and $\left.\frac{dR}{dt}\right|^{n+1}$ are defined by (2.21b), and ${\bar{\phi}}^{n+1}$ is defined by (2.21c). $\tilde{\phi}^{n+1}$, $\tilde{\phi}^{n+3/2}$ and $\tilde{\mu}^{n+1}$ are second-order approximations of $\phi^{n+1}$, $\phi^{n+3/2}$ and $\mu^{n+1}$, respectively, to be specified later in (4.25)–(4.27). $\left.\frac{\partial\phi}{\partial t}\right|^{*,n+1}$ is an approximation of $\left.\frac{\partial\phi}{\partial t}\right|^{n+1}$ to be specified later in (4.26).

Theorem 4.1**.**

In the absence of the external force $f=0,$ and with zero boundary conditions $d_{a}=d_{b}=0,$ the scheme consisting of (4.10)-(4.11) is unconditionally energy stable in the sense that

[TABLE]

if the approximation of $R(t)$ at time step $\frac{1}{2}$ is positive.

Proof.

Multiplying $-\lambda\nabla^{2}\phi^{n+1}+h(\phi^{n+1})$ to equation (4.10a), integrating over the domain, and adding the resultant equation to equation (4.11), we obtain the energy balance relation as follows:

[TABLE]

where we have used the relation (2.28). If $f=0$ and $d_{a}=d_{b}=0$, then

[TABLE]

If $R^{n+1/2}|_{n=0}>0$, one can conclude by induction that $\xi>0$ for any $n\geqslant 0$. This leads to (4.13). ∎

The method from the Appendix A can be employed to compute the first time step, which can ensure that the approximation of $R(t)$ at the step $\frac{1}{2}$ is positive.

To implement the scheme we note that equation (4.10a) can be transformed into

[TABLE]

where we have used the notation in equation (2.38). This equation can be reformulated into the following two Helmholtz type equations that are de-coupled from each other (barring the unknown scalar number $\xi$), (see e.g. DongS2012 ; YangLD2019 for details)

[TABLE]

where $\psi^{n+1}$ is an auxiliary field variable defined by (4.17b), and the constant $\alpha$ is given by and the chosen constant $S$ must satisfy

[TABLE]

In light of (4.17b) and (4.10c), the boundary condition (4.10b) can be transformed into

[TABLE]

To solve equations (4.17a)-(4.17b) together with the boundary conditions (4.19) and (4.10c), we take advantage of the fact that $\xi$ is a scalar number and introduce two sets of field functions $(\psi_{i}^{n+1},\phi_{i}^{n+1})$ $(i=1,2)$ as solutions of the following equations:

For $\psi_{1}^{n+1}$:

[TABLE]

For $\psi_{2}^{n+1}$:

[TABLE]

For $\phi_{1}^{n+1}$:

[TABLE]

For $\phi_{2}^{n+1}$:

[TABLE]

Then for given scalar number $\xi,$ the following field functions solve the system consisting of equations (4.17), (4.19) and (4.10c):

[TABLE]

where $(\psi_{i}^{n+1},\phi_{i}^{n+1})$ $(i=1,2)$ are given by equations (4.20a)-(4.23).

Now we are ready to determine the unknown scalar $\xi.$ Following equations (2.43), we define

[TABLE]

where equation (4.17b) has been used. Accordingly, in light of equations (4.1b) and (2.38), we define

[TABLE]

We further define

[TABLE]

Combining equations (4.10d) and (4.14), we obtain the formula for $\xi$,

[TABLE]

where $S_{0}$ is given by

[TABLE]

Once $\xi$ is known, $\phi^{n+1}$ and $\psi^{n+1}$ can be obtained directly by equation (4.24) and $R^{n+1}$ can be computed based on equation (2.46).

Equations (4.17a)-(4.23) are Helmholtz type equations with Neumann type boundary conditions. They can be implemented with $C^{0}$ spectral elements in a straightforward fashion.

Remark 4.1**.**

In equation (4.10a), we have treated the nonlinear term explicitly by $h(\bar{\phi}^{n+1})$. When $\Delta{t}$ becomes large, $\bar{\phi}^{n+1}$ can no longer approximate $\phi^{n+1}$ well. Thus, although the scheme (4.10)-(4.11) is unconditionally stable, the simulation will lose accuracy for large time steps. One possible approach to improve the accuracy is to replace $\xi h(\bar{\phi}^{n+1})$ in equation (4.10a) by

[TABLE]

where $\phi_{0}$ is a chosen field function close to $\phi^{n+1}$, e.g. a snapshot of the $\phi$ field in the recent past. The first term in the above equation serves as a linearized approximation of $h(\phi^{n+1})$ and the second term serves as a correction to this approximation. By doing so, equation (4.10a) with the mentioned modification is still linear, but can no longer be decoupled straightforwardly. One needs to solve either a fourth-order linear equation or a coupled linear system. However, this treatment can result in improved accuracy besides unconditional stability. We will demonstrate this in the forthcoming case for the Cahn-Hilliard equation with variable mobility.

4.2 Variable Mobility

Next, we consider the case with a variable mobility, $m(\phi)=\max(m_{0}(1-\phi^{2}),0)$. We reformulate the equations (4.1a)–(4.1c) into

[TABLE]

In these equations, $\mu$ is given by (4.1b), $\phi_{0}$ is a chosen field distribution corresponding to $\phi(\bm{x},t)$ at a certain time instant or at some time instants, and

[TABLE]

where $S\geqslant 0$ is a chosen constant. By incorporating the following zero terms into the RHS of (4.6),

[TABLE]

we can transform this equations into,

[TABLE]

Following equations (2.25a)-(2.25e), we propose the following scheme:

[TABLE]

and

[TABLE]

together with the boundary condition (4.10c) and the initial condition (4.12). In these equations, $\left.\frac{\partial\phi}{\partial t}\right|^{n+1}$ and $\left.\frac{dR}{dt}\right|^{n+1}$ are defined in (2.21b), $\bar{\phi}^{n+1}$ is given by (2.21c), and $\bar{C}^{n+1}$ and $\bar{\mu}^{n+1}$ are computed by

[TABLE]

$\tilde{\phi}^{n+1}$, $\tilde{\phi}^{n+3/2}$, $\tilde{\mu}^{n+1}$, and $\left.\frac{\partial\phi}{\partial t}\right|^{*,n+1}$ are approximations to be specified later.

Theorem 4.2**.**

In the absence of the external source term ($f=0$), and with zero boundary conditions ($d_{a}=d_{b}=0$), the scheme consisting of (4.34)-(4.35) is unconditionally energy stable in the sense that

[TABLE]

if the approximation of $R(t)$ at time step $\frac{1}{2}$ is positive.

Proof.

We take the $L^{2}$ inner product between \big{(}-\lambda\nabla^{2}\phi^{n+1}+h(\phi^{n+1})\big{)} and equation (4.34a), and add the resultant equation to equation (4.35). This leads to

[TABLE]

By the same arguments as in the proof of Theorem 4.1, we arrive at the relation (4.37) based on the above equation. ∎

For implementation of the scheme, one notes that equation (4.34a) can be transformed into

[TABLE]

Barring the unknown scalar $\xi,$ equations (4.39), (4.34b), (4.34e) and (4.10c) can be solved as follows. Introduce two pairs of field functions $(\phi_{i}^{n+1},C_{i}^{n+1})$ $(i=1,2),$ as the solution of the following equations:

For $(\phi_{1}^{n+1},C_{1}^{n+1})$:

[TABLE]

For $(\phi_{2}^{n+1},C_{2}^{n+1})$:

[TABLE]

Then for given scalar value $\xi,$ the following field functions solve the system consisting of equations (4.34a)-(4.34e) and (4.10c):

[TABLE]

where $(C_{i}^{n+1},\phi_{i}^{n+1})$ $(i=1,2)$ are given by equations (4.40)-(4.41).

The unknown scalar value $\xi$ remains to be determined. Following equation (2.43), $\tilde{\phi}^{n+1},$ $\tilde{\mu}^{n+1}$ and $\left.\frac{\partial\phi}{\partial t}\right|^{*,n+1}$ are again given by equations (4.25) and (4.26), where based on equation (4.34b) we compute $\nabla^{2}\tilde{\phi}^{n+1}$ by

[TABLE]

The approximation $\tilde{\phi}^{n+\frac{3}{2}}$ is given by (4.27). As a result, $\xi$ can be computed by,

[TABLE]

where $S_{0}$ is given by (4.29), and $\phi^{n+1}$ and $R^{n+1}$ can be evaluated by equations (4.42) and (2.46), respectively.

Equations (4.40)-(4.41) can be discretized in space by $C^{0}$ spectral elements, and their weak forms are:

For $(\phi_{1}^{n+1},C_{1}^{n+1})$: Find $\phi_{1}^{n+1}\;,C_{1}^{n+1}\in H^{1}(\Omega)$ such that

[TABLE]

for all $\varphi\in H^{1}(\Omega).$

For $(\phi_{2}^{n+1},C_{2}^{n+1})$: Find $\phi_{2}^{n+1}\;,C_{2}^{n+1}\in H^{1}(\Omega)$ such that

[TABLE]

for all $\varphi\in H^{1}(\Omega).$

Remark 4.2**.**

If one chooses $\kappa(\phi_{0})=0$ and $m_{c}(\phi_{0})=m_{0}>0$, then the scheme (4.34a)–(4.35) can also be implemented by solving four de-coupled Helmholtz type equations in a way similar to the constant mobility case in Section 4.1.

4.3 Numerical Results

We next provide numerical examples to demonstrate the accuracy and unconditional stability of the proposed schemes (4.10)-(4.11) and (4.34)-(4.35) for Cahn-Hilliard equation with constant and variable mobilities. For cases with variable mobility we employ $m_{c}(\phi_{0})=m(\phi_{0})=\max(m_{0}(1-\phi_{0}^{2}),0)$ in the algorithm with these tests, where $m_{0}$ and $\phi_{0}$ will be specified below.

4.3.1 Convergence Rates

Consider domain $\Omega=[0,2]\times[-1,1]$ and a contrived solution in this domain:

[TABLE]

The external force and boundary source terms $f(\bm{x},t)$, $d_{a}(\bm{x},t)$ and $d_{b}(\bm{x},t)$ in (4.1a), (4.1c) and (4.1d) are chosen such that the analytic expression (4.49) satisfies (4.1).

The computational domain $\Omega$ is discretized with two equal-sized quadrilateral elements. The algorithms (4.10)-(4.11) for the constant-mobility case and (4.34)-(4.35) for the variable-mobility case are employed to numerically integrate the Cahn-Hilliard equation from $t=t_{0}$ to $t=t_{f}$. The initial field function $\phi_{in}$ is obtained by setting $t=t_{0}$ in the contrived solution (4.49). The numerical errors are computed by comparing the numerical solution against the analytic solution (4.49) at $t=t_{f}.$ In the following convergence tests, we fix $\mathscr{F}(R)=R,$ $C_{0}=1,$ and $\phi_{0}=\phi_{in}(\bm{x})$ in (4.34). The values for the simulation parameters are summarized in Table 1.

In the spatial convergence test, we fix $\Delta{t}=0.001,$ $t_{0}=0.1$ and $t_{f}=0.2$, and vary the element order systematically from 2 to 20. The numerical errors in $L^{\infty}$ and $L^{2}$ norms at $t=t_{f}$ are then recorded. For the algorithm with constant mobility, $S$ in equation (4.10) is chosen as $S=\sqrt{\frac{4\gamma_{0}\lambda}{m_{0}\Delta t}},$ while for the algorithm with variable mobility we use $S=1.$ Figures 4.1(a) and (b) show the numerical errors as a function of the element order from these tests. It can be observed that the errors decrease exponentially with increasing element order and that the error curves level off at around $10^{-5}$ and $10^{-6}$ beyond element order 8 and 10, respectively for these two solvers, due to the saturation of temporal errors.

In the temporal convergence test, we fix the element order at a large value 18, $t_{0}=0.1,$ and $t_{f}=1.1$, and vary $\Delta{t}$ systematically from $0.2$ to $1.953125\times 10^{-4}$ to study the behavior of numerical errors. For the constant-mobility case, $S=\sqrt{\frac{4\gamma_{0}\lambda}{m_{0}\Delta t_{\min}}}$ (where $\Delta t_{\min}=10^{-4}$), while for the variable-mobility case $S=1.$ Figures 4.1(c) and (d) show the numerical errors as a function of $\Delta t$ for these cases. We observe a second-order convergence rate in time for both cases.

4.3.2 Constant Mobility: Coalescence of Two Drops

We next consider the coalescence of two drops to demonstrate the numerical properties of the proposed scheme (4.10)-(4.11) for problems with constant mobility. Consider a square domain $\Omega=[0,1]^{2}$ and two materials contained in this domain. It is assumed that the dynamics of the material regions is governed by the Cahn-Hilliard equation with a constant mobility, $m(\phi)=m_{0}>0$, and that $\phi=1$ and $\phi=-1$ correspond to the bulk of the first and second materials, respectively. We assume that at $t=0$ the first material occupies two circular regions that are right next to each other and the the rest of the domain is filled by the second material.

To be more specific, the initial distribution of the material takes the form

[TABLE]

where $\bm{x}_{0}=(x_{0},y_{0})=(0.3,0.5)$ and $\bm{x}_{1}=(0.7,0.5)$ are the centers of the circular regions for the first material, and $R_{0}=0.19$ is the radius of these circles. The external force and the boundary source terms in (4.1) are set to $f(\bm{x},t)=d_{a}(\bm{x},t)=d_{b}(\bm{x},t)=0.$ We discretize the domain using 400 equal-sized quadrilateral elements with element order 10. We employ a mapping function $\mathscr{F}(R)=R^{2}$ for this problem. The simulation parameters are listed as follows:

[TABLE]

Figure 4.2 shows the evolution of the two material regions with a temporal sequence of snapshots of the interfaces between these two materials visualized by the contour level $\phi=0.$ It can be observed that the two separate regions of the first material gradually coalescence with each other to form a single drop under the Cahn-Hilliard dynamics.

To investigate the effect of time step size on the accuracy of the simulation results, in Figure 4.3 we compare the distributions of the material interfaces at $t=50$ obtained with several time step sizes, ranging from $\Delta{t}=10^{-1}$ to $\Delta{t}=10^{-4}.$ The distribution computed with $\Delta{t}=10^{-2},$ $10^{-3}$ and $10^{-4}$ are essentially the same. With the larger time step size $\Delta{t}=10^{-1},$ some difference can be noticed in the material distribution compared with those obtained using smaller $\Delta{t}$ values. This suggests the simulation is starting to lose accuracy with time step sizes $\Delta t=10^{-1}$ and larger.

Figure 4.4 shows the time histories of the total energy $E_{tot}(t)$ (see equation (4.3)) and the ratio $\xi=\frac{\mathscr{F}(R)}{E}$ obtained using time step sizes $\Delta{t}=10^{-2}$ to $\Delta{t}=10^{-4}.$ It can be observed that the history curves essentially overlap with one another for different time step sizes. The computed values for $\xi=\frac{\mathscr{F}(R)}{E}$ are very close to 1 for each $\Delta{t}$, suggesting that $\mathscr{F}(R)$ is a good approximation for $E(t)$ and the numerical approximation is accurate with these time steps.

Thanks to the energy stability property of the current method, we can use fairly large time step sizes for the simulations. In Figure 4.5, we depict some longer time histories (up to $t=10000$) of the total energy $E_{tot}(t)$ and the ratio $\xi=\frac{\mathscr{F}(R)}{E}$ obtained using several large time step sizes $\Delta{t}=0.1,1,10.$ At these large $\Delta{t}$ values we can no longer expect the results to be accurate. Indeed, in Figure 4.5(a), $E_{tot}$ increases initially, and levels off over time at around $E_{tot}\approx 2000.$ Meanwhile, $\xi$ decreases rapidly to a smaller number close to 0, suggesting that there is a large discrepancy between $\mathscr{F}(R)$ and $E(t).$ While these computation results are not accurate, they nonetheless demonstrate the proposed method is stable and robust with large time steps.

As discussed in previous sections, the current scheme guarantees the positivity of the computed $\xi$ and $R(t)$ values, regardless of the time step size or the external forces. In Figure 4.6, we compare the time histories of the computed auxiliary variable $R(t)$ obtained using the current method and the scalar auxiliary variable (SAV) method from ShenXY2018 . In the SAV method, the auxiliary variable $R(t)$ is computed by a dynamic equation stemming from the relation $R(t)=\sqrt{E_{1}(t)}$, where $E_{1}(t)=\int_{\Omega}H(\phi)d_{\Omega}+C_{0}>0$. Therefore, $R(t)$ is expected to be positive on the continuous level. In reality, however, the discrete solutions for $R(t)$ computed by the SAV method can become negative. This is evident from Figure 4.6(b), where the result obtained using the SAV method with a large $\Delta{t}=1$ is shown. On the other hand, the discrete solutions for $R(t)$ from the current method are guaranteed to be positive, which is evident from Figure 4.6(a).

4.3.3 Variable Mobility: Evolution of a Drop

We next consider the evolution of a square drop governed by the Cahn-Hilliard equation with a variable mobility. The computational domain and the settings follow those for the coalescence of two drops discussed above. The difference lies in the initial distribution of the materials. To be precise, the initial distribution of field function is set as follows:

[TABLE]

where $(x_{0},y_{0})=(0.5,0.5)$ is the center of the domain and $h_{0}=0.2.$

Figure 4.7 shows the evolution of the system with a temporal sequence of snapshots of the interfaces between the two materials. These results are computed with a time step size $\Delta{t}=0.01$, $S=1$, $C_{0}=10^{6}$, and the mapping function $\mathscr{F}(R)=R^{2}$. The $\phi_{0}$ in the algorithm is taken as the field $\phi(\bm{x},t)$ at every fifth time step, i.e. $\phi_{0}(\bm{x})=\phi^{5k}(\bm{x})$ ($k=0,1,2\dots$). In other words, the $\phi_{0}$ field and also the coefficient matrices of the system are updated every $5$ time steps in this set of tests. These results illustrate the process for the evolution of the initial square region into a circular region under the Cahn-Hilliard dynamics.

In Figure 4.8, we show the time histories of the total energy $E_{tot}(t)$ and $\xi=\frac{\mathscr{F}(R)}{E(t)}$ obtained with several time step sizes ranging from $\Delta{t}=10^{-2}$ to $\Delta{t}=10^{-4}.$ Note that the variable mobility is $m(\phi)=\max(m_{0}(1-\phi^{2}),0)$. Here we have considered two ways to simulate the problem:

• •

by setting $\phi_{0}=0$ in the algorithm. This leads to $m_{c}(\phi_{0})=m(\phi_{0})=m_{0}$ and $\kappa(\phi_{0})=-\frac{\lambda}{\eta^{2}}$, and a time-independent coefficient matrix for the system, which can be pre-computed. We refer to this setting as the standard way.

• •

by setting $\phi_{0}=\phi^{5k}$ ($k=0,1,2,\dots$) in the algorithm. The $\phi_{0}$ field and the coefficient matrix are quasi time-independent, and they are updated every $5$ time steps.

With the smaller time step sizes $\Delta{t}=10^{-3}$ and $10^{-4}$, we set $\phi_{0}=0$ in the algorithm (the standard way) when performing simulations. With the larger $\Delta{t}=10^{-2}$, we have conducted simulations in both ways with the algorithm. In Figure 4.8(a) the results from these two settings are marked by “no update” (standard way) and “update” (second way) in the legend corresponding to $\Delta t=10^{-2}$. It is observed that the energy histories corresponding to $\Delta{t}=10^{-4}$ and $10^{-3}$, and $\Delta t=10^{-2}$ with $\phi_{0}$ updated periodically, essentially overlap with each other. However, the energy history corresponding to $\Delta{t}=10^{-2}$ with $\phi_{0}=0$ exhibits a pronounced discrepancy compared with the other cases. These results indicate that with the standard way (by setting $\phi_{0}=0$) in the algorithm the simulation result would cease to be accurate when the time step size increases to $\Delta t=10^{-2}$. However, if one uses the second way (by updating $\phi_{0}$ periodically), accurate simulation result can be obtained even with $\Delta t=10^{-2}$. In other words, by updating $\phi_{0}$ in the algorithm from time to time, one can improve the accuracy of the simulations even at larger time step sizes. We depict in Figure 4.8(b) the time histories of $\xi=\frac{\mathscr{F}(R)}{E}$ corresponding to these time step sizes. Shown for $\Delta{t}=10^{-2}$ in this plot is the result with $\phi_{0}$ updated periodically. It is observed that the computed $\xi$ is essentially 1 with $\Delta{t}=10^{-3}$ and $10^{-4}.$ With $\Delta{t}=10^{-2}$ (and $\phi_{0}$ updated periodically), the computed $\xi$ is substantially smaller than 1. But interestingly, the simulation results for the field function $\phi$ are still quite accurate with this larger $\Delta t$. This group of tests suggests that one possible way to improve the accuracy of the proposed energy-stable scheme is to update the $\phi_{0}$ in the algorithm periodically, e.g. every $N$ time steps. By choosing an appropriate $N$ for a given problem, one can enhance the simulation accuracy even at large or fairly large time step sizes. Because $\phi_{0}$ and the coefficient matrix for the system only needs to be updated infrequently, the cost associated with updating the coefficient matrix can be manageable. There is a drawback with this, however. The computations using the second way (updating $\phi_{0}$ periodically) seems not as robust as the standard way (by setting $\phi_{0}=0$) for large $\Delta t$. Because of the non-zero $\phi_{0}$ field in the algorithm, the conditioning of the system coefficient matrix using the second way seems to become worse for large $\Delta t$. We observe that for larger $\Delta{t}\geqslant 0.1$ the system coefficient matrix using the second way can become singular and the computation may break down.

5 Nonlinear Klein-Gordon Equation

We consider an energy-conserving system, the nonlinear Klein-Gordon equation, in this section and apply the gPAV method to this system. Consider the nonlinear Klein-Gordon equation Strauss1978 on a domain $\Omega$ (with boundary $\Gamma$)

[TABLE]

where $\varepsilon$, $\alpha$ and $\varepsilon_{1}$ are positive constants. These equations are supplemented by the initial conditions

[TABLE]

In these equations $g(u)=G^{\prime}(u)$ and $G(u)$ is a potential energy function with $G(u)\geqslant 0.$ The above system satisfies the following energy balance law:

[TABLE]

We define a shifted total energy according to equation (2.11),

[TABLE]

where $C_{0}$ is chosen such that $E(t)>0$. Choose $\mathscr{F}$ and $\mathscr{G}$, and define the auxiliary variable $R(t)$ based on equation (2.13a). Following equation (2.14), we have

[TABLE]

where integration by part has been used.

Following equations (2.17)-(2.19), we reformulate equations (5.2) and (5.7) into

[TABLE]

Note that when deriving (5.8b) we have incorporated the following zero terms to the RHS,

[TABLE]

The reformulated system consists of equations (5.1), (5.8a)-(5.8b) and (5.3)-(5.4), which is equivalent to the original system (5.1)-(5.4).

Since the Klein-Gordon equation is conservative (in the absence of external source term and with appropriate boundary condition), we will employ the Crank-Nicolson method for time discretization of the field variables, by enforcing the discretized equations at step $(n+1/2)$. This corresponds to the approximations (2.33a)–(2.34) with $\theta=\frac{1}{2}$ and $\beta=0$. So the method here is slightly different than the one presented in Section 2.2, which corresponds to $\theta=1$ and $\beta=\frac{1}{4}$ in the approximations (2.33a)–(2.34). The energy-stable scheme for the nonlinear Klein-Gordon equation is then as follows:

[TABLE]

together with

[TABLE]

These equations are supplemented by the initial conditions

[TABLE]

where $E^{0}$ is evaluated by

[TABLE]

In the above equations, $D_{\mathscr{F}}(R)|^{n+\frac{1}{2}}$ is defined by (2.34) with $\theta=1/2$, and

[TABLE]

$\tilde{u}^{n+1},$ $\tilde{v}^{n+1},$ $\tilde{u}^{n+\frac{1}{2}}$ and $\tilde{v}^{n+\frac{1}{2}}$ are second-order approximations of $u^{n+1},$ $v^{n+1},$ $u^{n+\frac{1}{2}}$ and $v^{n+\frac{1}{2}},$ respectively, defined later in (5.22)-(5.23).

Theorem 5.1**.**

In the absence of the external force $f=0,$ and with homogeneous boundary condition ($d_{a}=0$) and suppose that the initial condition $v_{\rm in}$ satisfies the compatibility condition $v_{in}|_{\Gamma}=0,$ the scheme consisting of (5.9)-(5.11) conserves the modified energy $\mathscr{F}(R)$ in the sense that:

[TABLE]

Proof.

Multiplying \big{(}-\alpha^{2}\nabla^{2}u^{n+\frac{1}{2}}+\varepsilon_{1}^{2}u^{n+\frac{1}{2}}+g(u^{n+\frac{1}{2}})\big{)} to equation (5.9a), $\varepsilon^{2}v^{n+\frac{1}{2}}$ to equation (5.9b), taking the $L^{2}$ integrals, and summing up the resultant equations with equation (5.10), we arrive at the relation,

[TABLE]

where we have used equations (2.21b)-(2.22). If $d_{a}=0$, then $u^{n}|_{\Gamma}=0$ and $v^{n}|_{\Gamma}=0$ for all $n>0$. Based on the definition of $\tilde{v}^{n+\frac{1}{2}}$ in the equation (5.23) below, it is straightforward to verify that $\tilde{v}^{n+\frac{1}{2}}|_{\Gamma}=0$ as long as $v^{0}|_{\Gamma}=0$. Furthermore, if $f=0,$ the volume integrals in equation (5.15) vanish. This leads to equation (5.14). ∎

Remark 5.1**.**

Since $\mathscr{F}(R)$ is an approximation of $E(t),$ the discrete conservation for $\mathscr{F}(R)$ in equation (5.14) does not imply the conservation for $E(t)$ on the discrete level. However, it does lead to an unconditionally energy stable scheme for long time simulations.

Despite the complication caused by the unknown scalar variable $\xi,$ the proposed scheme can be solved in a decoupled fashion. Combining equations (5.9a) and (5.13), we get

[TABLE]

Inserting equation (5.16) into (5.9b) leads to

[TABLE]

To solve this equations, we introduce $u_{1}^{n+1}$ and $u_{2}^{n+1}$ as solutions of the following two equations:

[TABLE]

and

[TABLE]

Then the solution to equation (5.17), together with the boundary condition (5.9e), is given by

[TABLE]

where $\xi$ is to be determined.

We define

[TABLE]

By combining equations (5.9c) and (5.15), we can determine $\xi$,

[TABLE]

With $\xi$ known, $u^{n+1}$ and $v^{n+1}$ can be computed by equations (5.21) and (5.16), respectively. $R^{n+1}$ can be computed by,

[TABLE]

The weak formulations for equations (5.18) and (5.20) are: Find $(u_{1}^{n+1},u_{2}^{n+1})\in H^{1}(\Omega)$ such that

[TABLE]

[TABLE]

These can be implemented with $C^{0}$ spectral elements in a straightforward fashion.

5.1 Numerical Results

We next provide numerical examples to demonstrate the accuracy and unconditional stability of the proposed scheme to the Klein-Gordon equation (5.1)-(5.3). Specifically, we fix the parameters therein and the potential energy function as

[TABLE]

This corresponds to the dimensionless relativistic Sine-Gordon equation (DRSG) (see e.g. Bao2012 ).

5.1.1 Convergence Rates

To study the convergence rates in space and time of the proposed method, we employ the following manufactured analytic solution

[TABLE]

The external force $f(\bm{x},t)$ in (5.2) and the external boundary source term $d_{a}(\bm{x},t)$ are chosen such that the above expression (5.29) satisfies equations (5.1)-(5.3).

The computational domain $\Omega=[0,2]\times[-1,1]$ is discretized using two equal-sized quadrilateral elements, with the element order and the time step size $\Delta{t}$ varied systematically in the spatial and temporal tests. The algorithm presented in this section is employed to numerically integrate the DRSG equation from $t=0$ to $t=t_{f}.$ The mapping $\mathscr{F}(R)=R$ and $C_{0}=1$ are used in these computations. The initial condition $u_{in}$ and $v_{in}$ are obtained by setting $t=0$ in the analytic expression (5.29) and using (5.1). We then record the numerical errors in different norms by comparing the numerical solution with the analytic solution at $t=t_{f}.$

To conduct the spatial convergence test, we vary systematically the element order from 2 to 20 and depict in Figure 5.1(a) the $L^{\infty}$ and $L^{2}$ errors of $u$ as a function of the element order with a fixed $\Delta{t}=0.001$ and $t_{f}=0.1.$ It is observed that the numerical errors decay exponentially with increasing element order, and levels off beyond element order 12, caused by the saturation of temporal errors.

To study the temporal convergence rate, we fix the element order at a large value 18 and $t_{f}=1.0$. The time step size $\Delta{t}$ is varied systematically from 0.2 to $7.8125\times 10^{-4}$ and the numerical errors in $L^{\infty}$ and $L^{2}$ norms are depicted in Figure 5.1(b). A second-order convergence rate in time is clearly observed.

5.1.2 Study of Method Properties

We next study the remarkable stability of the proposed method with the DRSG equation. Consider the DRSG equation on the domain $\Omega=[0,14]^{2}$, with zero external force $f(\bm{x},t)=0$ and zero boundary source term $d_{a}(\bm{x},t)=0$ in (5.3). The initial conditions are set to

[TABLE]

With these initial and boundary conditions, the DRSG equation is energy conserving.

The domain $\Omega$ is discretized with 400 equal-sized quadrilateral elements with a fixed element order 10. We employ a mapping function $\mathscr{F}(R)=\frac{e_{0}}{2}\ln(\frac{\kappa_{0}+R}{\kappa_{0}-R})$ ($e_{0}=10,\,\kappa_{0}=100$) and the energy constant $C_{0}=1$ in the algorithm. Figure 5.2 illustrates the evolution of $u$ by a sequence of snapshots of its contour levels. One can observe a circular wave pattern starting from the center of the domain and propagating outward toward the boundaries. As the wave reaches the boundaries, the interaction with the Dirichlet boundary ($u=0$) gives rise to an extremely complicated wave pattern; see Figure 5.2(d).

Figure 5.3(a) shows the time histories of the energy errors, $|E(t)-E(0)|$, obtained using several time step sizes ($\Delta t=10^{-4}$, $10^{-3}$ and $10^{-2}$). One can observe oscillations in the history curves about their respective mean values that are consistent with a second order accuracy in time. It should again be noted that the current algorithm conserves the modified energy $\mathscr{F}(R)$ discretely, not the original energy $E(t)$. Figure 5.3(b) shows time histories of the ratio $\xi=\frac{\mathscr{F}(R)}{E}$ corresponding to these $\Delta t$ values. The computed $\xi$ values are essentially $1$, indicative of the accuracy of these simulations.

We then increase the time step size to $\Delta{t}=0.1,1$ and $10$, and depict in Figure 5.4(a) the time histories of $E(t)$ and $\mathscr{F}(R)$ for a long time simulation to $t=1000$. Large discrepancies between the energy $E(t)$ and $\mathscr{F}(R)$ can be observed, especially for $\Delta{t}=1$ and $10$, suggesting that $\mathscr{F}(R)$ no longer approximates well the energy $E(t)$ with these time step sizes. Note that the $\mathscr{F}(R)$ histories obtained by different large $\Delta t$ values overlap with one another. This is consistent with Theorem 5.1 that the current scheme conserves the modified energy $\mathscr{F}(R)$. It can be observed from Figure 5.4(b) that the computed $\xi=\frac{\mathscr{F}(R)}{E}$ becomes significantly smaller than 1, indicative of large errors in the simulations with these large time step sizes. However, the computations are evidently stable, even with these large $\Delta t$ values.

6 Concluding Remarks

In this paper we have presented a framework (gPAV) for developing unconditionally energy-stable schemes for general dissipative systems. The scheme is based on a generalized auxiliary variable (which is a scalar number) associated with the energy functional of the system. We find that the square root function, which is critical to previous auxiliary-variable approaches, is not essential to devising energy-stable schemes. In the current method, the auxiliary variable can be defined by a rather general class of functions, not limited to the square-root function. The gPAV method is applicable to general dissipative systems, and a unified procedure for discretely treating the dissipative governing equations and the generalized auxiliary variable has been presented. The discrete energy stability of the proposed scheme has been proven for general dissipative systems. The presented method has two attractive properties:

• •

The scheme requires only the solution of linear algebraic equations within a time step, and no nonlinear solver is needed. Furthermore, with appropriate choice of the $\bm{F}_{L}$ operator in the algorithm, the resultant linear algebraic systems upon discretization involve only constant and time-independent coefficient matrices, which only need to be computed once and can be pre-computed. In terms of computational cost, the scheme is computationally very competitive and attractive.

• •

The generalized auxiliary variable can be computed directly by a well-defined explicit formula. The computed values for the auxiliary variable are guaranteed to be positive, regardless of the time step size or the external forces or source terms.

Three specific dissipative systems (a chemo-repulsion model, Cahn-Hilliard equation with constant and variable mobility, and the nonlinear Klein-Gordon equation) have been studied in relative detail to demonstrate the gPAV framework developed herein. Ample numerical experiments have been presented for each system to demonstrate the performance of the method, the effects of algorithmic parameters, and the stability of the scheme with large time step sizes.

All physically meaningful systems in the real world are energy dissipative (or conserving) due to the second law of thermodynamics, and these systems are typically nonlinear. The design of energy-stable and computationally-efficient schemes for such systems is critical to their numerical simulations, and this is in general a very challenging task. The gPAV framework presented here lays out a roadmap for devising discretely energy-stable schemes for general dissipative systems. The computational efficiency (e.g. involving linear equations with pre-computable coefficient matrices) and the guaranteed positivity of the computed auxiliary variable of the method are particularly attractive, in the sense that the gPAV method is not only unconditionally energy-stable but also can be computationally efficient and competitive. We anticipate that the gPAV method will be useful and instrumental in numerical simulations of a number of computational science and engineering disciplines.

Acknowledgement

This work was partially supported by NSF (DMS-1522537).

Appendix A. Approximation for the First Time Step

We present a method on how to deal with the first time step such that the approximation for the auxiliary variable $R(t)$ at time step $\frac{1}{2}$ shall be positive. We consider below only the formulation based on $\frac{\mathscr{F}(R)}{E}$. It is noted that for the alternative formulation based on $\frac{R}{\mathscr{G}(E)}$ (see Section 2.4) one can modify the following scheme in a straightforward fashion to achieve the same property. The notations here follow those employed in the main text.

Consider the system consisting of equations (2.17), (2.19), the boundary condition (2.2), and the initial conditions (2.3) and (2.20). Define

[TABLE]

One notes that $E^{0}>0$ and $R^{0}>0$.

We compute the first time step in two substeps. In substep one we compute an approximation of ($\bm{u}^{1},R^{1}$), denoted by ($\bm{u}_{a}^{1},R_{a}^{1}$), and in substep two we compute the final ($\bm{u}^{1},R^{1}$). More specifically, the scheme is as follows:

Substep One:

[TABLE]

Substep Two:

[TABLE]

Note that in the above equations the superscript of a variable such as $(\cdot)^{1/2}$ and $(\cdot)^{3/2}$ denotes the time step index. In (6.2b) and (6.2e) $\tilde{\bm{u}}_{a}^{1}$ is an approximation of $\bm{u}_{a}^{1}$ and will be specified later in (6.12). In (6.3e) $\tilde{\bm{u}}^{1}$ is an approximation of $\bm{u}^{1}$ and will be specified later also in (6.12). In (6.3b), (6.3c) and (6.3e), $\tilde{\bm{u}}^{3/2}$, $R^{1/2}$ and $R^{3/2}$ are defined by

[TABLE]

It can be noted that the above scheme represents a first-order approximation of ($\bm{u}^{1},R^{1}$) for the first time step.

Combine equations (6.2a) and (6.2e) and we have

[TABLE]

where $S_{a}=\int_{\Omega}V_{s}(\bm{f}^{1},\tilde{\bm{u}}_{a}^{1})d\Omega+\int_{\Gamma}B_{s}(\bm{f}_{b}^{1},\tilde{\bm{u}}_{a}^{1})d\Gamma.$ In light of (6.2b), this leads to

[TABLE]

Since $R^{0}>0$, we conclude that $\xi_{a}>0$ and $R_{a}^{1}>0$ based on these equations. It follows that $R^{1/2}=\frac{1}{2}(R_{a}^{1}+R^{0})>0$ in light of equation (6.4).

Similarly, combining equations (6.3a) and (6.3e) gives rise to

[TABLE]

where $S_{0}=\int_{\Omega}V_{s}(\bm{f}^{1},\tilde{\bm{u}}^{1})d\Omega+\int_{\Gamma}B_{s}(\bm{f}_{b}^{1},\tilde{\bm{u}}^{1})d\Gamma.$ In light of (6.3b) and (6.4), we have

[TABLE]

We therefore conclude that $\xi>0$, $R^{3/2}>0$ and $R^{1}>0$.

We still need to determine $\bm{u}_{a}^{1}$ and $\bm{u}^{1}$, and specify $\tilde{\bm{u}}_{a}^{1}$ and $\tilde{\bm{u}}^{1}$. Note that $\bm{F}_{L}(\bm{u})$ and $\bm{B}(\bm{u})$ are linear operators. Equations (6.2a) and (6.2d), and also equations (6.3a) and (6.3d), can be solved as follows. Define two variables $\bm{u}_{1}^{1}$ and $\bm{u}_{2}^{1}$ as solutions to the following systems, respectively:

For $\bm{u}_{1}^{1}$:

[TABLE]

For $\bm{u}_{2}^{1}$:

[TABLE]

Then it is straightforward to verify that, for given $\xi_{a}$ and $\xi$, the following functions respectively solve the equations (6.2a) and (6.2d), and equations (6.3a) and (6.3d),

[TABLE]

We then specify $\tilde{\bm{u}}_{a}^{1}$ and $\tilde{\bm{u}}^{1}$ as follows,

[TABLE]

The solution for ($\bm{u}^{1},R^{1}$) at the first time step consists of the following procedure:

• •

Solve equations (6.9a)–(6.9b) for $\bm{u}_{1}^{1}$;

Solve equations (6.10a)–(6.10b) for $\bm{u}_{2}^{1}$.

• •

Compute $\tilde{\bm{u}}_{a}^{1}$ and $\tilde{\bm{u}}^{1}$ by equation (6.12);

Compute $\xi_{a}$ and $R_{a}^{1}$ by equation (6.6);

Compute $\bm{u}_{a}^{1}$ by equation (6.11a).

• •

Compute $\tilde{\bm{u}}^{3/2}$ and $R^{1/2}$ based on equation (6.4);

Compute $\xi$ and $R^{1}$ based on equation (6.8);

Compute $\bm{u}^{1}$ by equation (6.11b).

We can make the following conclusion based on the above discussions.

Theorem 6.1**.**

The scheme represented by (6.2a)–(6.3e) for computing the first time step has the property that

[TABLE]

where $R^{1/2}$ and $R^{3/2}$ are given by (6.4), regardless of the time step size $\Delta t$ and the external forces $\bm{f}$ and $\bm{f}_{b}$.