 # Analytical travelling vortex solutions of hyperbolic equations for validating very high order schemes

Testing the order of accuracy of (very) high order methods for shallow water (and Euler) equations is a delicate operation and the test cases are the crucial starting point of this operation. We provide a short derivation of vortex-like analytical solutions in 2 dimensions for the shallow water equations (and, hence, Euler equations) that can be used to test the order of accuracy of numerical methods. These solutions have different smoothness in their derivatives (up to 𝒞^∞) and can be used accordingly to the order of accuracy of the scheme to test.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Moving vortex solutions and regularity requirements

### 1.1 Shallow water equations and moving profiles

We consider the shallow water equations (SWEs) on a flat bathymetry reading

 {∂th+∇⋅(h→u)=0,∂t(h→u)+∇⋅(h→u⊗→u)+gh∇h=0, (1)

, , .

Following  we consider solutions of the form and , with , and constant. Replacing in the SWEs we obtain

 ⎧⎨⎩→U0⋅∇→ζH0+H0∇→ζ⋅→U0=0,→U0⋅∇→ζ→U0+∇→ζH0=0, (2)

which justifies looking for stationary solutions with solenoidal velocity fields with depth variations uniquely in the cross-stream direction. These conditions are easily met for solutions with cylindrical symmetry.

### 1.2 SWEs in cylindrical coordinates

We consider cylindrical coordinates defined in 2D by the distance from the origin , and the counter-clockwise angle measured from the positive axis so that

 x=rcosθ,y=rsinθ. (3)

We also introduce the direction vectors

, and . This notation can be used to write the shallow water equations in polar coordinates as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂th+∂r(hur)+1r(∂θ(huθ)+hur)=0∂t(hur)+∂r(hu2r)+gh∂rh+1r(∂θ(huruθ)+h(u2r−u2θ))=0∂t(huθ)+∂r(huruθ)+1r(gh∂θh+∂θ(hu2θ)+2huruθ)=0 (4)

where

 ur=→u⋅^r=cos(θ)ux+sin(θ)uy,uθ=→u⋅^r⊥=−sin(θ)ux+cos(θ)uy,ux=cos(θ)ur−sin(θ)uθ,uy=sin(θ)ur+cos(θ)uθ. (5)

### 1.3 Stationary vortex ODE

To mimic (2), we consider the particular case with all zero time derivatives and only radial velocity, i.e.,

 ∂th=0,∂tur=∂tuθ=0,h=h(r),∂θh=0,uθ=uθ(r),∂θuθ=0,ur=0. (6)

With the above hypotheses we can readily check that the first and the last in (4) are identically satisfied, while the second reduces to

 h′(r)=u2θgr. (7)

Assuming further that we end up with

 h′(r)=rω2(r)g. (8)

Given a law for the angular velocity this ODE can be integrated to obtain closed form expressions for the depth and, conversely, given a law for , and can be obtained differentiating .

### 1.4 Regularity requirements for validating high order methods

To validate higher order methods, exact solutions of the shallow water system should have enough regularity to allow the validity of high order approximation results. The classical interpolation estimate for finite element approximations

[3, §1.5] for a function with , for a multi-index , is the following

 ∥v−vh∥Lp(Ω)≤chl+1|v|l+1,Lp(Ω), (9)

with the mesh size. This means that to benchmark an order accurate method, we need an exact solution with integrable derivatives. In practice, for a system different variables may have different regularity. This may lead in practice to convergence rates somewhat in between those of the different variables, depending on how the error is defined and measured.

### 1.5 Extension to Euler equations

Similarly to SWEs other moving vortexes can solve exactly Euler equations. They read

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩∂tρ+∂x(ρux)+∂y(ρuy)=0,∂t(ρux)+∂x(ρu2x+p)+∂y(ρuxuy)=0,∂t(ρuy)+∂x(ρuxuy)+∂y(ρu2y+p)=0,∂t(ρE)+∂x[ux(ρE+p)]+∂y[uy(ρE+p)]=0, (10)

where the perfect gas equation of state closes the system, i.e.,

 ρE=pγ−1+12ρ(u2x+u2y). (11)

Here is the adiabatic constant. Equivalently, we can write Euler equations in polar coordinates

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩∂tρ+1r∂r(rρur)+1r∂θ(ρuθ)=0,∂t(ρur)+1r∂r(r(ρu2r+p))+1r∂θ(ρuruθ)−1r(ρu2θ+p)=0,∂t(ρuθ)+1r∂r(rρuruθ)+1r∂θ(ρu2θ+p)+1rρuθur=0,∂t(ρE)+1r∂r[rur(ρE+p)]+1r∂θ[uθ(ρE+p)]=0. (12)

Using the vortex form (6), we have that is set to 0 and is equal to 0 for all the unknowns. Hence, the system reduces for steady vortexes to

 r∂rp=ρu2θ. (13)

This form can be easily solved in two situations.

#### Isentropic case

The isentropic case (constant) for leads to

 γγ−1∂r(ργ−1)=u2θr. (14)

Following the approach of (7), if we define any angular velocity and consequentially

 ur=0,uθ(r):=rω(r),ρ(r):=(ρ0+∫γ−1γrω2(r)dr)1γ−1,p=ργ, (15)

where with the integral is meant up to a constant, that can be set with at . Clearly, this formulation is equivalent to the SWE ones setting and , and all the following derivation can be used in Euler equations.

#### Isochoric vortex

The second option to obtain a steady vortex is to set a constant density , and (13) becomes a differential equation for . Again, we obtain

 ur=0,uθ(r):=rω(r),ρ(r)=ρ0,p(r):=∫ρ0rω2(r)dr. (16)

## 2 The vortex solution of 

The following is an example often used in literature obtained by setting

 ω(r)={Γ(1+cos(πrr0))if r≤r0,0otherwise. (17)

This expression can be fed into (8) and integrated backwards from to with initial condition to obtain

 h=h0−Γ2r20gπ2{H(π)−H(πrr0)if r≤r0,0otherwise, (18)

with

 H(x)=2cosx+2xsinx+cos(2x)8+xsin(2x)4+12x216. (19)

Note that

 H′(π)=0,u′θ(π)=0,H(2)(π)=0,u(2)θ(π)=π≠0,H(3)(π)=0,u(3)θ(π)≠0,H(4)(π)=0,u(4)θ(π)≠0,H(5)(x)=6π≠0,u(5)θ(π)≠0. (20)

According to (9) this solution should allow to obtain at most second order of accuracy, due to the limited regularity of the velocity. In practice, some computations have shown a little extra convergence for the depth, for which, at least on relatively coarse resolutions, one can manage to obtain toughly third order convergence. The convergence study obtained with a WENO5 on Cartesian grids on with a RK(6,5) is shown in table 1. The tests are run on the vortex defined with and centered in , and such that . It is clear that the order of accuracy for cannot reach 3 and that the velocities can be approximated with only order 2.

## 3 Non compact supported vortexes

There are many other vortexes for Euler’s equation, Euler himself in  suggested a similar simplification to obtain a steady solution, and in 1998 Shu  proposed a vortex to study the accuracy of WENO5 schemes from a qualitatively point of view. In the following years, this tests and some modifications of it became a real benchmark for testing the accuracy of many schemes [5, 10, 9, 1, 2]. The main idea follows form the definition of a vortex where the base function is a Gaussian function , where is a rescaling factor for the width and rescales the amplitude. Despite being the benchmark test for many problems, this test does not have a compact support and boundaries are a real issue, in particular when dealing with very high order of accuracy methods, as shown, for instance, in . Nevertheless, it was often used with periodic boundary conditions even with nonzero background speed. This leads to traveling discontinuities (on the derivatives) from the boundaries to all over the domain. The techniques to cure this issue were various starting from enlarging the domain (or equivalently reducing ) to treating with nonperiodic boundary conditions (inflow/outflow with steady vortex). Of course increasing the domain means that the computational costs increase too and using inflow/outflow does not allow to make the vortex travel outside the domain.

We can easily compute

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ω(r)=Γe−(r/r0)2uθ(r)=rω(r)h(r)=h0−Γ2r204ge−2(r/r0)2=h0−r204gω2(r). (21)

We plot the profile of and for , such that and different in figs. 2 and 1.

On the pictures one finds also the values of the maximum absolute values of the derivatives and the values of speed and height at . As one can see, increasing the derivatives decrease, but the values at are diverging from for and 0 for . This means that boundary effects might disturb the convergence process. For the numerical part, we will test only , as all the other parameters leads to either very high derivatives or boundary effects. This is also visible in these tests, but this is the range where they are more bounded. Figure 1: h profile and its derivatives for vortexes (21) with r0=0.2,0.3,0.4,0.6,0.8,1 Figure 2: uθ profile and its derivatives for vortexes (21) with r0=0.2,0.3,0.4,0.6,0.8,1

## 4 Compact supported traveling vortexes of arbitrary smoothness

### 4.1 Iterative correction of the RB-vortex to obtain arbitrary smoothness

In order to obtain more regularity in the vortex, and to be able to test higher order accuracy methods, we can generalize the previous solution. We remark that for (17) is equivalent to

 ω=Γ(1+cos(ρ))=2Γcos2(ρ/2),ρ=πrr0. (22)

A natural way to improve the regularity of this definition for is to increase the exponent of the cosinus. So we look into definitions of the type

 ω=2pΓcos2p(ρ/2)=Γ(1+cosρ)p,p≥1 (23)

which allows to increase the regularity, keeping bounded values of the th derivative.

In practice, we need to integrate the ODE

 h′(r)=4pΓ2grcos4p(πr2r0),r∈[0,r0], (24)

with the condition . The solution to this problem within can be written as

 h(r)=h0−1g(2pΓr0π)2(Hp(π/2)−Hp(ρ/2)),ρ=πrr0, (25)

having set

 Hp(x)=∫xycos4p(y)dy (26)

up to an additional constant. Comparing to (18) we can deduce that

 H1(x)=116H(2x) (27)

with defined in (19). For we can use the iterative integration rule

 ∫cosn(x)dx=cosn−1(x)sin(x)n+n−1n∫cosn−2(x)dx (28)

to show that

 Hp(x):=∫xcos4p(x)=x∫cos4p(x)−∫∫cos4p(x)=xcos4p−1(x)sin(x)4p+x4p−14p∫cos4p−2(x)−∫{cos4p−1(x)sin(x)4p+4p−14p∫cos4p−2(x)}=xcos4p−1(x)sin(x)4p−∫cos4p−1(x)sin(x)4p+4p−14p{x∫cos4p−2(x)−∫∫cos4p−2(x)}=xcos4p−1(x)sin(x)4p+cos4p(x)(4p)2+4p−14p{xcos4p−3(x)sin(x)4p−2+cos4p−2(x)+4p−34p−2∫xcos4p−4(x)}=4p−14p4p−34p−2Hp−1(x)+xcos4p−3(x)sin(x)4p{4p−14p−2+cos2(x)}+cos4p−2(x){cos2(x)(4p)2+4p−14p(4p−2)2}.

This formula can be used to compute the exact solution iteratively for any with initial value given by (27). In appendix A, we provide some values of . Figure 3: h profile and its derivatives for vortexes (25) with p=1,…,5 Figure 4: uθ profile and its derivatives for vortexes (25) with p=1,…,5

These vortexes suited for testing high order methods, as they are . Nevertheless, it is not recommended to use a very high for not so high order methods. Indeed, we see that the estimation on the error (9) depends not only on the smoothness of the solution, but also on the seminorms. Essentially, if the derivatives are too large, we might need very fine meshes to reach the expected order of accuracy. To assess this information we plot for various vortexes the profile of and of and their derivatives (up to the 5th derivative). We fix for all the vortexes , the values and , by setting and . In figs. 4 and 3 the plot related to the vortexes of type (25) for variables and respectively.

In fig. 3 we barely see that the fifth derivative of the case is not zero at , but we observe that the amplitude of the derivative functions increases with , in particular the fifth derivative of is 10 times larger than . In fig. 4 we immediately see that the second derivative of for is discontinuous in , and the fourth derivative is discontinuous also for . Again, the higher we choose the larger the derivatives become.

### 4.2 A few C∞ examples

In this case we start from (8) and we reverse it to obtain a definition of the angular velocity given the depth:

 ω=√gh′r. (29)

To avoid the singularity in , we define the depth as a function of . An example of compactly supported function is obtained with the definition

 h=h0−Γ2⎧⎪⎨⎪⎩e−1(1−ρ)2 if ρ<1,0 else,ρ=(rr0)2. (30)

The angular velocity, and thus the tangential linear velocity, can be readily computed from (29):

 ω=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩2Γe−12(1−ρ)2√gr0(1−ρ)3 if r

Though being , the previous function has very large derivatives and need very fine mesh to see the order of the mesh.

Alternatives could be obtained, for instance, using higher power in the exponential, for example

 h=h0−Γ2⎧⎪⎨⎪⎩e−1(1−ρ)p if ρ<1,0 else,ρ=(rr0)2, (32)

 ω=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩Γ√2pgr0(1−ρ)p+1e−12(1−ρ)p if r

Or one can make the derivatives a bit smaller with an additional function, i.e.,

 h=h0−Γ2⎧⎪⎨⎪⎩e−1arctanp(1−ρ) if ρ<1,0 else,ρ=(rr0)2, (34)

 ω=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩Γe−12arctan(1−ρ)p√2pgr0arctanp+1(1−ρ)11+(1−ρ)2 if % r Figure 5: h profile and its derivatives for vortexes (32) with p=2,…,4 Figure 6: uθ profile and its derivatives for vortexes (32) with p=2,…,4 Figure 7: h profile and its derivatives for vortexes (34) with p=2,…,4 Figure 8: uθ profile and its derivatives for vortexes (34) with p=2,…,4

As for the previous examples, though being all these vortexes suited for testing high order methods, as they are , the amplitude of their derivatives changes considerably. To assess this information we plot for various vortexes the profile of and of and their derivatives (up to the 5th derivative). We fix for all the vortexes , the values and , by setting and . In figs. 6 and 5 we plot the vortexes of type (32) for variables and respectively, while in figs. 8 and 7 we plot the vortexes for (34).

With the vortexes (32) and (34), the profile of grows to more slowly, see fig. 5, but the high derivatives are larger with respect to the ones of the previous examples, but they do not vary much with respect to the parameter . Also for we observe similar behaviors: comparing the 5th derivatives of these vortexes with the previous ones, we have values at least 100 times larger. With the test cases we actually can slightly decrease the amplitude of high derivatives for large , as shown in fig. 8.

## 5 Numerical tests

We test a 5th order method with the presented vortexes. The method consists of a finite volume discretization on a Cartesian grid with WENO5 reconstruction and Rusanov numerical flux. We use 4 Gauss-Legendre points for one-dimensional quadrature rules. The time discretization is carried out with Butcher’s RK(6,5) 6 stages 5th order method, see appendix B for its Butcher’s tableau. The CFL number used is 0.95. The domain is and it is discretized with uniform Cartesian grids and periodic boundary conditions. The vortexes are set with , and , while for the non compact supported one we test different parameters. Final time is set to , but it might be useful to run the code for larger or smaller times. In  it is well observed how this impact on the convergence of the method as the machine precision error accumulates exponentially in time, one might see its effect before reaching the desired accuracy. On the other side, small final time leads to very small errors also for coarser meshes, not letting appreciate the quality of the high order schemes.

With the non compact supported vortex we obtain the results in table 2. We observe for large values of that the error from the boundary effect is predominant already for very coarse meshes. For (corresponding to in figs. 2 and 1) we have a discontinuity of the order of for the velocity variables on the boundary of the domain (periodic BCs) and this prevent reaching more than 3rd order in the coarse regime. For we see again boundary effects but for finer meshes, hence it is possible to reach an order of accuracy 4 for few steps in the mesh refinement process. On the other side, we have for (corresponding to in figs. 2 and 1) that at the boundaries is up to machine precision and . Indeed, the convergence in this test is smooth enough for , even for not too fine meshes. Nevertheless, one should really be careful using this vortex in hitting machine precision at the boundaries.

For the tests run with the vortexes (25), we observe that for we reach order 3 for and around for , a little better than expected, for the order reaches 4 and almost 5 for , and finally for we essentially obtain the expected 5th order of the scheme already at the mesh .

For the cases we observe that a finer mesh is needed to reach the expected accuracy, in particular for the (32) tests, where we barely reach it for with the mesh with . We observe in this case that we obtain better convergence orders for smaller times, e.g. .

In the vortex (34) we have smaller amplitude of higher derivatives, in particular for large , so with we see an order of accuracy larger than 4 for , which is a bit better than with the only exponential test, still with a quite fine mesh. The same reasoning holds for these tests and with smaller final times we can observe better convergence order (around 4.8 for ).

## 6 Conclusion

Steady and moving vortex solutions are useful tools to test the order of accuracy of high order methods for hyperbolic balance laws (shallow water and Euler equations). The vortexes known in literature are, anyway, either discontinuous in some derivatives , hence not suitable for very high order methods, or non compact supported , hence presenting troubles with boundary conditions if the parameters are not carefully chosen.

In this work we propose a class of compactly supported vortexes with arbitrary , and some compactly supported vortexes. The latter are very attractive as they do not need to compute integrals, but only derivatives, and have all the derivatives continuous. On the other side, their derivatives are particularly large and the expected order of accuracy is reached only for very fine meshes. The one that obtain better results in this class is the one based on the with for final time . For larger times the order is not so neat.

Among the based vortexes, the one with gives better results (for order 5 schemes) even for coarser meshes , as its derivative are not too large. On the other side, this vortex can only be used to test methods up to order 6 of accuracy.

Overall, we suggest the use of vortexes for testing arbitrarily high order methods, while, for a fixed method of order the recipe (25) builds a vortex where the order of accuracy can be observed in coarser meshes.

## Acknowledgments

D.T. acknowledges Wasilij Barsukow for a fruitful discussion on vortexes and their origin.

## Appendix A Values of Hp

Here, we provide the explicit form of for , the corresponding can be computed from (25). The value of can be set given , using

 Γ=π2pr0√g(h0−hmin)Hp(π/2)−Hp(0). (36)

We provide in the following also the value of .

 H1(x)= cos(2x)8+xsin(2x)4+cos(2x)264+3x216+xcos(2x)sin(2x)16, H1(π/2)− H1(0)=3π264−14,
 H2(x)= 35cos(2x)384+35xsin(2x)192+cos(x)6(cos(x)264+7288)+35cos(2x)23072+ 35x2256+35xcos(2x)sin(2x)768+xcos(x)5sin(x)(cos(x)2+7/6)8, H2(π/2)− H2(0)=35π21024−29,
 H3(x)= 77cos(2x)1024+77xsin(2x)512+33cos(x)6(cos(x)2/64+7/288)40+cos(x)10(cos(x)2144+111200)+ 77cos(2x)28192+231x22048+77xcos(2x)sin(2x)2048+33xcos(x)5sin(x)(cos(x)2+7/6)320+ xcos(x)9sin(x)(cos(x)2+11/10)12, H3(π/2)− H3(0)=231π28192−3591800.

## Appendix B Runge Kutta (6,5)

The Butcher Tableau of Butcher’s Runge Kutta (6,5) method used is the following

 0141414181812001234316−38389161−378767−12787790016452151645790. (37)

## References

•  R. Abgrall, P. Bacigaluppi, and S. Tokareva (2019) High-order residual distribution scheme for the time-dependent Euler equations of fluid dynamics. Computers & Mathematics with Applications 78 (2), pp. 274–297. Cited by: §3.
•  R. Abgrall and D. Torlo (2020) High order asymptotic preserving deferred correction implicit-explicit schemes for kinetic models. SIAM Journal on Scientific Computing 42 (3), pp. B816–B845. Cited by: §3.
•  A. Ern and J. Guermond (2013) Theory and practice of finite elements. Vol. 159, Springer Science & Business Media. Cited by: §1.4.
•  L. Euler (1755) Principes généraux du mouvement des fluides. Académie Royale des Sciences et des Belles Lettres de Berlin Mémoires (11), pp. 274–315. Cited by: §3.
•  J. S. Hesthaven and T. Warburton (2007) Nodal discontinuous galerkin methods: algorithms, analysis, and applications. Springer Science & Business Media. Cited by: §3.
•  M. Ricchiuto and A. Bollermann (2009) Stabilized residual distribution for shallow water simulations. Journal of Computational Physics 228 (4), pp. 1071–1115. Cited by: §1.1, §2, §6.
•  C. Shu (1998) Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In Advanced numerical approximation of nonlinear hyperbolic equations, pp. 325–432. Cited by: §3, §6.
•  S. C. Spiegel, H. Huynh, and J. R. DeBonis (2015) A survey of the isentropic Euler vortex problem using high-order methods. In 22nd AIAA Computational Fluid Dynamics Conference, pp. 2444. Cited by: §3, §5.
•  Z. J. Wang and H. Gao (2009) A unifying lifting collocation penalty formulation including the discontinuous galerkin, spectral volume/difference methods for conservation laws on mixed grids. Journal of Computational Physics 228 (21), pp. 8161–8186. Cited by: §3.
•  Z. J. Wang, Y. Liu, G. May, and A. Jameson (2007) Spectral difference method for unstructured grids ii: extension to the Euler equations. Journal of Scientific Computing 32 (1), pp. 45–71. Cited by: §3.