Scalar auxiliary variable approache for the surface quasi-geostrophic equation

Abstract

We present a numerical method for approximating the two-dimensional surface quasi-geostrophic equation. We first reformulate the equation into an equivalent system by using the scalar auxiliary variable approach. Then we propose first-order and second-order discretization schemes for solving the new surface quasi-geostrophic system in time. Furthermore, we prove that both of these schemes based on the scalar auxiliary variable approach satisfy an unconditionally energy stability property. Since the method gives two linear equations with constant coefficients in each time step. We have an efficient approach and accurate for solving these spacial fractional diffusion equations. Ample numerical experiments are carried out to validate the correctness of these schemes and their accuracy for inviscid problems and the problems with small viscosity.

Keywords

Scalar auxiliary variable unconditionally stable fractional reaction-diffusion equation

Introduction

In this article, we consider a numerical approximation of the two-dimensional (2D) surface quasi-geostrophic (SQG) equation. The numerical schemes are based on a modified SQG equations by using the scalar auxiliary variable (SAV) approach. The general three-dimensional (3D) quasi-geostrophic equation, first derived by Charney in 1940s,¹ have been very successful in describing major features of large-scale motions in the atmosphere and the oceans in the mid-latitudes.^2,3 These 3D equations can be reduced to the SQG equation with uniform potential, modeling the potential temperature on the 2D boundaries^4,5 as

θ_{t} + u \cdot \nabla θ + κ (- Δ)^{α} θ = 0, x = (x, y) \in Ω

(1)

where

κ \geq 0

and

α > 0

are parameters,

Ω \in R^{2}

is a bounded periodic domain,

θ (x, t)

is a scalar representing the potential temperature, and

u = (u_{1}, u_{2})

is the velocity field determined from

θ (x, t)

by the stream function

ψ (x, t)

via the auxiliary relations

u = (u_{1}, u_{2}) = (- \partial_{y} ψ, \partial_{x} ψ), (- Δ)^{\frac{1}{2}} ψ = - θ

(2)

The fractional Laplacian

(- Δ)^{α}

in this paper is defined as follows

(- Δ)^{α} θ (x, t) = \sum_{i = 1}^{\infty} λ_{i}^{α} c_{i} (t) ϕ_{i} (t)

(3)

where

(λ_{i}, ϕ_{i})_{i}^{\infty}

are the eigenpairs of the standard Laplacian

- Δ

and

θ

has the expansion

θ (x, t) = \sum_{i = 1}^{\infty} c_{i} (t) ϕ_{i} (t)

. Alternatively, the fractional operator

(- Δ)^{α}

,⁶ with

α \in R

can be defined by

\hat{(- Δ)^{α}} f (k) = | k |^{2 α} \hat{f (k)}

where

\hat{f}

denotes the Fourier transform of

f

When the fractional order $α = \frac{1}{2}$ , equation (1) derived from the more general quasi-geostrophic models⁷ describes the evolution of the temperature on the 2D boundary of a rapidly rotating half-space with small Rossby and Ekman numbers. The turbulence theory for these special quasi-geostrophic flows has been studied earlier by Blumen⁸ and more recently by Pierrehumbert et al.⁹ Some qualitative features of the solutions of these equations in a geophysical context are developed by Held et al.⁵ The equations are an important example of a 2D active scalar.¹⁰ Furthermore, there are certain similar features between the level sets of solutions of the 2D SQG equation and the vortex lines of the 3D Euler equation.⁴ The study of the singular front formation of SQG shed lights on the behavior of the 3D Euler equation.

Quasi-geostrophic equation has been intensively investigated because of both its mathematical importance and its potential for applications in meteorology and oceanography. On application, 2D and 3D quasi-geostrophic equations are important mathematical models for simplifying 2D and 3D atmospheres and oceans. These models focus on highlighting certain features of the atmospheric and ocean. Theoretically, the existence of weak and strong solutions of subcritical quasi-geographic equations has been solved, but the uniqueness of weak solutions is still an open problem. Constantin et al.⁴ demonstrated the smooth solution of the SQG equation is unique but, if $κ = 0$ , it is known to exist only for a finite time. They also proved the global existence of a smooth solution of critical quasi-geostrophic equation when $‖ θ_{0} ‖_{L^{\infty}}$ is sufficiently small.¹¹ Hmidi and Keraani¹² gave the existence and uniqueness of the solution of the supercritical quasi-geostrophic equation. Whether solutions of the SQG equation can develop singularities with smooth initial conditions were concerned to many researchers. The global regularity for general data remains an open problem.^13,14 Numerical methods are fundamental importance to assess the behavior of the solutions for the SQG system. There are existing numerical algorithms for solving the SQG system. Constantin et al.^4,13 presented a pseudo-spectral method with the high wavenumber Fourier modes filtered out to avoid the heavy oscillations created by the singularity of the solution of the SQG equations. Meerschaert proved an implicit Euler method based on a modified Grünwald-Letnikov is unconditionally stable.¹⁵ Song introduced the concept of fractional spectral vanishing viscosity (fSVV) to conservations laws that govern the evolution of steep fronts. Bonito raised an approach based on standard finite element discretization with nonlinear stabilization,¹⁶ but this approach needs to satisfy the Courant–Friedrichs–Lewy condition. In spite of these schemes having a low computational cost and desired accuracy, they are either conditionally stable or require to solve nonlinear problems in each time step.

