EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 12, No. 4, 2019, 1464-1482 ISSN 1307-5543 – www.ejpam.com Published by New York Business Global The higher-Order CESE Method for two-dimensional Shallow Water Magnetohydrodynamics Equations Sidrah Ahmed1,∗, Saqib Zia2 1 Mathematical Sciences, Sukkur IBA Uniersity, Sukkur, Pakistan 2 Department of Mathematics, COMSATS Unviversity Islamabad, Park Road Chak Shehzad Islamabad, Pakistan Abstract. The numerical solution of one and two dimensional shallow water magnetohydro- dynamics model is obtained using the 4th-order conservation element solution element method (CESE). The method is based on unified treatment of spatial and temporal dimensions contrary to the finite difference and finite volume methods. The higher-order CESE scheme is constructed us- ing same definitions of conservation and solution elements that are used for 2nd-order CESE scheme formulation. Hence it is more convenient to increase accuracy of CESE methods as compared to the finite difference and finite volume methods. Moreover the scheme is developed using the con- servative formulation and does not require change in the source term for treating the degenerate hyperbolic nature of shallow water magnetohydrodynamics system due to divergence constraint. The spatial and temporal derivatives have been obtained by incorporating 3rd-order Taylor expan- sion and the projection method is used to handle the divergence constraint. The accuracy and robustness of the extended method is tested by performing benchmark numerical tests taken from the literature. Numerical experiments revealed the accuracy and computational efficiency of the scheme. 2010 Mathematics Subject Classifications: 35L60, 35L65, 35R05 Key Words and Phrases: Shallow Wanter Equations, Hyperbolic Conservation Laws, Nonlinear Partial Differential Equations 1. Introduction Gilman was the pioneer to introduce the shallow water magnetohydrodynamics (here- after referred as SWMHD) system to mathematically model the physical phenomena in the solar tachocline [9]. The solar tachocline is a thin region inside a star where heat radi- ates from the interior zone to the outer turbulent zone predominantly. Its thickness ranges between 2 to 5 percentage of the solar radius. Because of predominantly radiative heat transfer and its thin vertical thickness the tachocline can be described similarly to ocean ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v12i4.3538 Email addresses: dr.sidrah@iba-suk.edu.pk (A. Sidrah), saqibzia81@hotmail.com (Z. Saqib) http://www.ejpam.com 1464 c© 2019 EJPAM All rights reserved. A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1465 and atmosphere. This urges the introduction of MHD analogue of shallow water system (SW) to study physical processes in solar tachocline. Although the SWMHD system has been derived with one layer assumption, it is possible to add little changes to link more than one layer as in hydrodynamics case. This indeed allows to study several kind of solar dynamics within solar tachocline reference their in [5–7, 9]. Following the derivations by [9] and Rossmanith [16], the SWMHD system is written below for incompressible flow under magnetohydrodtatic equilibrium in the z-direction. ρ ≡ constant , (1) ∂ ∂z ( p+ ‖B‖2 2 ) = ρg. (2) The reader is referred to [10] for complete derivation of the model. The one dimensional SWMHD equations in the conservative formulation is given as below: ∂q ∂t + ∂f ∂x + ∂g ∂y = 0, (3) q =  h hu1 hu2 hB1 hB2  t , f =  hu1 hu2 1 + 1 2gh 2 − hB2 1 hu1u2 − hB1B2 0 hu1B2 − hu2B1  , g =  hu2 hu1u2 − hB1B2 hu2 2 + 1 2gh 2 − hB2 2 hu2B1 − hu1B2 0  , (4) where g is the gravitational constant, the variables v1, v2, B1, B2 are the velocity and magnetic field components respectively and h is the layer depth. This system is degenerate hyperbolic due to the last equation ∇ · (hB) = 0. Hence in order to develop a Riemann solver for numerical solution of SWMHD model one must needs to relax the equation ∇ · (hB) = 0. This will need an addition of extra source term as can be seen in [22]. This is one of the reasons to present my work based on developing a numerical scheme without change in source term. The one and two dimensional SWMHD system has been solved numerically with dif- ferent schemes [12, 15–17]. Due to the eigenvalue structure of MHD analogue of shallow water system, it is obvious to design a robust and accurate numerical method particularly for one dimensional case to obtain a higher order extension that captures all waves region. This is the second reason to carry out the present work. The aim of this paper is to build a novel method for one dimensional SWMHD model with the following features: a. The method should be accurate for one dimensional case b. It is easy to extend the method to higher dimensions c. There is no need of Riemann solvers d. There is no need to change the source term The present work follows [23], in which authors have presented a fourth order CESE method for MHD system. CESE method is widely applied to solve Riemann problems associated to hyperbolic conservation laws. It was first proposed by Chang et al. [3]. A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1466 Now a large number of studies are available in the literature that demonstrate the CESE method’s robustness, easy implementation process and accuracy [4, 11–15, 23]. The ap- plications of CESE method for MHD problems can be found in [18, 19, 24, 26]. This paper is organized in four sections. Section 1 gives an introduction of SWMHD model and the numerical schemes that have been applied to this system. I have tried to emphasize on the causes that urges to develop a higher order CESE scheme for this model. Section 2 comprises of detailed formulation of fourth order CESE scheme for one dimensional SWMHD equations. Section 3 includes the numerical test cases that demonstrate the reliability of present work. In Section 5, concluding remarks are given about this work and some future extensions. 2. Fourth-Order CESE Scheme The fourth-order CESE scheme is based on the same basic concept of the second- order CESE scheme in [25]. Let us consider the space-time Cartesian coordinates system E3(x, y, t), where x and y are the spatial coordinates combined with the time coordinate t. Eq.3 can be rewritten as ∇ · hm = 0, (5) where hm = (fm, gm, qm), (m = 1, · · ·, 5) is the spacetime flux vector and ∇ = ( ∂ ∂x , ∂ ∂y , ∂ ∂t). Integrating equation refe:3 over the arbitrary control volume v and applying the Gausss divergence theorem in the 3D spacetime domain E3(x, y, t), allows to write the following form: ∮ A(v) hm · ds = 0, (6) where A(v) denotes the closed boundary surface of any spacetime region v in E3(x, y, t), da = (dδ)n, with the area dδ and the unit outward normal vector n. Moreover hm · ds denotes the total spacetime flux leaving the surface element da. Figure 1 shows the projection of conservation element on 2D Cartesian coordinates grid in (x, y) plane at time level n− 1 2 as the polygon V1V5V7V9. Each vertex Vi is marked as a filled circle. The centroid Ci, i = 1, 2, 3, 4, associated to each polygon is marked as hollow circle. The 3D space-time region is divided into non-overlapping quadrilaterals called conservation elements hereafter referred as CE. Each CE is divided into four basic CEs. The associated BCEs are V1V2V3V4, V3V4V5V6, V6V7V8V3 and V8V3V2V9. The cen- troid of each BCE is marked as Ci and the centroid of CE is marked as G. Due to uniform grid, the points G and V3 coincide. The staggered CESE scheme is designed in a way that the solution at BCE cell centers Ci are obtained during the first half time step and the solution at the vertices (in our case G) is updated in the second half time step. For each solution point G there is an A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1467 Figure 1: Projection of CE onto (x, y) plane associated solution element (SE). The role of solution element is to get Taylor approxima- tion of flux variables qm, fm and gm. Consider for example the solution point C / i . The corresponding CE is the polygon V1V2V3V4V / 1 V / 2 V / 3 V / 4 while the SE is the union of three planes V / 1 V / 2 V / 3 V / 4 , A1A3A // 1 A // 3 and A2A4A // 2 A // 4 . This space-time discretization implies that the boundaries of CEs are parallel to the coordinate planes and also their outward unit normal vectors are along the coordinate axis. It results only one component of total space-time flux vector (f, g, q) passing through boundary planes of CEs and BCEs. The conservation law form of shallow water MHD equations holds for an arbitrary closed spacetime domain and the defined CE is also closed space-time region so the conservation law can be transformed to discrete form. Continuing the example of solution point C / i and impose the conservation law 6 on its CE. Contrary to the finite difference schemes, the CESE method is extended to higher order schemes based on same grid and associated CEs and SEs which is one of the reason for convenient extension of CESE scheme to higher order as compared to other numerical methods. Equation 6 implies that ∮ S(CE(C1)) hm · ds = ∮ S(V1V2V3V4V ′1V ′2V ′3V ′4) hm · ds. (7) The solution at C ′i is updated by carrying out integration over the bottom surface V1V2V3V4, side surfaces V1V ′ 1V2V ′ 2 , V2V ′ 2V3V ′ 3 , V3V ′ 3V4V ′ 4 and V4V ′ 4V1V ′ 1 and the top surface V ′1V ′ 2V ′ 3V ′ 4 . The unit outward normal to the bottom surface is (0, 0,−1). Hence only space-time flux q passes through this surface. As seen in Figure 2, the bottom surface is divided into A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1468 Figure 2: Schematic diagram of CE and SE four parts. On each part qm is approximated by third-order Taylor expansion at points V1, V2, V3, V4. The integration over V1A2C1A1 is carried out using the third-order Taylor expansion of qm at the point V1. It gives following result:∫ i∆x (i− 1 2)∆x(j− 1 2)∆y { (qm)V1 + (qmx)V1 [ x− ( i− 1 2 ) ∆x ] + (qmy)V1 [ y − ( j − 1 2 ) ∆y ] + 1 2 (qmxχ)V1 [ x− ( i− 1 2 ) ∆x ]2 + 1 2 (qmyy)V1 [ y − ( j − 1 2 ) ∆y ]2 + 1 2 [ (qmxy)V1 + (qmyx)V1 ] [ x− ( i− 1 2 ) ∆x ] [ y − ( j − 1 2 ) ∆y ] + 1 6 (qmxxx)V1 [ x− ( i− 1 2 ) ∆x ]3 + 1 6 (qmyyy)V1 [ y − ( j − 1 2 ) ∆y ]3 + 1 6 [ (qmyyx)V1 + (qmyxy)V1 + (qmxyy)V1 ] [ y − ( j − 1 2 ) ∆y ]2 [ x− ( i− 1 2 ) ∆x ] + 1 6 [ (qmxxy)V1 + (umxyx)V1 + (qmyxx)V1 ] [ x− ( i− 1 2 ) ∆x ]2 [ y − ( j − 1 2 ) ∆y ] }dxdy = ∆x∆y 4 { (qm)V1 + (qmx)V1 ∆x 4 + (qmy)V1 ∆y 4 + (qmxx)V1 ∆x2 24 + (qmyy)V1 ∆y2 24 + [ (umxy)V1 + (qmyx)V1 ] ∆x∆y 32 + (qmxxx)V1 ∆x3 192 + (qmyyy)V1 ∆y3 192 + [ (qmyyx)V1 + (qmyxy)V1 + (qmxyy)V1 ] ∆y2∆x 288 + [ (qmxxy)V1 + (qmxyx)V1 + (qmyxx)V1 ] ∆x2∆y 288 }. (8) A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1469 The integration over the surface V2A2C1A3 results as follows:∫ (i+ 1 2)∆x i∆x ∫ j∆y (j− 1 2)∆y { (qm)V2 + (qmx)V2 [ x− ( i+ 1 2 ) ∆x ] + (qmy)V2 [ y − ( j − 1 2 ) ∆y ] + 1 2 (qmxχ)V2 [ x− ( i+ 1 2 ) ∆x ]2 + 1 2 (qmyy)V2 [ y − ( j − 1 2 ) ∆y ]2 + 1 2 [ (qmxy)V2 + (qmyx)V2 ] [ x− ( i+ 1 2 ) ∆x ] [ y − ( j − 1 2 ) ∆y ] + 1 6 (qmxχX)V2 [ x− ( i+ 1 2 ) ∆x ]3 + 1 6 (qmyyy)V2 [ y − ( j − 1 2 ) ∆y ]3 + 1 6 [ (qmyyx)V2 + (qmyxy)V2 + (qmxyy)V2 ] [ x− ( i+ 1 2 ) ∆x ] [ y − ( j − 1 2 ) ∆y ]2 + 1 6 [ (qmxxy)V2 + (qmxyx)V2 + (qmy×x)V2 ] [ x− ( i+ 1 2 ) ∆x ]2 [ y − ( j − 1 2 ) ∆y ] }dxdy = ∆x∆y 4 { (qm)V2 − (qmx)V2 ∆x 4 + (qmy)V2 ∆y 4 + (qmxx)V2 ∆x2 24 + (qmyy)V2 ∆y2 24 − [ (qmxy)V2 + (qmyx)V2 ] ∆x∆y 32 − (qmxxx)V2 ∆x3 192 + (qmyyy)V2 ∆y3 192 − [ (qmyyx)V2 + (qmyxy)V2 + (qmxyy)V2 ] ∆y2∆x 288 + [ (qmxxy)V2 + (qmxyx)V2 + (qmyxx)V2 ] ∆x2∆y 288 }. (9) The integration over polygon V3A3C1A4 gives:∫ (i+ 1 2)∆x(j+ 1 2)∆y i∆x { (qm)V3 + (qmx)V3 [ x− ( i+ 1 2 ) ∆x ] + (qmy)V3 [ y − ( j + 1 2 ) ∆y ] + 1 2 (qmxχ)V3 [ x− ( i+ 1 2 ) ∆x ]2 + 1 2 (qmyy)V3 [ y − ( j + 1 2 ) ∆y ]2 + 1 2 [ (qmxy)V3 + (qmyx)V3 ] [ x− ( i+ 1 2 ) ∆x ] [ y − ( j + 1 2 ) ∆y ] + 1 6 (qmxxx)V3 [ x− ( i+ 1 2 ) ∆x ]3 + 1 6 (qmyyy)V3 [ y − ( j + 1 2 ) ∆y ]3 + 1 6 (qmxxx)V3 [ x− ( i+ 1 2 ) ∆x ]3 + 1 6 (qmyyy)V3 [ y − ( j + 1 2 ) ∆y ]3 + 1 6 [ (qmyyx)V3 + (qmyxy)V3 + (qmxyy)V3 ] [ x− ( i+ 1 2 ) ∆x ] [ y − ( j + 1 2 ) ∆y ]2 + 1 6 [ (qmxxy)V3 + (qmxyx)V3 + (qmyxx)V3 ] [ x− ( i+ 1 2 ) ∆x ]2 [ y − ( j + 1 2 ) ∆y ] }dxdy A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1470 = ∆x∆y 4 { (qm)V3 − (qmx)V3 ∆x 4 − (qmy)V3 ∆y 4 + (qmxx)V3 1 24 ∆x2 + (qmyy)V3 1 24 ∆y2 + [ (qmxy)V3 + (qmyx)V3 ] ∆x∆y 32 − (qmxxx)V3 ∆x3 192 − (qmyyy)V3 ∆y3 192 − [ (qmyyx)V3 + (qmyxy)V3 + (qmxyy)V3 ] ∆y2∆x 288 − [ (qmxxy)V3 + (qmxyx)V3 + (qmyxx)V3 ] ∆x2∆y 288 }. (10) Finally, integration on the part V4A1C1A4 gives:∫ i∆x (i− 1 2)∆x ∫ (j+ 1 2)∆y j∆y { (qm)V4 + (qmx)V4 [ x− ( i− 1 2 ) ∆x ] + (qmy)V4 [ y − ( j + 1 2 ) ∆y ] + 1 2 (qmxx)V4 [ χ− ( i− 1 2 ) ∆x ]2 + 1 2 (qmyy)V4 [ y − ( j + 1 2 ) ∆y ]2 + 1 2 [ (qmxy)V4 + (qmyx)V4 ] [ x− ( i− 1 2 ) ∆x ] [ y − ( j + 1 2 ) ∆y ] + 1 6 (qmxxx)V4 [ x− ( i− 1 2 ) ∆x ]3 + 1 6 (qmyyy)V4 [ y − ( j + 1 2 ) ∆y ]3 + 1 6 [ (qmyyx)V4 + (qmyxy)V4 + (qmxyy)V4 ] [ x− ( i− 1 2 ) ∆x ] [ y − ( j + 1 2 ) ∆y ]2 + 1 6 [ (qmxxy)V4 + (qmxyx)V4 + (qmyxx)V4 ] [ x− ( i− 1 2 ) ∆x ]2 [ y − ( j + 1 2 ) ∆y ] }dxdy = ∆x∆y 4 { (qm)V4 + (qmx)V4 ∆x 4 − (qmy)V4 ∆y 4 + (qmxx)V4 ∆x2 24 + (qmyy)V4 ∆y2 24 − [ (qmxy)V4 + (qmyx)V4 ] ∆x∆y 32 + (qmxxx)V4 ∆x3 192 − (qmyyy)V4 ∆y3 192 + [ (qmyyx)V4 + (qmyxy)V4 + (qmxyy)V4 ] ∆y2∆x 288 − [ (qmxxy)V4 + (qmxyx)V4 + (qmyxχ)V4 ] ∆x2∆y 288 }. (11) The unit outward normal to the side surfaces V1V ′ 1V ′ 2V2 and V3V ′ 3V ′ 4V4 are (1, 0, 0) and (−1, 0, 0) respectively. Hence only space-time flux f passes through these planes. Each side surface is divided into two parts and the integration over each part is carried out using the third-order Taylor expansion of fluxes at the near solution point. Consider V1V ′ 1V4V ′ 4 as an example. It is the union of two polygons V1V ′ 1A1A ′ 1 andA1A ′ 1V4V ′ 4 . The integration over V1V ′ 1A1A ′ 1 uses the Taylor expansion of fm at V1 and the integration over A1A ′ 1V4V ′ 4 uses the Taylor expansion of fm at V4. The integration over the surface V1V ′ 1A1A ′ 1 implies following results:∫ j∆y (j− 1 2)∆y(k− 1 2)∆t { (fm)V1 + (fmy)V1 [ y − ( j − 1 2 ) ∆y ] + (fmt)V1 [ t− ( k − 1 2 ) ∆t ] A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1471 + 1 2 (fmyy)V1 [ y − ( j − 1 2 ) ∆y ]2 + 1 2 (fmtt)V1 [ t− ( k − 1 2 ) ∆t ]2 + (fmyt)V1 [ y − ( j − 1 2 ) ∆y ] [ t− ( k − 1 2 ) ∆t ] + 1 6 (fmyyy)V1 [ y − ( j − 1 2 ) ∆y ]3 + 1 6 (fmttt)V1 [ t− ( k − 1 2 ) ∆t ]3 + 1 2 (fmyyt)V1 [ y − ( j − 1 2 ) ∆y ]2 [ t− ( k − 1 2 ) ∆t ] + 1 2 (fmytt)V1 [ y − ( j − 1 2 ) ∆y ] [ t− ( k − 1 2 ) ∆t ]2 }dydt = ∆y∆t 4 { (fm)V1 + (fmy)V1 ∆y 4 + (fmt)V1 ∆t 4 + (fmyy)V1 ∆y2 24 + (fmtt)V1 ∆t2 24 + (fmyt)V1 ∆y∆t 16 + (fmyyy)V1 ∆y3 192 + (fmttt)V1 ∆t3 192 + (fmyt)V1 ∆y2∆t 96 + (fmytt)V1 ∆y∆t2 96 }, (12) and on the surface A1A ′ 1V4V ′ 4 , gives∫ (j+ 1 2)∆y j∆y ∫ k∆t (k− 1 2)∆t { (fm)V4 + (fmy)V4 [ y − ( j + 1 2 ) ∆y ] + (fmt)V4 [ t− ( k − 1 2 ) ∆t ] + 1 2 (fmyy)V4 [ y − ( j + 1 2 ) ∆y ]2 + 1 2 (fmtt)V4 [ t− ( k − 1 2 ) ∆t ]2 + (fmyt)V4 [ y − ( j + 1 2 ) ∆y ] [ t− ( k − 1 2 ) ∆t ] (13) + 1 6 (fmyyy)V4 [ y − ( j + 1 2 ) ∆y ]3 + 1 6 (fmttt)V4 [ t− ( k − 1 2 ) ∆t ]3 + 1 2 (fmyyt)V4 [ y − ( j + 1 2 ) ∆y ]2 [ t− ( k − 1 2 ) ∆t ] + 1 2 (fmytt)V4 [ y − ( j + 1 2 ) ∆y ] [ t− ( k − 1 2 ) ∆t ]2 }dydt = ∆y∆t 4 { (fm)V4 ∆y 4 + (fmt)V4 ∆t 4 + (fmyy)V4 ∆y2 24 + (fmtt)V4 ∆t2 24 − (fmyt)V4 ∆y∆t 16 − (fmyyy)V4 ∆y3 192 + (fmttt)V4 ∆t3 192 + (fmyyt)V4 ∆y2∆t 96 − (fmytt)V4 ∆y∆t2 96 }. (14) The integration on the surface V2V ′ 2A ′ 3A3 and A3A ′ 3V ′ 3V3 is similar to that on V1V ′ 1A ′ 1A1 and A1A ′ 1V ′ 1V1 with the only difference that the outward normal vector is on opposite A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1472 direction and fm is approximated at points V2 and V4 respectively. The outward unit normal to surface V1V ′ 1V ′ 2V2 and V3V ′ 3V ′ 4V4 are (0, 1, 0) and (0,−1, 0) respectively. Hence only space-time flux g passes through these surfaces. Therefore the integration over these surfaces is similar to the surfaces V2V ′ 2A ′ 3A3, A3A ′ 3V ′ 3V3, V1V ′ 1A ′ 1A1 and A1A ′ 1V ′ 1V1 with fm replaced by gm and y replaced by x. The outward unite vector to the top surface is (0, 0, 1) but the flux qm is approximated by Taylor expansion at the solution point C ′1. Hence the integration over the top surface results as follows:∫ (i+ 1 2)∆x(j+ 1 2)∆y (i− 1 2)∆x(j− 1 2)∆y { (qm)C′1 + ( qmx)C′1 [x− i∆x] + (qmy)C′1 [y − j∆y] + 1 2 (qmxχ)C′1 [x− i∆x]2 + 1 2 (qmyy)C′1 [y − j∆y]2 + 1 6 (qmxxx)C′1 [x− i∆x]3 + 1 6 (qmyyy)C′1 [y − j∆y]3 + 1 6 (qmxxx)C′1 [x− i∆x]3 + 1 6 (qmyyy)C′1 [y − j∆y]3 + 1 6 [ (qmyyx)C′1 + (qmyxy)C′1 + (qmxyy)C′1 ] [x− i∆x][y − j∆y]2 + 1 6 [ (qmyxχ)C′1 + (qmxyx)C′1 + (qmxxy)C′1 ] [x− i∆x]2[y − j∆y]}dxdy = ∆x∆y { (qm)C′1 + (qmxx)C′1 ∆x2 24 + (qmyy)C′1 ∆y2 24 } . (15) Combining all of above integration gives the following time marching scheme: (qm)ni,j = 1 ∆x∆y 4∑ l=1 〈 ∆x∆y 4 { (qm) n−1/2 l + (qmx) n−1/2 l ( kx ∆x 4 ) + (qmy) n−1/2 l ( ky ∆y 4 ) + 1 24 [ (qmxx) n−1/2 l (kx∆x)2 + (qmyy) n−1/2 l (ky∆y)2 ] + 1 32 [ (qmxy) n−1/2 l + (qmyx) n−1/2 l ] (kx∆x) (ky∆y) + 1 192 [ (qmxxx) n−1/2 l (kx∆x)3 + (qmyyy) n−1/2 l (ky∆y)3 ] + (kx∆x)2 (ky∆y) 288 [ (qmxxy) n−1/2 l + (qmxyx) n−1/2 l + (qmyxx) n−1/2 l ] + (ky∆y)2 (kx∆x) 288 [ (qmyyx) n−1/2 l + (qmyxy) n−1/2 l + (qmxyy) n−1/2 l ] } + kx ∆y∆t 4 { (fm) n−1/2 l + (fmy) n−1/2 l ( ky ∆y 4 ) + (fmt) n−1/2 l ( ∆t 4 ) + 1 24 [ (fmtt) n−1/2 l (∆t)2 + (fmyy) n−1/2 l (ky∆y)2 ] + 1 16 (fmyt) n−1/2 l (ky∆y) ∆t A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1473 + 1 192 [ (fmyyy) n−1/2 l (ky∆y)3 + (fmttt) n−1/2 l (∆t)3 ] + 1 96 [ (fmyyt) n−1/2 l (ky∆y)2 ∆t+ (fmytt) n−1/2 l (ky∆y) (∆t)2 ] } + ky ∆x∆t 4 { (gm) n−1/2 l + (gmx) n−1/2 l ( kx ∆x 4 ) + (gmt) n−1/2 l ( ∆t 4 ) + 1 24 [ (gmtt) n−1/2 l (∆t)2 + (gmxχ) n−1/2 l (kx∆x)2 ] + 1 16 (gmxt) n−1/2 l (kx∆x) ∆t + 1 192 [ (gmxxx) n−1/2 l (kx∆x)3 + (gmttt) n−1/2 l (∆t)3 ] + 1 96 [ (gmxxt) n−1/2 l (kx∆x)2 ∆t+ (gmxtt) n−1/2 l (kx∆x) (∆t)2 ] }) − 1 24 [ (qmxx)nC1 (∆x)2 + (qmyy) n C1 (∆y)2 ] . (16) where m = 1, 2, · · · , 5 is the number of equations, l = 1, 2, 3, 4 represents the four nodes V1, V2, V3, V4 and kx, ky are 1 or -1 depending upon the position of solution point C1 relatives four nodes V1, V2, V3, V4. 2.1. Higher order derivatives of flux function fm and gm The space and time derivatives (fmξ, fmξξ, fmξξξ, gmξ, gmξξ, gmξξξ, ξ = x, y, t; ) etc are obtained according to the chain rule. The detail of these derivatives is presented below only for the flux function fm. The higher order derivatives for gm are obtained in the same manner and hence the details are omitted. ∂fm ∂ξ = 5∑ r=1 ∂fm ∂qr ∂qr ∂ξ , (17) Differentiating again and applying the chain rule gives the following second order flux component’s derivatives: ∂2fm ∂ξ∂φ = 5∑ r=1 ∂fm ∂qr ∂2qr ∂ξ∂φ + 5∑ r=1 5∑ s=1 ∂2fm ∂qr∂qs ∂qr ∂ξ ∂qs ∂φ , ξ = x, y, t;φ = x, y, t, (18) and Similarly the third-order derivatives can be calculated using the chain rule. The third order derivatives are given below. fmxxx = 5∑ p=1 ∂fm ∂qp .qpξφχ + 5∑ r=1 5∑ s=1 ∂2fm ∂qr∂qs (qrξφ.qsχ + qrξχ.qsφ + qrφχqsξ) (19) + 5∑ r=1 5∑ s=1 5∑ t=1 ∂3fm ∂qr∂qs∂qt qrξ.qsφ.qtχ, ξ = x, y, t;φ = x, y, t;χ = x, y, t, (20) A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1474 The partial derivatives of flux components ∂fm ∂qr , ∂2fm ∂qr∂qs and ∂3fm ∂qr∂qs∂qt are the entries of the Jacobian matrices. For shallow water MHD equations, these derivatives are obtained by Matlab symbolic computations specially for the third order case. The computation of time derivatives of flux variables are obtained by using the conservation law. qmt = −fmx − gmy. (21) Taking derivative with respect to space variable x on both sides and treating the mixed space and time derivatives same allows to write down qmxt = −fmxx − gmyx. (22) Similarly, for time variable t qmtt = −fmxt − gmyt. (23) Using the values of partial derivatives of fm and gm provides the values of qmtt which gives the values for fmtt according to above equations. Again taking the time and space derivatives on both sides of equation 23 follows qmxtt = −fmxxt − gmyxt, (24) and qmttt = −fmxtt − gmytt. (25) Using the values of partial derivatives of fm and gm provides the values of qmtt which gives the values for fmttt according to above equations. Similarly we can find the higher order time derivatives for gm only be changing the differentiation with respect to y in place of x. 2.2. Spatial Derivatives of Conserved variables At this point, all the time derivatives of conservative variables are available through the derivatives of flux fm and the flux derivatives. Now at the end it is only left to get the space derivative variables of qm such as qmξ, qmξφ, and qmξφχ ξ = x, y; φ = x, y; χ = x, y. The first derivative qmξ is updated at the ith node using the updated values of qm obtained from Equation 21. In general each higher order space derivatives of conserved variables is updated using the updated values of next higher order derivative. Now, for the fourth- order scheme, the values of two more derivatives qmξφ and qmξφχ are required to update. Proceeding in following steps allows to get the updated values of all the derivatives: 1: updating qmξφχ at Ci; 2: updating qmξφ at Ci; 3: updating qm at Ci; 4: updating qmξ at Ci. . The details are omitted and the reader is referred to the work of Yang Yun et al [23] for details. A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1475 2.3. Divergence Constraint The SWMHD system is augmented by the divergence constraint ∇ · (hB). The nu- merical solution obtained from the CESE method must satisfy this equation. Hence a divergence cleaning procedure is required to embed in the CESE method. The projection method has been adopted for this purpose as seen in the work [2],[16] and [20]. Con- sider the updated values (hn+1, (hu1)n+1, (hu2)n+1, (hB1)n+1, (hB2)n+1). It readily gives B∗1 = (hB1)n+1 h(n+1) and B∗2 = (hB2)n+1 h(n+1) . Define B∗ = (B∗1 , B ∗ 2) as the vector of magnetic field components. This is decomposed into curl of magnetic potential A and divergence of potential function φ as B∗ = 5×A+5φ. (26) Taking the divergence on both sides of above equation gives 52φ = 5 ·B∗. (27) Solving equation 27 for φ gives the divergence free magnetic field at time level t = tn+1 by Bn+1 = B∗ −5φ. 3. Numerical Test caes This section presents one and two dimensional test problems taken from literature. Both 2nd-order and 4th-order CESE schemes have been applied to these problem. 3.1. One Dimensional Riemann Problem A numerical test is performed to demonstrate the robustness and the accuracy of 4th- order CESE scheme over 2nd-order CESE schem. A strong Riemann problem has been solved here. This problem is taken from [22]. The initial conditions are given as: h v1 v2 B1 B2  = { [1, 0, 0, 1, 0]T , if x ≤ 0 [2, 0, 0, 0.5, 1]T , if x > 0 (28) The problem is solved with periodic boundary conditions and domain Ω = [−1, 1] at t = 0.4. The results are shown in Figures 3- 8. The L2 errors and experimental order of convergence for 4th-order CESE scheme have been listed in table have been listed in table 1. As expected, it is clear seen that the 4th-order CESE scheme with 500 grid points shows remarkable details in the solution compared to the 2nd-order CESE scheme at grid of 1000 nodes. The zoomed plot reveals better performance of 4th order CESE scheme near discontinuities as compared to the second order CESE scheme. The execution time for both schemes is observed on intel core i7 8th Generation machine. It comes out that A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1476 the second order CESE scheme takes 8 seconds and the 4th order CESE scheme takes 12 seconds to produce these results. Hence the 4th order CESE scheme is computationally more efficient with higher order accuracy. Table 1: Error Analysis of 4th-order CESE Scheme Number of Grid points L2-Error Experimental Order of Convergence 50 0.00435 - 100 0.000220 4.3 200 0.000011932 4.23 400 0.000001015 4.19 Hence it concludes that the 4th-order CESE scheme advantages over the 2nd-order CESE scheme by showing much better results. It follows that the 4th-order scheme is computationally more efficient and accurate. 4. Two Dimensional Riemann Problem Simulation The numerical test has been performed to demonstrate the robustness and the accuracy of 4th-order CESE scheme over 2nd-order CESE schem. Problem 1: Multidimensional Steady-State Shock This problem is taken from [22] in which authors have applied the 2nd-order CESE method to this problem. The computational domain (x, y) ∈ [−1, 1] × [−1, 1] is divided into 200 × 200 grid with g = 1 and the initial data is given as: h u1 u2 B1 B2  = { [1, 4.5, 0, 2, 0]T , if y ≤ 0 [2, 5.5, 0, 0.5, 0]T , if y > 0 (29) The problem is solved with outflow boundary conditions at the top, bottom and right boundaries whereas inflow boundary condition at the left boundary specified with above mentioned initial data. The results at time t = 0.25 are shown in Figures 9 - 11. As expected, it is clear seen that the fourth-order scheme shows remarkable details in the solution compared to the second -order scheme at same grid. The results from higher order scheme are much more dense as compared to the 2nd-order scheme. 5. Conclusions The fourth-order CESE method is presented in detail for SWMHD equations and tested for a problem taken from literature. The benchmark test revealed excellent characteristics A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1477 [H] time 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 H e ig h t 0.95 1 1.05 1.1 4th order CESE method 2nd order CESE method Figure 4: Zoomed Plot for Height at t = 0.4 of the method. The scheme is based on unified treatment of space and time variables that allows a systematic and easy extension to higher dimensions as well as to arbitrary higher order scheme. The scheme is computationally efficient and accurate and outperforms the second-order CESE method even with less grid. Acknowledgements This work has been carried out as a part of project under the grant by Higher Education Commission of Pakistan vide office letter No. 6255/Sindh/NRPU/R&D/HEC/2016. I also acknowledge Sukkur IBA University for providing state of the art research environment. I would like to acknowledge my MS student Irfan Hyder for his help in editing some part of this manuscript. x-axis 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 H e ig h t 0.8 1 1.2 1.4 1.6 1.8 2 4th order CESE 2nd order CESE scheme Figure 3: Height at t = 0.4. A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1478 [H] x-axis 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 v 1 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 4th order CESE scheme 2nd order CESE scheme Figure 5: v1 at time t = 0.4. [H] x-axis 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 v 2 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 4th-order CESE scheme 2nd-order CESE scheme Figure 6: v2 at time t = 0.4. x-axis -1 1 y- ax is -1 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 (a) Results of 4th-order CESE method x-axis -1 1 y- ax is -1 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 (b) Results of 2nd-order CESE method Figure 9: Height h computed at t = 0.25 A. Sidrah, Z. Saqib / Eur. J. Pure Appl. Math, 12 (4) (2019), 1464-1482 1479 [H] x-axis 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 B 1 0.5 0.6 0.7 0.8 0.9 1 1.1 4th-order CESE scheme 2nd-order CESE scheme Figure 7: B1 at time t = 0.4. [H] x-axis 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 B 2 0 0.2 0.4 0.6 0.8 1 1.2 4th-order CESE scheme 2nd-order CESE scheme Figure 8: B2 at time t = 0.4. x-axis -1 1 y- ax is -1 1 5 6 7 8 9 10 (a) Results of 4th-order CESE method x-axis -1 1 y- ax is -1 1 5 6 7 8 9 10 (b) Results of 2nd-order CESE method Figure 10: ux computed at t = 0.25 REFERENCES 1480 x-axis -1 1 y- ax is -1 1 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 (a) Results of 4th-order CESE method x-axis -1 1 y- ax is -1 1 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 (b) Results of 2nd-order CESE method Figure 11: bx computed at t = 0.25 References [1] Balsara, Dinshaw S., Michael Dumbser, and Remi Abgrall. ”Multidimensional HLLC Riemann solver for unstructured mesheswith application to Euler and MHD flows.” Journal of Computational Physics 261 (2014): 172-208. [2] Brackbill, Jeremiah U., and Daniel C. Barnes. ”The effect of nonzero 5·B on the nu- merical solution of the magnetohydrodynamic equations.” Journal of Computational Physics 35, no. 3 (1980): 426-430. [3] Chang, Sin-Chung. ”The method of space-time conservation element and solution elementA new approach for solving the Navier-Stokes and Euler equations.” Journal of computational Physics 119, no. 2 (1995): 295-324. [4] Chang, S.C.; Yu, S.T.; Himansu, A.; Wang, X.Y.; Chow, C.Y.; Loh, C.Y. The method of space-time conservation element and solution elementA new paradigm for numerical solution of conservation laws. Comput. Fluid Dyn. Rev. 1998, 2, 206–240. [5] De Sterck, H. Hyperbolic theory of the shallow water magnetohydrodynamics equa- tions. Phys. Plasmas 2001, 8, 3293–3304. [6] Dikpati, M.; Gilman, P.A. Flux-transport dynamos with α-effect from global insta- bility of tachocline differential rotation: A solution for magnetic parity selection in the Sun. Astrophys. J. 2001, 559, 428. [7] Dikpati, M.; Gilman, P.A. Prolateness of the solar tachocline inferred from latitudinal force balance in a magnetohydrodynamic shallow-water model. Astrophys. J. 2001, 552, 348. [8] Fuchs, Franz Georg, Andrew D. McMurry, Siddhartha Mishra, Nils Henrik Rise- bro, and Knut Waagan. ”Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations.” Communications in Computational Physics 9, no. 2 (2011): 324-362. [9] Gilman, P.A. Magnetohydrodynamic shallow water equations for the solar tachocline. Astrophys. J. Lett. 2000, 544, L79. REFERENCES 1481 [10] Sidrah, A. ”Numerical Approximation of Classical and Relativistic Magnetohydro- dynamics.” PhD diss., COMSATS Institute of Information Technology, Islamabad- Pakistan, 2011. [11] Jiang, C.; Feng, X.; Zhang, J.; Zhong, D. AMR simulations of magnetohydrodynamic problems by the CESE method in curvilinear coordinates. Sol. Phys. 2010, 267, 463– 491. [12] Krger, T.; Lukov-Medvidov, M. An evolution Galerkin scheme for the shallow water magnetohydrodynamic equations in two space dimensions. J. Comput. Phys. 2005, 206, 122–149. [13] Massimi, H.S.; Shen, H.; Wen, C.Y. Study of Hypersonic Dissociating Flows over Spheres Using the Space-Time CE/SE Method. In Proceedings of the 30th Inter- national Symposium on Shock Waves, Tel Aviv, Israel, 19–24 July 2015; Springer: Cham, Switzerland, 2017; pp. 145–148. [14] Qamar, S.; Mudasser, S. On the application of a variant CE/SE method for solving two-dimensional ideal MHD equations. Appl. Numer. Math. 2010, 60, 587–606. [15] Qamar, S.; Warnecke, G. Application of spacetime CE/SE method to shallow water magnetohydrodynamic equations. J. Comput. Appl. Math. 2006, 196, 132–149. [16] Rossmanith, J.A. A wave propagation method with constrained transport for ideal and shallow water magnetohydrodynamics. Ph.D. Thesis, University of Washington, Washington, DC, USA, 2002. [17] Rossmanith, J.A. A constrained transport method for the shallow water MHD equations. In Hyperbolic Problems: Theory, Numerics, Applications; Springer: Berlin/Heidelberg, Germany, 2003; pp. 851–860. [18] Shen, H., Wen, C.Y. and Zhang, D.L. A characteristic spacetime conservation element and solution element method for conservation laws. Comput. Phys. J. 2015, 288, 168– 182. [19] Shen, H. and Wen, C.Y. A characteristic spacetime conservation element and solution element method for conservation laws II. Multidimensional extension. Comput. Phys. J. 2016, 305, 775–792. [20] Tth, Gbor. ”The 5·B constraint in shock-capturing magnetohydrodynamics codes.” Journal of Computational Physics 161, no. 2 (2000): 605-652. [21] Wang, Gang, Deliang Zhang, Kaixin Liu, and Jingtao Wang. ”An improved CE/SE scheme for numerical simulation of gaseous and two-phase detonations.” Computers and Fluids 39, no. 1 (2010): 168-177. [22] Winters, A.R.; Gassner, G.J. An entropy stable finite volume scheme for the equations of shallow water magnetohydrodynamics. J. Sci. Comput. 2016, 67, 514–539. REFERENCES 1482 [23] Yang, Y.; Feng, X.S.; Jiang, C.W. A high-order CESE scheme with a new divergence- free method for MHD numerical simulation. J. Comput. Phys. 2017, 349, 561–581. [24] Yang, Y., Feng, X.S. and Jiang, C.W. An upwind CESE scheme for 2D and 3D MHD numerical simulation in general curvilinear coordinates. Comput. Phys. J. 2018, 371, 850–869. [25] Yu, Sheng-Tao John, Lixiang Yang, Robert L. Lowe, and Stephen E. Bechtel. ”Nu- merical simulation of linear and nonlinear waves in hypoelastic solids by the CESE method.” Wave Motion 47, no. 3 (2010): 168-182. [26] Zhang, M., Yu, S.T.J., Lin, S.C.H., Chang, S.C. and Blankson, I. Solving the MHD equations by the spacetime conservation element and solution element method. Com- put. Phys. J. 2006, 214, 599–617.