In this work, we present numerical schemes for the SQG equation. The significant feature of these schemes is an auxiliary variable is introduced into the SQG equation as a scalar number. This idea is inspired by Shen et al.,^17–19 working for gradient-type dynamical systems or a general dissipative system and has also been used for solving incompressible Navier-Stokes flow.^20,21 Since its introduction, the SAV approach has received much attention due to its efficiency and flexibility. The original SAV approach¹⁷ has been expanded numerous forms, including the Runge-Kutta SAV (RK-SAV) method²² in which it is shown that the SAV approach coupled with diagonally implicit high-order Runge-Kutta method is unconditionally energy stable, generalized SAV (G-SAV)²³ allow a range of functions in the definition of the SAV variable, relaxed-SAV (RSAV)²⁴ method penalizes the numerical errors of the auxiliary variables by a relaxation technique, and the Lagrange multiplier approach²⁵ which can conserve the original energy instead of a modified energy in the SAV approach. The SQG equation is reformulated into an equivalent system employing the auxiliary energy variable, and numerical schemes are devised to discrete the reformulated system. The algorithm requires only solving two linear equations with constant coefficients in each time step. In addition, these schemes satisfy discrete energy stability properties and accuracy.

The remainder of this article is organized as follows. In the “Numerical methods and stability analysis” section, we introduce an auxiliary energy variable and reformulate the SQG equation into a new system based on this variable. First-order and second-order numerical schemes are presented for approximating the reformulated equivalent system. We then prove the unconditional stability property of the proposed schemes and show the implementation details of the schemes. In the “Numerical results” section, we present the numerical results. We consider the inviscid and viscous SQG equation with smooth initial conditions, and investigate the effectiveness and accuracy of the SAV method. Finally, we provide a short summary in the “Summary” section. In Appendix A, we provide the details for the artificial term (fSVV term) of the schemes.

Numerical methods and stability analysis

In this section, we present the SAV approach for the SQG equation, and prove the first-order and second-order schemes based on the SAV approach are both unconditional stable. In addition, we also show that the algorithm can be implemented efficiently and easily.

Auxiliary variable approach

The schemes are based on an auxiliary variable approach for the SQG equation system. To this end, we first reformulate the equation into the equivalent system by introducing the following auxiliary variable

q (t) = e^{- \frac{t}{T}}, t \in [0, T], T > 0

(4)

The above scalar function will serve as the scalar auxiliary variable. Furthermore, by using the fact

\int_{Ω} (u \cdot \nabla) θ \cdot θ d x = 0

We can rewrite equation (1) into the following equivalent form

{\begin{matrix} \frac{\partial θ}{\partial t} + κ (- Δ)^{α} θ + \frac{q (t)}{e^{-} \frac{t}{T}} u \cdot \nabla θ = 0 \\ \frac{d q}{d t} = - \frac{1}{T} q + \frac{1}{e^{-} \frac{t}{T}} \int_{Ω} (u \cdot \nabla) θ \cdot θ d x \end{matrix}

(5)

Taking

L^{2}

-inner products of the above equations with

θ

and

2 q (t)

, respectively, gives

\frac{d}{d t} {‖ θ^{2} ‖}_{0}^{2} = - 2 κ {‖ (- Δ)^{\frac{α}{2}} θ ‖}_{0}^{2}

(6)

This implies that the new system can be obtained the same energy dissipation law as the original SQG equation.

From the new reformulation equation (5), we are able to construct various time discretization schemes for equations (1) and (2), and analyze the stability of the schemes. Let $n \geq 0$ denotes the time step index, $Δ t$ is the size of time step, and $t_{n} = n Δ t$ . In particular, we use $‖ \cdot ‖$ to denote the norm in $L^{2} (Ω)$ . Besides, $(\cdot, \cdot)$ is used to denote the inner product in $L^{2} (Ω)$ . The vector functions and vector spaces will be indicated by bold type.

Time discretization schemes and their stability analysis

In this subsection, we will propose two time discretization schemes for the new SQG system. We will also prove that both of these two schemes follow the discrete energy law. $S c h e m e I :$ First-order SAV-BDF scheme

The first-order backward differentiation scheme for equations (2) and (5) can be written as follows: Find $(θ^{n + 1}, q^{n + 1})$ by

\begin{aligned} \frac{θ^{n + 1} - θ^{n}}{Δ t} & + κ (- Δ)^{α} θ^{n + 1} + \frac{q^{n + 1}}{e^{-} \frac{t^{n + 1}}{T}} u^{n + 1} \\ \cdot \nabla θ^{n} + ϵ S^{β} θ^{n + 1} = 0 \end{aligned}

(6a)

\frac{q^{n + 1} - q^{n}}{Δ t} = - \frac{1}{T} q^{n + 1} + \frac{1}{e^{-} \frac{t^{n + 1}}{T}} (u^{n + 1} \cdot \nabla θ^{n}, θ^{n + 1})

(6b)

u^{n + 1} = (- \partial_{y} ψ^{n + 1}, \partial_{x} ψ^{n + 1}), (- Δ)^{\frac{1}{2}} ψ^{n + 1} = - θ^{n}

(6c)

where

θ^{n}

is an approximation of

θ (t_{n})

at time

t = t_{n}

Remark 1

The extra term $ϵ S^{β} θ^{n + 1}$ add in equation (6a) is a crucial linear stabilization term to enhance the stability and keep the required accuracy.^26,24 Relevant definitions are shown in Appendix A.

Now we consider the stability of the above time scheme with periodic boundary conditions. The discrete energy law of the first-order scheme is then given in the following theorem.

Theorem 1

The scheme (6) with periodic boundary conditions is unconditionally stable in the sense that the following discrete energy law holds

\begin{aligned} E^{n + 1} & - E^{n} ⩽ - 2 κ Δ t {‖ (- Δ)^{\frac{α}{2}} θ^{n + 1} ‖}^{2} - 2 ϵ Δ t ‖ S^{\frac{β}{2}} θ^{n + 1} ‖^{2}, \\ n \geq 0 \end{aligned}

(7)

where the energy is denoted as follows

E^{n + 1} = {‖ θ^{n + 1} ‖}^{2} + {| q^{n + 1} |}^{2}

(9)

Proof: Taking $L^{2}$ -inner products of (6a) with $2 Δ t θ^{n + 1}$ and using the identity

2 a^{k + 1} (a^{k + 1} - a^{k}) = {| a^{k + 1} |}^{2} + {| a^{k + 1} - a^{k} |}^{2} - {| a^{k} |}^{2}

(8)

we can obtain

\begin{aligned} {‖ θ^{n + 1} ‖}^{2} - {‖ θ^{n} ‖}^{2} + {‖ θ^{n + 1} - θ^{n} ‖}^{2} + 2 Δ t \frac{q^{n + 1}}{e^{-} \frac{t^{n + 1}}{T}} \\ (u^{n + 1} \cdot \nabla & θ^{n}, θ^{n + 1}) = - 2 κ Δ t {‖ (- Δ)^{\frac{α}{2}} θ^{n + 1} ‖}^{2} - 2 ϵ Δ t ‖ S^{\frac{β}{2}} θ^{n + 1} ‖^{2} \end{aligned}

(9)

Multiplying (6b) by

2 q^{n + 1} Δ t

and using the above identity, we have

{| q^{n + 1} |}^{2} - {| q^{n} |}^{2} + {| q^{n + 1} - q^{n} |}^{2} = - \frac{2}{T} Δ t {| q^{n + 1} |}^{2} + 2 Δ t \frac{q^{n + 1}}{e^{-} \frac{t^{n + 1}}{T}} (u^{n + 1} \cdot \nabla θ^{n}, θ^{n + 1})

(10)

Then combining (9) and (10), we can get

\begin{aligned} {‖ θ^{n + 1} ‖}^{2} & - {‖ θ^{n} ‖}^{2} + {| q^{n + 1} |}^{2} - {| q^{n} |}^{2} + {‖ θ^{n + 1} - θ^{n} ‖}^{2} \\ + {| q^{n + 1} - q^{n} |}^{2} + \frac{2}{T} Δ t {| q^{n + 1} |}^{2} \\ = - 2 κ Δ t {‖ (- Δ)^{\frac{α}{2}} θ^{n + 1} ‖}^{2} - 2 ϵ Δ t ‖ S^{\frac{β}{2}} θ^{n + 1} ‖^{2} \end{aligned}

(11)

and dropping the influential positive terms, we obtain (7).

Now we present the details for solving the time discretization scheme (7). Firstly, we denote $λ^{n + 1} = e^{\frac{t^{n + 1}}{T}} q^{n + 1}$ and let

θ^{n + 1} = θ_{1}^{n + 1} + λ^{n + 1} θ_{2}^{n + 1}

(12)

Plugging (12) into the scheme (6a), we can find that

θ_{i}^{n + 1} (i = 1, 2)

satisfy

{\begin{matrix} θ_{1}^{n + 1} + Δ t κ (- Δ)^{α} θ_{1}^{n + 1} + Δ t ϵ S^{β} θ_{1}^{n + 1} = θ^{n} \\ θ_{2}^{n + 1} + Δ t κ (- Δ)^{α} θ_{2}^{n + 1} + Δ t ϵ S^{β} θ_{2}^{n + 1} = - Δ t u^{n + 1} \cdot \nabla θ^{n} \end{matrix}

(13)

Obviously, we only need to solve a sequence of linear equations with constant coefficients in each time step. We can implement this method easily and very efficiently for solving the SAV system.

Once $θ_{i}^{n + 1} (i = 1, 2)$ are solved, we can determine explicitly $λ^{n + 1}$ from (6b) as follows

\begin{aligned} (\frac{T + Δ t}{T Δ t} - e^{\frac{2 t^{n + 1}}{T}} (u^{n + 1} \cdot \nabla θ^{n}, θ_{2}^{n + 1})) e^{- \frac{t^{n + 1}}{T}} λ^{n + 1} \\ = e^{\frac{t^{n + 1}}{T}} (u^{n + 1} \cdot \nabla θ^{n}, θ_{1}^{n + 1}) + \frac{q^{n}}{Δ t} \end{aligned}

(14)

According to the above description, we are able to calculate the numerical result of the subsequent time with the initial condition

θ^{0} (x)

It is clear that scheme (7) is only of first-order. Moreover, we can also construct a second-order SAV Crank–Nicolson scheme as follows.

$S c h e m e I I :$ Second-order SAV-CN scheme

The second-order semi-implicit Crank–Nicolson scheme of equation (2) and (5) can be written as follows: Find $(θ^{n + 1}, q^{n + 1})$ by

\begin{aligned} \frac{θ^{n + 1} - θ^{n}}{Δ t} & + κ (- Δ)^{α} θ^{n + \frac{1}{2}} + \frac{q^{n + \frac{1}{2}}}{{e \bar{}}^{\frac{t^{n + \frac{1}{2}}}{T}}} u^{n + \frac{1}{2}} \\ \cdot \nabla θ^{*, n + \frac{1}{2}} + ϵ S^{β} θ^{n + \frac{1}{2}} = 0 \end{aligned}

(15a)

\frac{q^{n + 1} - q^{n}}{Δ t} = - \frac{1}{T} q^{n + \frac{1}{2}} + \frac{1}{{e \bar{}}^{\frac{t^{n + \frac{1}{2}}}{T}}} (u^{n + \frac{1}{2}} \cdot \nabla θ^{*, n + \frac{1}{2}}, θ^{n + \frac{1}{2}})

(15b)

u^{n + \frac{1}{2}} = (- \partial_{y} ψ^{n + \frac{1}{2}}, \partial_{x} ψ^{n + \frac{1}{2}}), (- Δ)^{\frac{1}{2}} ψ^{n + \frac{1}{2}} = - θ^{*, n + \frac{1}{2}}

(15c)

where

θ^{n + \frac{1}{2}} : = \frac{1}{2} (θ^{n + 1} + θ^{n})

and

θ^{*, n + \frac{1}{2}} : = \frac{1}{2} (3 θ^{n} - θ^{n - 1})

. Then we have the following result.

Theorem 2

The solution of scheme (15) satisfies the following energy law:

\begin{aligned} E^{n + 1} - E^{n} ⩽ - 2 κ Δ t {‖ (- Δ)^{\frac{α}{2}} θ^{n + 1} ‖}^{2} - 2 ϵ Δ t ‖ S^{\frac{β}{2}} θ^{n + 1} ‖^{2}, \\ n \geq 0 \end{aligned}

(16)

where the energy is denoted as follows

E^{n + 1} = {‖ θ^{n + 1} ‖}^{2} + {| q^{n + 1} |}^{2}

(19)

Proof: Taking $L^{2}$ -inner products of (15a) with $2 Δ t θ^{n + \frac{1}{2}}$ and using the identity (10), we obtain

\begin{aligned} {‖ θ^{n + 1} ‖}^{2} - {‖ θ^{n} ‖}^{2} + 2 Δ t \frac{q^{n + \frac{1}{2}}}{{e \bar{}}^{\frac{t^{n + \frac{1}{2}}}{T}}} (u^{n + \frac{1}{2}} \cdot \nabla θ^{*, n + \frac{1}{2}}, θ^{n + \frac{1}{2}}) \\ = - 2 κ Δ t {‖ (- Δ)^{\frac{α}{2}} θ^{n + \frac{1}{2}} ‖}^{2} - 2 ϵ Δ t ‖ S^{\frac{β}{2}} θ^{n + \frac{1}{2}} ‖^{2} \end{aligned}

(17)

Multipying (15b) by

2 Δ t q^{n + \frac{1}{2}}

, we have

{| q^{n + 1} |}^{2} - {| q^{n} |}^{2} + \frac{2}{T} Δ t {| q^{n + \frac{1}{2}} |}^{2} = 2 Δ t \frac{q^{n + \frac{1}{2}}}{{e \bar{}}^{\frac{t^{n + \frac{1}{2}}}{T}}} (u^{n + \frac{1}{2}} \cdot \nabla θ^{*, n + \frac{1}{2}}, θ^{n + \frac{1}{2}})

(18)

Then summing up (17) with (18) results in

\begin{aligned} {‖ θ^{n + 1} ‖}^{2} + {| q^{n + 1} |}^{2} - {‖ θ^{n} ‖}^{2} - {| q^{n} |}^{2} + \frac{2}{T} Δ t {| q^{n + 1} |}^{2} \\ = - 2 κ Δ t {‖ (- Δ)^{\frac{α}{2}} θ^{n + \frac{1}{2}} ‖}^{2} - 2 ϵ Δ t ‖ S^{\frac{β}{2}} θ^{n + \frac{1}{2}} ‖^{2} \end{aligned}

(19)

and dropping the uninfluential positive terms, we obtain (16).

The numerical procedure for the second-order scheme (15) is exactly the same for the first-order scheme (6). We denote $λ^{n + \frac{1}{2}} = e^{\frac{t^{n + \frac{1}{2}}}{T}} q^{n + \frac{1}{2}}$ and let

u^{n + \frac{1}{2}} = u_{1}^{n + \frac{1}{2}} + λ^{n + \frac{1}{2}} u_{2}^{n + \frac{1}{2}}

(20)

Since

θ^{n + \frac{1}{2}} : = \frac{1}{2} (θ^{n + 1} + θ^{n})

and we can rewrite (15a) as

\frac{θ^{n + \frac{1}{2}} - θ^{n}}{\frac{1}{2} Δ t} + k (- Δ)^{α} θ^{n + \frac{1}{2}} . \nabla θ *,^{n + \frac{1}{2}} + ε S^{β} θ^{n + \frac{1}{2}} = 0.

(21)

Plugging (2.2), we find that $θ_{i}^{n + \frac{1}{2}} (i = 1, 2)$ satisfy

{\begin{matrix} 2 θ_{1}^{n + \frac{1}{2}} + Δ t κ (- Δ)^{α} θ_{1}^{n + \frac{1}{2}} + Δ t ϵ S^{β} θ_{1}^{n + \frac{1}{2}} = 2 θ^{n} \\ 2 θ_{2}^{n + \frac{1}{2}} + Δ t κ (- Δ)^{α} θ_{2}^{n + \frac{1}{2}} + Δ t ϵ S^{β} θ_{2}^{n + \frac{1}{2}} = - Δ t u^{n + \frac{1}{2}} \cdot \nabla θ^{*, n + \frac{1}{2}} \end{matrix}

(22)

Once we get

θ_{i}^{n + \frac{1}{2}} (i = 1, 2)

, we can determine explicitly

λ^{n + \frac{1}{2}}

from (15b) as following

\begin{aligned} (\frac{2 T + Δ t}{T Δ t} - e^{\frac{2 t^{n + 1}}{T}} (u^{n + \frac{1}{2}} \cdot \nabla θ^{*, n + \frac{1}{2}}, θ_{2}^{n + \frac{1}{2}})) e^{- \frac{t^{n + 1}}{T}} λ^{n + 1} \\ = e^{\frac{t^{n + 1}}{T}} (u^{n + \frac{1}{2}} \cdot \nabla θ^{*, n + \frac{1}{2}}, θ_{1}^{n + \frac{1}{2}}) + \frac{2 q^{n}}{Δ t} \end{aligned}

(23)

Obviously, the implementation of the scheme should be accompanied by suitable

(θ^{0}, θ^{1})

. We can compute

θ^{1}

by employing the first step of the scheme (6).

Remark 2

We have shown that the first-order and second-order schemes based on the SAV approach satisfy an unconditionally energy stability. Furthermore, the total computational cost of the first-order and second-order schemes is equal to solving two linear equations with constant coefficients in each time step, which can be implemented in a very efficient way.

Numerical results

In this section, we present some computational experiments using the numerical schemes constructed in the last section to demonstrate their accuracy and efficiency. In the first example, we verify the convergence rates of the first-order and second-order schemes in time. We then test the accuracy and efficiency for the inviscid SQG equation (i.e. $κ = 0$ ). Finally, a comparative study in accuracy of the FSVV approach and the SAV approach is demonstrated to show the accuracy and efficiency. In all examples, we consider the periodic boundary condition and use a Legendre spectral element method (SEM) in space.

Accuracy test

In order to test the convergence rates in time, we consider the inviscid SQG equation (i.e. $κ = 0$ )

θ_{t} + u \cdot \nabla θ = 0

(24)

in the square domain

Ω = (0, 2 π)^{2}

with the constant velocity

u = (U_{1}, U_{2})

. The exact solution of the above equation is given by

θ (x, y, t) = θ_{0} (x - U_{1} t, y - U_{2} t)

The velocity is given as

U_{1} = 1

U_{2} = 1

and the initial condition is given as follows

θ (x, 0) = s i n (x) s i n (y) + c o s (y)

(25)

The parameters are set as

E L = 400

β = 0.9

ϵ = 1 / N

, and

m_{N} = \frac{2 N}{3}

. We used

141^{2}

Legendre nodes in the whole domain where the polynomial degree

N = 7

in each element for space discretization so that the spatial discretization error is negligible compared with the time discretization error.

In Figure 1, we plot the convergence rate of the $L^{2}$ error at $T = 1$ for the SQG equation by using first-order and second-order schemes. We observe that, the first-order SAV scheme loses accuracy with larger time steps, it can present good approximations and first-order accuracy only with time step is smaller than $Δ t \leq 10^{- 2}$ .

Figure 1.

Convergence test for the inviscid SQG equation using the SAV-BDF and SAV-CN schemes. (a) and (b) Plot $L^{2}$ errors of $θ$ as a function of $Δ t$ . (a) SAV-BDF error of $θ$ . (b) SAV-CN error of $θ$ . SAV-BDF: scalar auxiliary variable-backward differentiation formula; SAV-CN: scalar auxiliary variable-Crank–Nicolson.

The inviscid SQG equation

Next, we consider another inviscid SQG equation with streamwise velocity in the square domain $Ω = [0, 2 π]^{2}$ . The initial condition is given as follows

θ (x, 0) = e x p (- (x - π)^{2} - 16 (y - π)^{2})

(26)

We set

N = 7

E L = 2500

β = 0.9

κ = 0

, and we discretize the space using

351^{2}

nodes. The other parameters are chosen to be the same as the first test.

In Figure 2, we present the evolution of the numerical solutions from the SAV-BDF, SAV-CN schemes, and fSVV method,²⁶ respectively. The time step is set $Δ t = 10^{- 3}$ for all schemes. The numerical results show that the cyclone develops thin arms during the rotated process. The solutions spread outwards and become unstable because of their strong shear. Meanwhile, smaller-scale vortices can be generated from a mesoscale vortex. These visualizations are qualitatively consistent with results presented by Held et al.,⁵ Bonito and Nazarov¹⁶ and Song and Karniadakis.²⁶ In Figure 3, we present evolution of parameter $λ$ by using SAV-BDF and SAV-CN schemes with $Δ t = 10^{- 3}$ . It can be observed that the difference of $λ$ from $1$ is small, this means the simulation $θ$ tends to be accurate.

Figure 2.

The dynamic evolution of vortex growth driven by the SQG equation by using SAV-BDF, SAV-CN schemes, and fSVV method. Contours of the numerical solutions $θ$ at $T = 0, 8, 16, 22, 28$ . (a) SAV-BDF, (b) SAV-CN, and (c) Reference of FSVV. SQG: surface quasi-geostrophic; SAV-BDF: scalar auxiliary variable-backward differentiation formula; SAV-CN: scalar auxiliary variable-Crank–Nicolson; fSVV: fractional spectral vanishing viscosity.

Figure 3.

Numerical result of $λ$ .

The dissipative SQG equation

In this subsection, we consider the dissipative SQG equation which the solution corresponding to the initial condition given as

θ (x, 0) = s i n (x) s i n (y) + c o s (y)

(27)

We solve the system (5) and (2) based on the SAV approach by using the second-order SAV-CN scheme (17) with

α = 0.4

κ = 0.001

Ω = [0, 2 π]^{2}

. We set

N = 7

E L = 2500

Δ t = 10^{- 3}

ϵ_{N} = 1 / N

m_{N} = \frac{2 N}{3}

and

β = 0.9

. We discretize the space by the Legendre spectral method with

351 \times 351

modes in this test.

In Figure 4, we plot the evolution of the contours of $θ$ at various times by using the SAV-CN scheme and fSVV method, no visible difference is observed. Meanwhile, we found two schemes generate same sharp front along the hyperbolic saddle. All these numerical results are qualitatively consistent with the simulations shown by Constantin et al.,¹³ Bonito and Nazarov,¹⁶ and Song and Karniadakis.²⁶

Figure 4.

Evolution comparison of dissipative SQG equation by using the SAV-CN and fSVV schemes at $t = 1, 8, 16, 20$ . (a) SAV-CN. (b) Reference of fSVV.SQG: surface quasi-geostrophic; SAV-CN: scalar auxiliary variable-Crank–Nicolson; fSVV: fractional spectral vanishing viscosity.

Summary

We have presented an algorithm for approximating the 2D SQG equation based on the SAV approach. By introducing a scalar equation as a variable, we reformulated this equation into an equivalent system. Further, we constructed first-order and second-order schemes for the reformulation system. We carried out the stability analysis for these schemes. We show that the methods are unconditionally energy stable. Such that, the schemes allow an efficiently implemented for solving the SQG equation system. In each time step, the algorithm only requires solving two linear equations with constant coefficient, and there is no need to solve nonlinear algebraic equation when we compute the parameter $λ$ . Finally, we demonstrated our numerical results for inviscid and dissipative SQG equation by computing the numerical solution and comparing with the analytical solutions or the numerical solutions from fSVV method wherever possible. We provided several numerical results to demonstrate the robustness and accuracy of the SAV-BDF and SAV-CN schemes for the SQG equation. The method can be extended to the scheme of adaptive time steps by using an efficient time scheme and turbulence simulation.

Supplemental Material

sj-cls-1-act-10.1177_17483026231176203 - Supplemental material for Scalar auxiliary variable approache for the surface quasi-geostrophic equation

Supplemental material, sj-cls-1-act-10.1177_17483026231176203 for Scalar auxiliary variable approache for the surface quasi-geostrophic equation by Shengtao Shi and Fangying Song in Journal of Algorithms & Computational Technology

Supplemental Material

sj-ins-2-act-10.1177_17483026231176203 - Supplemental material for Scalar auxiliary variable approache for the surface quasi-geostrophic equation

Supplemental material, sj-ins-2-act-10.1177_17483026231176203 for Scalar auxiliary variable approache for the surface quasi-geostrophic equation by Shengtao Shi and Fangying Song in Journal of Algorithms & Computational Technology

Supplemental Material

sj-dtx-3-act-10.1177_17483026231176203 - Supplemental material for Scalar auxiliary variable approache for the surface quasi-geostrophic equation

Supplemental material, sj-dtx-3-act-10.1177_17483026231176203 for Scalar auxiliary variable approache for the surface quasi-geostrophic equation by Shengtao Shi and Fangying Song in Journal of Algorithms & Computational Technology

Supplemental Material

sj-bst-4-act-10.1177_17483026231176203 - Supplemental material for Scalar auxiliary variable approache for the surface quasi-geostrophic equation

Supplemental material, sj-bst-4-act-10.1177_17483026231176203 for Scalar auxiliary variable approache for the surface quasi-geostrophic equation by Shengtao Shi and Fangying Song in Journal of Algorithms & Computational Technology

Supplemental Material

sj-cls-5-act-10.1177_17483026231176203 - Supplemental material for Scalar auxiliary variable approache for the surface quasi-geostrophic equation

Supplemental material, sj-cls-5-act-10.1177_17483026231176203 for Scalar auxiliary variable approache for the surface quasi-geostrophic equation by Shengtao Shi and Fangying Song in Journal of Algorithms & Computational Technology

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship and/or publication of this article: This work was supported by the NSF of China (Grant No. 11901100) and the Natural Science Foundation of Fujian Province of China (Grant No. 2020J05111).

ORCID iDs

Shengtao Shi

Fangying Song

Supplemental material

Supplemental material for this article is available online.

Appendix

Artificial term

In the following reference,²⁷ authors have developed a numerical method for computing fractional Laplacians on complex-geometry domains, by considering the following eigenvalue problem (EVP) for the Laplacian: .0Δ028 (A.1)

- Δ u - λ u = 0, x \in Ω

(A.2)

proper boundary conditions .

The spectral element method (SEM)^28,29 is used for solving equations (A.1)–(A.2). Then, equations (A.1)–(A.2) can be written in the discretized form (A.3)

A_{N} U - λ M_{N} U = 0

where

N

represents the number of the degrees-of-freedom (DoF) of the linear system (A.3) for the given number of elements

E l

(see Figure A1) and polynomial degree N in each element.

A_{N}

is the corresponding matrix of the Laplacian operator under certain boundary conditions,

M_{N}

is the mass matrix, and

U

is the numerical solution of

u

. The continuous EVP is approximated by the numerical solution of the eigenpairs

(λ_{i}, ϕ_{i})_{i = 1}^{N}

of the matrix

K = M_{N}^{- 1} A_{N}

, and

λ_{1} \leq λ_{2} \leq λ_{3} \leq \dots \leq λ_{N}

Using the numerical eigenpairs $(λ_{i}, ϕ_{i})$ of the Laplace operator $- Δ$ from the SEM solution. Then the fractional operators $(- Δ_{N})^{α}$ and $S_{N}^{β}$ are defined as follows (A.4)

(- Δ_{N})^{α} u = \sum_{i = 1}^{N} λ_{i}^{α} u_{i} ϕ_{i}, S_{N}^{β} u = \sum_{i = 1}^{N} {\tilde{λ}}_{i}^{β} u_{i} ϕ_{i}

where

u_{i} = (u, ϕ_{i})_{N}

and (A.5)

{\tilde{λ}}_{i} = {\begin{matrix} 0, i \leq m_{N}, \\ \exp (- (\frac{d - i}{m_{N} - i})^{2}) λ_{i}, i > m_{N} \end{matrix}

Here

m_{N}

can have different forms

m_{N} = {\sqrt{N}, \frac{N}{2}, or \frac{2 N}{3} etc.}

.^30,31 Of course, the usual spectral approximations of equations (7) and (17) are recovered when

ϵ_{N} = 0

m_{N} = N

References

Charney

. On the scale of atmospheric motions. Boston: Springer, 1990. pp. 251–265.

Gill

. Atmosphere-ocean dynamics. 30. San Diego: Academic press, 1982.

Majda

. Introduction to PDEs and Waves for the Atmosphere and Ocean. 9. New York: American Mathematical Soc, 2003.

Constantin

Majda

Tabak

. Formation of strong fronts in the 2-D quasigeostrophic thermal active scalar. Nonlinearity 1994; 7(6): 1495.

Held

Pierrehumbert

Garner

et al. Surface quasi-geostrophic dynamics. J Fluid Mech 1995; 282: 1–20.

Constantin

. Behavior of solutions of 2d quasi-geostrophic equations. SIAM J Math Anal 1999; 30(5): 937–948.

Pedlosky

. Geophysical fluid dynamics. 710. New York: Springer, 1987.

Blumen

. Uniform potential vorticity flow: Part i. Theory of wave interactions and two-dimensional turbulence. J Atmos Sci 1978; 35(5): 774–783.

Pierrehumbert

Held

Swanson

. Spectra of local and nonlocal two-dimensional turbulence. Chaos, Solitons Fract 1994; 4(6): 1111–1116.

10.

Constantin

. Geometric statistics in turbulence. SIAM Rev 1994; 36(1): 73–98.

11.

Constantin

Cordoba

. On the critical dissipative quasi-geostrophic equation. Indiana Univ Math J 2001; 50(1): 97–107.

12.

Hmidi

Keraani

. Global solutions of the super-critical 2d quasi-geostrophic equation in Besov spaces. Adv Math (N Y) 2007; 214(2): 618–638.

13.

Constantin

Lai

M-C

Sharma

et al. New numerical results for the surface quasi-geostrophic equation. J Sci Comput 2012; 50(1): 1–28.

14.

Ohkitani

Yamada

. Inviscid and inviscid-limit behavior of a surface quasigeostrophic flow. Phys Fluids 1997; 9(4): 876–882.

15.

Meerschaert

Tadjeran

. Finite difference approximations for fractional advection–dispersion flow equations. J Comput Appl Math 2004; 172(1): 65–77.

16.

Bonito

Nazarov

. Numerical simulations of surface quasi-geostrophic flows on periodic domains. SIAM J Sci Comput 2021; 43(2): B405–B430.

17.

Shen

Yang

. The scalar auxiliary variable (SAV) approach for gradient flows. J Comput Phys 2018; 353: 407–416.

18.

Shen

Zhang

. A generalized sav approach with relaxation for dissipative systems. J Comput Phys 2022; 464: 111311.

19.

Shen

Yang

. A new class of efficient and robust energy stable schemes for gradient flows. SIAM Rev 2019; 61(3): 474–506.

20.

Lin

Yang

Dong

. Numerical approximation of incompressible Navier-Stokes equations based on an auxiliary energy variable. J Comput Phys 2019; 388: 1–22.

21.

Shen

Liu

. New SAV-pressure correction methods for the Navier-Stokes equations: Stability and error analysis. Math Comput 2022; 91(333): 141–167.

22.

Akrivis

. Energy-decaying extrapolated RK–SAV methods for the Allen–Cahn and Cahn–Hilliard equations. SIAM J Sci Comput 2019; 41(6): A3703–A3727.

23.

Cheng

Liu

Shen

. Generalized SAV approaches for gradient systems. J Comput Appl Math 2021; 394: 113532.

24.

Jiang

Zhang

Zhao

. Improving the accuracy and consistency of the scalar auxiliary variable (SAV) method with relaxation. J Comput Phys 2022; 456: 110954.

25.

Cheng

Liu

Shen

. A new Lagrange multiplier approach for gradient flows. Comput Methods Appl Mech Eng 2020; 367: 113070.

26.

Song

Karniadakis

. Fractional spectral vanishing viscosity method: application to the quasi-geostrophic equation. Chaos, Solitons Fract 2017; 102: 327–332.

27.

Song

Karniadakis

. Computing fractional Laplacians on complex-geometry domains: Algorithms and simulations. SIAM J Sci Comput 2017; 39(4): A1320–A1344.

28.

Patera

. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. J Comput Phys 1984; 54(3): 468–488.

29.

Karniadakis

. Spectral element simulations of laminar and turbulent flows in complex geometries. Appl Numer Math 1989; 6(1): 85–105.

30.

Karamanos

Karniadakis

. A spectral vanishing viscosity method for large-eddy simulations. J Comput Phys 2000; 163(1): 22–50.

31.

Maday

Kaber

SMO

Tadmor

. Legendre pseudospectral viscosity method for nonlinear conservation laws. SIAM J Numer Anal 1993; 30(2): 321–342.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.03 MB

0.02 MB

0.05 MB

0.03 MB

0.00 MB