# Parallel Pseudo Arc-Length Moving Mesh Schemes for Multidimensional Detonation

Authors: Jianguo Ning; Xinpeng Yuan; Tianbao Ma; Jian Li
Journal: Scientific Programming
 (2017)

---

## Abstract

We have discussed the multidimensional parallel computation for pseudo arc-length moving mesh schemes, and the schemes can be used to capture the strong discontinuity for multidimensional detonations. Different from the traditional Euler numerical schemes, the problems of parallel schemes for pseudo arc-length moving mesh schemes include diagonal processor communications and mesh point communications, which are illustrated by the schematic diagram and key pseudocodes. Finally, the numerical examples are given to show that the pseudo arc-length moving mesh schemes are second-order convergent and can successfully capture the strong numerical strong discontinuity of the detonation wave. In addition, our parallel methods are proved effectively and the computational time is obviously decreased.


---

## Body

## 1. Introduction

Detonation is a type of combustion involving a supersonic exothermic front accelerating through a medium that eventually drives a shock front propagating directly in front of it. Detonations occur in both conventional solid and liquid explosives, [1] as well as in reactive gases [2]. Due to the strong discontinuity, one of the most important challenges in detonation simulations is to capture the precursor shock wave. To capture the strong discontinuity in detonation waves, mesh adaptation is an indispensable tool for use in the efficient numerical solution of this type of problem. The moving mesh method is one of the mesh adaptation methods, which relocate mesh point positions while maintaining the total number of mesh points and the mesh connectivity [3–5]. Particularly, Tang et al. proposed a moving mesh method that contained two parts: physical PDE time evolution and mesh redistribution [6–8]. The physical PDE time evolution and mesh redistribution were alternating, and conservative interpolation was used to transfer solutions from old mesh to new mesh. The method is shown to work well generally for hyperbolic conservation. After that, Ning et al. have improved this method and proposed the pseudo arc-length moving mesh schemes, which can deal with the multidimensional chemical reaction detonation problem [9].In this paper, we will discuss the parallel computation of pseudo arc-length moving mesh schemes for multidimensional detonation. There are some works about the parallel schemes for Euler numerical scheme. Knepley and Karpeev developed a new programming framework, called Sieve, to support parallel numerical partial differential equation (PDE) algorithms operating over distributed meshes [10]. Morris et al. demonstrated how to build a parallel application that encapsulates the Message Passing Interface (MPI) without requiring the user to make direct calls to MPI except for startup and shutdown [11]. Mubarak et al. present us with a parallel algorithm of creating and deleting data copies, referred to as ghost copies, which localize neighborhood data for computation purposes while minimizing interprocess communication [12]. Wang et al. use the parallel scheme to reduce the cost for adaptive mesh refinement WENO scheme of multidimensional detonation [13]. However, there are no discussions about the parallel computation for moving mesh schemes, which will be of concern in this paper. Different from the traditional Euler numerical scheme, the data communications of moving mesh schemes between processors are more complex, which include physical values and mesh points. Besides, the processor communications for pseudo arc-length moving mesh schemes include adjacent processor and diagonal processor. Here, we adopt the software architecture of MPI. The most important advantages of this model are twofold: achievable performance and portability. Performance is a direct result of available optimized MPI libraries and full user control in the program development cycle. Portability arises from the standard API and the existence of MPI libraries on a wide range of machines. In general, an MPI program runs on distributed memory machines. The processor communications for the parallel computation of pseudo arc-length moving mesh schemes are more complex than the traditional Euler scheme, and it includes adjacent processor and diagonal processor. Here, we will consider three kinds of processor partitions to show our parallel schemes. This article is organized as follows. Section 2 introduces the chemical and physical model. Section 3 presents the numerical scheme. Section 4 is devoted to the parallel computation. Section 5 conducts several numerical experiments to demonstrate our schemes. The paper ends with a conclusion and discussion in Section 6.


## 2. Governing Equations

Instead of using many real elementary reactions, a two-step model was utilized as the testing model. Two-step reaction model considers a complicated chemical reaction to be an induction and an exothermic reaction. For both induction reaction and exothermic reaction, the progress parametersα and β are unity at first, then α decreases to zero, and β decreases until an equilibrium state is reached. The rates ωα and ωβ are given as follows [14].(1)ωα=dαdt=-kαρexp⁡-EαRT,ωβ=dβdt=0α>0,-kβp2β2exp⁡-EβRT-1-β2exp⁡-Eβ+QRTα≤0,where ρ is the mass density, p the pressure, T the temperature, R the gas constant, Q the heat release parameter, kα and kβ the constants of reaction rates, and Eα and Eβ the activation energies.In deriving fundamental equations, the gas is assumed to be perfect, nonviscous, and non-heat-conducting. In Cartesian coordinates, governing equations for gaseous detonation problem, including the above chemical reaction, are(2)∂∂tw+∇xFw=sw,t≥0,where x is the multidimensional vector, Fw is the multidimensional matrix function, and sw is the chemical reaction source term. For the one-dimensional space,(3)w=ρ,ρu,E,ρα,ρβT,F=ρu,ρu2+p,E+pu,ρuα,ρuβT,s=0,0,0,ω1,ω2T.For the two-dimensional space,(4)w=ρρuρvEραρβ,F=ρuρvρu2+pρuvρuvρv2+pE+puE+pvρuαρvαρuβρvβ,s=0000ω1ω2.In the case of three-dimensional space,(5)w=ρρuρvρwEραρβ,F=ρuρvρwρu2+pρuvρuwρvuρv2+pρvwρwuρwvρw2+pE+puE+pvE+pwρuαρvαρwαρuβρvβρwβ,s=00000ω1ω2,where u, v, and w are the Cartesian components of the fluid velocity in the x, y, and z directions, respectively. Total energy density E is defined as(6)E=pγ-1+12ρu2+v2+w2+ρQβ.Here γ is the specific heat ratio.


## 3. Numerical Method

Firstly, we present the framework of pseudo arc-length moving mesh schemes for the gaseous detonation problem (2). Our adaptive scheme is formed by two independents parts: the evolution of the governing equation and the iterative mesh redistribution.


### 3.1. Time Evolution of Governing Equations

The physical domain for computation isΩp. Given a partition of the mesh domain Kj¯∣Kj¯∈Ωp,Kj¯=xj,xj+1,⋃j¯Kj¯=Ωp and a partition of the time interval 0,T. The average on the cell Kj¯ is(7)w¯j¯=1Kj¯∫Kj¯wj¯dσ.Here, the sign Kj¯ denotes the length for the region Kj¯Kj¯ in one dimension, the area for the region Kj¯Ki¯,j¯ in two dimensions, and the volumes for the region Kj¯Ki¯,j¯,k¯ in three dimensions. Integrating (2) over Kj¯, we have(8)∫Kj¯∂w∂tdσ+∫Kj¯∇xFwdσ=∫Kj¯swdσ.Applying the Green-Gauss’s theorem and rewriting, we have(9)Kj¯dwj¯dt+∮∂Kj¯Fw·nj¯dσ=Kj¯swj¯,where nj¯ is the outward unit normal vector of boundary external surface ∂Kj¯. The Lax-Friedrichs flux is defined by(10)Hu,v,n=12Fu·n+Fv·n-av-u,where a=maxu,vF′u·n. It satisfies the conservation and consistency(11)Hu,v,n=-Hu,v,-n,Hu,u,n=Fu·n.Let ej¯l∣ej¯l∈∂Kj¯,⋃l=1Nlej¯l=∂Kj¯ be a partition of boundary external surface ∂Kj¯ and nj¯ be the outward unit normal vector of surface ej¯l. Then we have a semidiscrete scheme of (2)(12)Kj¯dwj¯dt=-∑l=1Nl∫ej¯lHwj¯intl,wj¯extl,nj¯ldσ+Kj¯swj¯,where wj¯intl or wj¯extl, l=1,2,…,Nl, is the internal or external approximate value at ej¯l.For the one-dimensional space,Nl=2, and ej¯1ej¯2 is the left (right) point of line cell Kj¯. Thus, (12) becomes(13)dwj¯dt=1Kj¯-Hwj¯int1,wj¯ext1+Hwj¯int2,wj¯ext2+swj¯.By the initial reconstruction technique [15] to reset wj¯intl, wj¯extl, l=1,2, on the edge ej¯l, we can obtain the second-order accurate spatial discretization:(14)wj¯int1=wj¯+12ψwj¯-wj¯+1,wj¯-1-wj¯,wj¯ext1=wj¯-1+12ψwj¯-wj¯-1,wj¯-1-wj¯-2,wj¯int2=wj¯+1ext1,wj¯ext2=wj¯+1int1,where ψ is a nonlinear limiter function which is used to suppress the possible pseudo oscillation. In our computations, we use van Leer’s limiter [15](15)ψa,b=signa+signbaba+b+ε.For the two-dimensional space, Nl=4, and the diagram for edges ei¯,j¯l, l=1,2,3,4, of the quadrilateral element Ki¯,j¯ is shown in Figure 1.Figure 1
Schematic diagram of the quadrilateral elementKi¯,j¯.The valueswi¯,j¯intl, wi¯,j¯extl at edge ei¯,j¯l, l=1,2,3,4, can be defined by(16)wi¯,j¯int1=wi¯,j¯+12ψwi¯,j¯-wi¯,j¯+1,wi¯,j¯-1-wi¯,j¯,wi¯,j¯ext1=wi¯,j¯-1+12ψwi¯,j¯-wi¯,j¯-1,wi¯,j¯-1-wi¯,j¯-2,wi¯,j¯ext3=wi¯,j¯+1int1,wi¯,j¯int3=wi¯,j¯+1ext1wi¯,j¯int2=wi¯,j¯+12ψwi¯,j¯-wi¯-1,j¯,wi¯+1,j¯-wi¯,j¯,wi¯,j¯ext2=wi¯+1,j¯+12ψwi¯+1,j¯-wi¯+2,j¯,wi¯,j¯-wi¯+1,j¯,wi¯,j¯ext4=wi¯-1,j¯int2,wi¯,j¯int4=wi¯-1,j¯ext2.For the three-dimensional hexahedron element, the situation is complex. Due to the movement of the mesh points, the hexahedron element will change to polyhedron element. For calculating the element volume and the boundary external surface area, each surface of hexahedron element is divided into two surfaces and calculated. The control element Ki¯,j¯,k¯ is illustrated in Figure 2. The boundary surfaces ei¯,j¯,k¯l, l=1,2,…,12, are defined in Table 1.Table 1
Definition of boundary surfacesei¯,j¯,k¯l, l=1,2,…,12.
Surfaces
Vertexes
ei-,j-,k-1
x
i
+
1
,
j
,
k
+
1
,
x
i
+
1
,
j
,
k
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-2
x
i
+
1
,
j
,
k
,
x
i
+
1
,
j
+
1
,
k
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-3
x
i
,
j
,
k
,
x
i
,
j
,
k
+
1
,
x
i
,
j
+
1
,
k
+
1
ei-,j-,k-4
x
i
,
j
,
k
,
x
i
,
j
+
1
,
k
+
1
,
x
i
,
j
+
1
,
k
ei-,j-,k-5
x
i
,
j
+
1
,
k
,
x
i
,
j
+
1
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-6
x
i
,
j
+
1
,
k
,
x
i
+
1
,
j
+
1
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
ei-,j-,k-7
x
i
,
j
,
k
+
1
,
x
i
,
j
,
k
,
x
i
+
1
,
j
,
k
+
1
ei-,j-,k-8
x
i
,
j
,
k
,
x
i
+
1
,
j
,
k
,
x
i
+
1
,
j
,
k
+
1
ei-,j-,k-9
x
i
,
j
,
k
+
1
,
x
i
+
1
,
j
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-10
x
i
,
j
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
+
1
,
x
i
,
j
+
1
,
k
+
1
ei-,j-,k-11
x
i
+
1
,
j
,
k
,
x
i
,
j
,
k
,
x
i
+
1
,
j
+
1
,
k
ei-,j-,k-12
x
i
,
j
,
k
,
x
i
,
j
+
1
,
k
,
x
i
+
1
,
j
+
1
,
kFigure 2
Schematic diagram of the quadrilateral elementKi¯,j¯,k¯.The second-order valueswi¯,j¯,k¯intl, wi¯,j¯,k¯extl at edge ei¯,j¯,k¯l, l=1,2,…,12, can be obtained by(17)wi¯,j¯,k¯int1=wi¯,j¯,k¯int2=w¯i¯,j¯,k¯+12ψw¯i¯,j¯,k¯-w¯i¯-1,j¯,k¯,w¯i¯+1,j¯,k¯-w¯i¯,j¯,k¯,wi¯,j¯,k¯ext1=wi¯,j¯,k¯ext2=w¯i¯+1,j¯,k¯+12ψw¯i¯+1,j¯,k¯-w¯i¯+2,j¯,k¯,w¯i¯,j¯,k¯-w¯i¯+1,j¯,k¯,wi¯,j¯,k¯int5=wi¯,j¯,k¯int6=w¯i¯,j¯,k¯+12ψw¯i¯,j¯,k¯-w¯i¯,j¯-1,k¯,w¯i¯,j¯+1,k¯-w¯i¯,j¯,k¯,wi¯,j¯,k¯ext5=wi¯,j¯,k¯ext6=w¯i¯,j¯+1,k¯+12ψw¯i¯,j¯+1,k¯-w¯i¯,j¯+2,k¯,w¯i¯,j¯,k¯-w¯i¯,j¯+1,k¯,wi¯,j¯,k¯int9=wi¯,j¯,k¯int10=w¯i¯,j¯,k¯+12ψw¯i¯,j¯,k¯-w¯i¯,j¯,k¯-1,w¯i¯,j¯,k¯+1-w¯i¯,j¯,k¯,wi¯,j¯,k¯ext9=wi¯,j¯,k¯ext10=w¯i¯,j¯,k¯+1+12ψw¯i¯,j¯,k¯+1-w¯i¯,j¯,k¯+2,w¯i¯,j¯,k¯-w¯i¯,j¯,k¯+1,wi¯,j¯,k¯int3=wi¯-1,j¯,k¯ext1,wi¯,j¯,k¯int4=wi¯-1,j¯,k¯ext2,wi¯,j¯,k¯ext3=wi¯-1,j¯,k¯int1,wi¯,j¯,k¯ext4=wi¯-1,j¯,k¯int2,wi¯,j¯,k¯int7=wi¯,j¯-1,k¯ext5,wi¯,j¯,k¯int8=wi¯,j¯-1,k¯ext6,wi¯,j¯,k¯ext7=wi¯,j¯-1,k¯int5,wi¯,j¯,k¯ext8=wi¯,j¯-1,k¯int6,wi¯,j¯,k¯int9=wi¯,j¯,k¯-1ext11,wi¯,j¯,k¯int9=wi¯,j¯,k¯-1ext11,wi¯,j¯,k¯ext10=wi¯,j¯,k¯-1int12,wi¯,j¯,k¯ext10=wi¯,j¯,k¯-1int12.Time discretization can be achieved by the strong stability preserving high order Runge-Kutta time discretization [16, 17]. The semidiscrete scheme (12) can be written as(18)w~t=Lw~.Here, the second-order TVD Runge-Kutta method for (18) in the time discretization is(19)w~1=w~n+ΔtLw~n,w~n+1=12w~n+12w~1+12ΔtLw~1.


### 3.2. Pseudo Arc-Length Moving Mesh Scheme

Letx=x1,x2,…,xNd and ξ=ξ1,ξ2,…,ξNd denote the physical and computational coordinates, respectively. Here, Nd denotes the number of spatial dimensions. A one-to-one coordinate transformation from the computational domain Ωc to the physical domain Ωp is denoted by(20)x:ξ⟶x,Ωc⟶Ωp.In the variational approach, the mesh map between the domains Ωc and Ωp is usually provided by(21)Ψξ=12∑r=1Nd∫Ωp∇xξrTGr-1∇xξrdx,where ∇x≔∂x1,∂x2,…,∂xNdT and Gr, r=1,2,…,Nd, are given symmetric positive definite matrices called monitor functions, depending on the underlying solution to be adaptive. The Euler-Lagrange equations of the functional (21) have the form(22)∇x·Gr-1∇xξr=0,1≤r≤Nd.In practice, Ωp may have a very complex geometry, and as a result it is not realistic to solve (22) directly. An alternative is to consider a functional(23)Ψ~x=12∑r=1Nd∫Ωc∇ξxrTGr∇ξxrdξ,where ∇ξ≔∂ξ1,∂ξ2,…,∂ξNdT. In addition, one of the choices of monitor functions is G=ωI, where I is the identity matrix and ω is a positive weight function. For the purpose to have more accuracy near the nonsmooth part of solutions, we introduce the monitor function of pseudo arc-length norm [18, 19](24)ω=1+α1W+α2∇W2,where α1 and α2 are some nonnegative constants. Thus, the corresponding Euler-Lagrange equations about (23) are(25)∇ξ·ω∇ξxr=0,1≤r≤Nd.In our computations, we will use the Gauss-Seidel type iteration method to solve the mesh equation (25)(26)∑r=1Ndϑj1,j2,…,jr+1/2,…,jNdrxj1,j2,…,jr+1,…,jNdκ-xj1,j2,…,jr,…,jNdκ+1-ϑj1,j2,…,jr-1/2,…,jNdrxj1,j2,…,jr,…,jNdκ+1-xj1,j2,…,jr-1,…,jNdκ=0,for 1≤jr≤Nξr, 1≤r≤Nd, and κ=0,1,…, where(27)ϑj1,j2,…,jr-1/2,…,jNdr=ωj1,j2,…,jr-1/2,…,jNd=12Nd-1∑sr=jr-1/2,sk=jk-1/2,jk+1/2ωs1,s2,…,sr,…,sNd,ϑj1,j2,…,jr+1/2,…,jNdr=ωj1,j2,…,jr+1/2,…,jNd=12Nd-1∑sr=jr+1/2,sk=jk-1/2,jk+1/2ωs1,s2,…,sr,…,sNd,where ωj¯1,j¯2,…,j¯Nd≔ωwj¯1,j¯2,…,j¯Ndκ. In addition, let us consider the solutions updating on new grids xκ+1 from the old grids xκ. Here, the solution updating should preserve conservative property for wj¯1,j¯2,…,j¯Nd in the sense of that(28)∑j¯1,j¯2,…,j¯NdKj¯1,j¯2,…,j¯Ndκwj¯1,j¯2,…,j¯Ndκ=∑j¯1,j¯2,…,j¯NdKj¯1,j¯2,…,j¯Ndκ+1wj¯1,j¯2,…,j¯Ndκ+1.Let Dl denote the region scanned by the boundary ej¯1,j¯2,…,j¯Ndl after one iterative step of (28). So we can have the following conservative interpolation scheme:(29)Kj¯1,j¯2,…,j¯Ndκ+1wj¯1,j¯2,…,j¯Ndκ+1=Kj¯1,j¯2,…,j¯Ndκwj¯1,j¯2,…,j¯Ndκ+∑l=1NlF⌢lwj¯1,j¯2,…,j¯Ndintl,wj¯1,j¯2,…,j¯Ndextl,where wj¯1,j¯2,…,j¯Ndintl and wj¯1,j¯2,…,j¯Ndextl are the reconstructed left and right states on the boundary ej¯1,j¯2,…,j¯Ndl by wj¯1,j¯2,…,j¯Ndκ (for more details see Section 3.2). F⌢l·,·, l=1,2,…,Nl, denotes the integration of w over the domain Dl. Here we have omitted the subscripts j¯r, r=1,2,…,Nd, of F⌢l·,· and Dl for simplicity. Since only cell averages of wj¯1,j¯2,…,j¯Nd over Kj¯1,j¯2,…,j¯Nd are known, we should give an approximation of F⌢l·,· instead of its exact value. Then F⌢l·,· can be approximately calculated as(30)F⌢lu,v=max⁡Dl,0·v+min⁡Dl,0·u.Next, we will show the examples to calculate Dl, l=1,2,…,Nd, according to different dimensions.In one-dimensional space,D1, D2 can be obtained from the following:(31)D1=xjκ-xjκ+1,D2=xj+1κ+1-xj+1κ.In two-dimensional space, the changed area Dl, l=1,2,…,4, are illustrated in Figure 3. As an example, we compute D1. Noticing xi,j=xi,j,yi,j, we have(32)D1≔12xi+1,jκ+1-xi,jκyi+1,jκ-yi,jκ+1-yi+1,jκ+1-yi,jκxi+1,jκ-xi,jκ+1.Obviously, D1 is the area, which is positive if xi-1,j-1z+1, xi,j-1z+1, xi,j-1z, and xi-1,j-1z are located in a counter-clockwise order; otherwise it is negative if those four points are located in a clockwise order.Figure 3
Schematic diagram of the two-dimensional changed areasDl, l=1,2,3,4.In the three-dimensional space, the control elementKi+1/2,j+1/2,k+1/2 is illustrated in Figure 4. Similarly, we compute D1 as an example. The changed area D1 is illustrated in Figure 5. D1 is composed by three tetrahedrons, so we have(33)D1=Vxi+1,j,k+1κxi+1,j,k+1κ+1xi+1,j+1,k+1κ+1xi,j,k+1κ+1+Vxi+1,j,k+1κxi+1,j+1,k+1κ+1xi,j,k+1κxi,j,k+1κ+1+Vxi+1,j,k+1κxi+1,j+1,k+1κ+1xi+1,j+1,k+1κxi,j,k+1κ.Each tetrahedron volume can be calculated by the volume formula(34)VOABC=16xOAyOAzOAxOByOBzOBxOCyOCzOC.Obviously, VOABC is the volume, which is positive if A, B, and C are located in a counter-clockwise order; otherwise it is negative if those three points are located in a clockwise order.Figure 4
Schematic diagram of the control elementKi+1/2,j+1/2,k+1/2.Figure 5
Three-dimensional changed areaD1.The solution procedure can be illustrated by the following flowchart.Algorithm 1.
Step  1. Give a uniform mesh ξj1,j2,…,jNd0 on the computational domain Ωc. Initialize the physical mesh xj1,j2,…,jNd0 on the physical domain Ωp, and compute the grid initial values wj¯1,j¯2,…,j¯Nd0. Then, let xj1,j2,…,jNd0≔xj1,j2,…,jNd0 and wj¯1,j¯2,…,j¯Nd0≔wj¯1,j¯2,…,j¯Nd0.
Step  2. With scheme (26), move the grid from xj1,j2,…,jNdκ to xj1,j2,…,jNdκ+1 and the corresponding grid values from wj¯1,j¯2,…,j¯Ndκ to wj¯1,j¯2,…,j¯Ndκ+1 based on scheme (29) for κ=0,1,2,…. Repeat the updating procedure for a fixed number of iterations or until xj1,j2,…,jNdκ+1-xj1,j2,…,jNdκ≤ε. Then, let x~j1,j2,…,jNd=xj1,j2,…,jNdκ+1 and w~j¯1,j¯2,…,j¯Nd=wj¯1,j¯2,…,j¯Ndκ+1 at time level tn.
Step  3. Evolve the underlying PDEs using finite volume scheme (12) and time discretization scheme (19) on the new mesh x~j1,j2,…,jNd to obtain the numerical approximations wj¯1,j¯2,…,j¯Ndn+1 at the time level tn+1.
Step  4. If tn+1⩽T, then let xj1,j2,…,jNd0=x~j1,j2,…,jNd and wj¯1,j¯2,…,j¯Nd0=wj¯1,j¯2,…,j¯Ndn+1, and go to Step  2.Remark 2.
In practice, it is common to use some temporal or spatial smoothing on the weight functionω to obtain smoother meshes in order to avoid very singular mesh and/or large approximation errors near those regions where the solution has a large gradient. A low pass filter should be used to smooth the weight function 2~3 times for each κ: in one-dimensional space:(35)ωj=14ωj-1+2ωj+ωj+1;in two-dimensional space:(36)ωi,j⟵14ωi,j+18ωi+1,j+ωi,j+1+ωi-1,j+ωi,j-1+116ωi+1,j+1+ωi+1,j-1+ωi-1,j+1+ωi-1,j-1;in three-dimensional space:(37)ωi,j,k⟵18ωi,j,k+116ωi-1,j,k+ωi+1,j,k+ωi,j-1,k+ωi,j+1,k+ωi,j,k-1+ωi,j,k+1+132ωi-1,j-1,k+ωi-1,j+1,k+ωi+1,j-1,k+ωi+1,j+1,k+132ωi-1,j,k-1+ωi-1,j,k+1+ωi+1,j,k-1+ωi+1,j,k+1+132ωi,j-1,k-1+ωi,j-1,k+1+ωi,j+1,k-1+ωi,j+1,k+1+164ωi-1,j-1,k-1+ωi+1,j-1,k-1+ωi-1,j+1,k-1+ωi-1,j-1,k+1+ωi-1,j+1,k+1+ωi+1,j-1,k+1+ωi+1,j+1,k-1+ωi+1,j+1,k+1.Remark 3.
Because discontinuities may initially exist on boundaries or move to boundaries at a later time, redistribution of boundary points should be made in order to improve the quality of the solution near boundary. For convenience, our attention is restricted to the case in which the physical domainΩp is straight boundary. Assume that a new set of grid points xj1,j2,…,jNdκ+1 is obtained in Ωp by solving the moving mesh equation (26). Then the speeds of the internal grid points xj1,j2,…,jNdκ are given by(38)cj1,j2,…,jNdj1,cj1,j2,…,jNdj2,…,cj1,j2,…,jNdjNd=xj1,j2,…,jNdκ+1-xj1,j2,…,jNdκ,where 1≤jr≤Nxr, 1≤r≤Nd. We assume that the points of boundaries are moving with the same speed as the tangential component of the speed for the internal points adjacent to those boundary point. For an example, when the boundary surface is the boundary surface of x0,j2,…,jNd, we have (39)c0,j2,…,jNdj1,c0,j2,…,jNdj2,…,c0,j2,…,jNdjNd=0,cj1,j2,…,jNdj2,…,cj1,j2,…,jNdjNd,where 1≤jr≤Nxr, 1≤r≤Nd. Thus new boundary points x0,j2,…,jNdκ+1 are defined by (40)x0,j2,…,jNdκ+1=x0,j2,…,jNdκ+c0,j2,…,jNdj1,c0,j2,…,jNdj2,…,c0,j2,…,jNdjNd,where 1≤jr≤Nxr, 1≤r≤Nd. The redistribution for other boundaries can be carried out in a similar way. Numerical experiments show that the above procedure for moving the boundary points is useful in improving the solution resolution.


## 3.1. Time Evolution of Governing Equations

The physical domain for computation isΩp. Given a partition of the mesh domain Kj¯∣Kj¯∈Ωp,Kj¯=xj,xj+1,⋃j¯Kj¯=Ωp and a partition of the time interval 0,T. The average on the cell Kj¯ is(7)w¯j¯=1Kj¯∫Kj¯wj¯dσ.Here, the sign Kj¯ denotes the length for the region Kj¯Kj¯ in one dimension, the area for the region Kj¯Ki¯,j¯ in two dimensions, and the volumes for the region Kj¯Ki¯,j¯,k¯ in three dimensions. Integrating (2) over Kj¯, we have(8)∫Kj¯∂w∂tdσ+∫Kj¯∇xFwdσ=∫Kj¯swdσ.Applying the Green-Gauss’s theorem and rewriting, we have(9)Kj¯dwj¯dt+∮∂Kj¯Fw·nj¯dσ=Kj¯swj¯,where nj¯ is the outward unit normal vector of boundary external surface ∂Kj¯. The Lax-Friedrichs flux is defined by(10)Hu,v,n=12Fu·n+Fv·n-av-u,where a=maxu,vF′u·n. It satisfies the conservation and consistency(11)Hu,v,n=-Hu,v,-n,Hu,u,n=Fu·n.Let ej¯l∣ej¯l∈∂Kj¯,⋃l=1Nlej¯l=∂Kj¯ be a partition of boundary external surface ∂Kj¯ and nj¯ be the outward unit normal vector of surface ej¯l. Then we have a semidiscrete scheme of (2)(12)Kj¯dwj¯dt=-∑l=1Nl∫ej¯lHwj¯intl,wj¯extl,nj¯ldσ+Kj¯swj¯,where wj¯intl or wj¯extl, l=1,2,…,Nl, is the internal or external approximate value at ej¯l.For the one-dimensional space,Nl=2, and ej¯1ej¯2 is the left (right) point of line cell Kj¯. Thus, (12) becomes(13)dwj¯dt=1Kj¯-Hwj¯int1,wj¯ext1+Hwj¯int2,wj¯ext2+swj¯.By the initial reconstruction technique [15] to reset wj¯intl, wj¯extl, l=1,2, on the edge ej¯l, we can obtain the second-order accurate spatial discretization:(14)wj¯int1=wj¯+12ψwj¯-wj¯+1,wj¯-1-wj¯,wj¯ext1=wj¯-1+12ψwj¯-wj¯-1,wj¯-1-wj¯-2,wj¯int2=wj¯+1ext1,wj¯ext2=wj¯+1int1,where ψ is a nonlinear limiter function which is used to suppress the possible pseudo oscillation. In our computations, we use van Leer’s limiter [15](15)ψa,b=signa+signbaba+b+ε.For the two-dimensional space, Nl=4, and the diagram for edges ei¯,j¯l, l=1,2,3,4, of the quadrilateral element Ki¯,j¯ is shown in Figure 1.Figure 1
Schematic diagram of the quadrilateral elementKi¯,j¯.The valueswi¯,j¯intl, wi¯,j¯extl at edge ei¯,j¯l, l=1,2,3,4, can be defined by(16)wi¯,j¯int1=wi¯,j¯+12ψwi¯,j¯-wi¯,j¯+1,wi¯,j¯-1-wi¯,j¯,wi¯,j¯ext1=wi¯,j¯-1+12ψwi¯,j¯-wi¯,j¯-1,wi¯,j¯-1-wi¯,j¯-2,wi¯,j¯ext3=wi¯,j¯+1int1,wi¯,j¯int3=wi¯,j¯+1ext1wi¯,j¯int2=wi¯,j¯+12ψwi¯,j¯-wi¯-1,j¯,wi¯+1,j¯-wi¯,j¯,wi¯,j¯ext2=wi¯+1,j¯+12ψwi¯+1,j¯-wi¯+2,j¯,wi¯,j¯-wi¯+1,j¯,wi¯,j¯ext4=wi¯-1,j¯int2,wi¯,j¯int4=wi¯-1,j¯ext2.For the three-dimensional hexahedron element, the situation is complex. Due to the movement of the mesh points, the hexahedron element will change to polyhedron element. For calculating the element volume and the boundary external surface area, each surface of hexahedron element is divided into two surfaces and calculated. The control element Ki¯,j¯,k¯ is illustrated in Figure 2. The boundary surfaces ei¯,j¯,k¯l, l=1,2,…,12, are defined in Table 1.Table 1
Definition of boundary surfacesei¯,j¯,k¯l, l=1,2,…,12.
Surfaces
Vertexes
ei-,j-,k-1
x
i
+
1
,
j
,
k
+
1
,
x
i
+
1
,
j
,
k
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-2
x
i
+
1
,
j
,
k
,
x
i
+
1
,
j
+
1
,
k
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-3
x
i
,
j
,
k
,
x
i
,
j
,
k
+
1
,
x
i
,
j
+
1
,
k
+
1
ei-,j-,k-4
x
i
,
j
,
k
,
x
i
,
j
+
1
,
k
+
1
,
x
i
,
j
+
1
,
k
ei-,j-,k-5
x
i
,
j
+
1
,
k
,
x
i
,
j
+
1
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-6
x
i
,
j
+
1
,
k
,
x
i
+
1
,
j
+
1
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
ei-,j-,k-7
x
i
,
j
,
k
+
1
,
x
i
,
j
,
k
,
x
i
+
1
,
j
,
k
+
1
ei-,j-,k-8
x
i
,
j
,
k
,
x
i
+
1
,
j
,
k
,
x
i
+
1
,
j
,
k
+
1
ei-,j-,k-9
x
i
,
j
,
k
+
1
,
x
i
+
1
,
j
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
+
1
ei-,j-,k-10
x
i
,
j
,
k
+
1
,
x
i
+
1
,
j
+
1
,
k
+
1
,
x
i
,
j
+
1
,
k
+
1
ei-,j-,k-11
x
i
+
1
,
j
,
k
,
x
i
,
j
,
k
,
x
i
+
1
,
j
+
1
,
k
ei-,j-,k-12
x
i
,
j
,
k
,
x
i
,
j
+
1
,
k
,
x
i
+
1
,
j
+
1
,
kFigure 2
Schematic diagram of the quadrilateral elementKi¯,j¯,k¯.The second-order valueswi¯,j¯,k¯intl, wi¯,j¯,k¯extl at edge ei¯,j¯,k¯l, l=1,2,…,12, can be obtained by(17)wi¯,j¯,k¯int1=wi¯,j¯,k¯int2=w¯i¯,j¯,k¯+12ψw¯i¯,j¯,k¯-w¯i¯-1,j¯,k¯,w¯i¯+1,j¯,k¯-w¯i¯,j¯,k¯,wi¯,j¯,k¯ext1=wi¯,j¯,k¯ext2=w¯i¯+1,j¯,k¯+12ψw¯i¯+1,j¯,k¯-w¯i¯+2,j¯,k¯,w¯i¯,j¯,k¯-w¯i¯+1,j¯,k¯,wi¯,j¯,k¯int5=wi¯,j¯,k¯int6=w¯i¯,j¯,k¯+12ψw¯i¯,j¯,k¯-w¯i¯,j¯-1,k¯,w¯i¯,j¯+1,k¯-w¯i¯,j¯,k¯,wi¯,j¯,k¯ext5=wi¯,j¯,k¯ext6=w¯i¯,j¯+1,k¯+12ψw¯i¯,j¯+1,k¯-w¯i¯,j¯+2,k¯,w¯i¯,j¯,k¯-w¯i¯,j¯+1,k¯,wi¯,j¯,k¯int9=wi¯,j¯,k¯int10=w¯i¯,j¯,k¯+12ψw¯i¯,j¯,k¯-w¯i¯,j¯,k¯-1,w¯i¯,j¯,k¯+1-w¯i¯,j¯,k¯,wi¯,j¯,k¯ext9=wi¯,j¯,k¯ext10=w¯i¯,j¯,k¯+1+12ψw¯i¯,j¯,k¯+1-w¯i¯,j¯,k¯+2,w¯i¯,j¯,k¯-w¯i¯,j¯,k¯+1,wi¯,j¯,k¯int3=wi¯-1,j¯,k¯ext1,wi¯,j¯,k¯int4=wi¯-1,j¯,k¯ext2,wi¯,j¯,k¯ext3=wi¯-1,j¯,k¯int1,wi¯,j¯,k¯ext4=wi¯-1,j¯,k¯int2,wi¯,j¯,k¯int7=wi¯,j¯-1,k¯ext5,wi¯,j¯,k¯int8=wi¯,j¯-1,k¯ext6,wi¯,j¯,k¯ext7=wi¯,j¯-1,k¯int5,wi¯,j¯,k¯ext8=wi¯,j¯-1,k¯int6,wi¯,j¯,k¯int9=wi¯,j¯,k¯-1ext11,wi¯,j¯,k¯int9=wi¯,j¯,k¯-1ext11,wi¯,j¯,k¯ext10=wi¯,j¯,k¯-1int12,wi¯,j¯,k¯ext10=wi¯,j¯,k¯-1int12.Time discretization can be achieved by the strong stability preserving high order Runge-Kutta time discretization [16, 17]. The semidiscrete scheme (12) can be written as(18)w~t=Lw~.Here, the second-order TVD Runge-Kutta method for (18) in the time discretization is(19)w~1=w~n+ΔtLw~n,w~n+1=12w~n+12w~1+12ΔtLw~1.


## 3.2. Pseudo Arc-Length Moving Mesh Scheme

Letx=x1,x2,…,xNd and ξ=ξ1,ξ2,…,ξNd denote the physical and computational coordinates, respectively. Here, Nd denotes the number of spatial dimensions. A one-to-one coordinate transformation from the computational domain Ωc to the physical domain Ωp is denoted by(20)x:ξ⟶x,Ωc⟶Ωp.In the variational approach, the mesh map between the domains Ωc and Ωp is usually provided by(21)Ψξ=12∑r=1Nd∫Ωp∇xξrTGr-1∇xξrdx,where ∇x≔∂x1,∂x2,…,∂xNdT and Gr, r=1,2,…,Nd, are given symmetric positive definite matrices called monitor functions, depending on the underlying solution to be adaptive. The Euler-Lagrange equations of the functional (21) have the form(22)∇x·Gr-1∇xξr=0,1≤r≤Nd.In practice, Ωp may have a very complex geometry, and as a result it is not realistic to solve (22) directly. An alternative is to consider a functional(23)Ψ~x=12∑r=1Nd∫Ωc∇ξxrTGr∇ξxrdξ,where ∇ξ≔∂ξ1,∂ξ2,…,∂ξNdT. In addition, one of the choices of monitor functions is G=ωI, where I is the identity matrix and ω is a positive weight function. For the purpose to have more accuracy near the nonsmooth part of solutions, we introduce the monitor function of pseudo arc-length norm [18, 19](24)ω=1+α1W+α2∇W2,where α1 and α2 are some nonnegative constants. Thus, the corresponding Euler-Lagrange equations about (23) are(25)∇ξ·ω∇ξxr=0,1≤r≤Nd.In our computations, we will use the Gauss-Seidel type iteration method to solve the mesh equation (25)(26)∑r=1Ndϑj1,j2,…,jr+1/2,…,jNdrxj1,j2,…,jr+1,…,jNdκ-xj1,j2,…,jr,…,jNdκ+1-ϑj1,j2,…,jr-1/2,…,jNdrxj1,j2,…,jr,…,jNdκ+1-xj1,j2,…,jr-1,…,jNdκ=0,for 1≤jr≤Nξr, 1≤r≤Nd, and κ=0,1,…, where(27)ϑj1,j2,…,jr-1/2,…,jNdr=ωj1,j2,…,jr-1/2,…,jNd=12Nd-1∑sr=jr-1/2,sk=jk-1/2,jk+1/2ωs1,s2,…,sr,…,sNd,ϑj1,j2,…,jr+1/2,…,jNdr=ωj1,j2,…,jr+1/2,…,jNd=12Nd-1∑sr=jr+1/2,sk=jk-1/2,jk+1/2ωs1,s2,…,sr,…,sNd,where ωj¯1,j¯2,…,j¯Nd≔ωwj¯1,j¯2,…,j¯Ndκ. In addition, let us consider the solutions updating on new grids xκ+1 from the old grids xκ. Here, the solution updating should preserve conservative property for wj¯1,j¯2,…,j¯Nd in the sense of that(28)∑j¯1,j¯2,…,j¯NdKj¯1,j¯2,…,j¯Ndκwj¯1,j¯2,…,j¯Ndκ=∑j¯1,j¯2,…,j¯NdKj¯1,j¯2,…,j¯Ndκ+1wj¯1,j¯2,…,j¯Ndκ+1.Let Dl denote the region scanned by the boundary ej¯1,j¯2,…,j¯Ndl after one iterative step of (28). So we can have the following conservative interpolation scheme:(29)Kj¯1,j¯2,…,j¯Ndκ+1wj¯1,j¯2,…,j¯Ndκ+1=Kj¯1,j¯2,…,j¯Ndκwj¯1,j¯2,…,j¯Ndκ+∑l=1NlF⌢lwj¯1,j¯2,…,j¯Ndintl,wj¯1,j¯2,…,j¯Ndextl,where wj¯1,j¯2,…,j¯Ndintl and wj¯1,j¯2,…,j¯Ndextl are the reconstructed left and right states on the boundary ej¯1,j¯2,…,j¯Ndl by wj¯1,j¯2,…,j¯Ndκ (for more details see Section 3.2). F⌢l·,·, l=1,2,…,Nl, denotes the integration of w over the domain Dl. Here we have omitted the subscripts j¯r, r=1,2,…,Nd, of F⌢l·,· and Dl for simplicity. Since only cell averages of wj¯1,j¯2,…,j¯Nd over Kj¯1,j¯2,…,j¯Nd are known, we should give an approximation of F⌢l·,· instead of its exact value. Then F⌢l·,· can be approximately calculated as(30)F⌢lu,v=max⁡Dl,0·v+min⁡Dl,0·u.Next, we will show the examples to calculate Dl, l=1,2,…,Nd, according to different dimensions.In one-dimensional space,D1, D2 can be obtained from the following:(31)D1=xjκ-xjκ+1,D2=xj+1κ+1-xj+1κ.In two-dimensional space, the changed area Dl, l=1,2,…,4, are illustrated in Figure 3. As an example, we compute D1. Noticing xi,j=xi,j,yi,j, we have(32)D1≔12xi+1,jκ+1-xi,jκyi+1,jκ-yi,jκ+1-yi+1,jκ+1-yi,jκxi+1,jκ-xi,jκ+1.Obviously, D1 is the area, which is positive if xi-1,j-1z+1, xi,j-1z+1, xi,j-1z, and xi-1,j-1z are located in a counter-clockwise order; otherwise it is negative if those four points are located in a clockwise order.Figure 3
Schematic diagram of the two-dimensional changed areasDl, l=1,2,3,4.In the three-dimensional space, the control elementKi+1/2,j+1/2,k+1/2 is illustrated in Figure 4. Similarly, we compute D1 as an example. The changed area D1 is illustrated in Figure 5. D1 is composed by three tetrahedrons, so we have(33)D1=Vxi+1,j,k+1κxi+1,j,k+1κ+1xi+1,j+1,k+1κ+1xi,j,k+1κ+1+Vxi+1,j,k+1κxi+1,j+1,k+1κ+1xi,j,k+1κxi,j,k+1κ+1+Vxi+1,j,k+1κxi+1,j+1,k+1κ+1xi+1,j+1,k+1κxi,j,k+1κ.Each tetrahedron volume can be calculated by the volume formula(34)VOABC=16xOAyOAzOAxOByOBzOBxOCyOCzOC.Obviously, VOABC is the volume, which is positive if A, B, and C are located in a counter-clockwise order; otherwise it is negative if those three points are located in a clockwise order.Figure 4
Schematic diagram of the control elementKi+1/2,j+1/2,k+1/2.Figure 5
Three-dimensional changed areaD1.The solution procedure can be illustrated by the following flowchart.Algorithm 1.
Step  1. Give a uniform mesh ξj1,j2,…,jNd0 on the computational domain Ωc. Initialize the physical mesh xj1,j2,…,jNd0 on the physical domain Ωp, and compute the grid initial values wj¯1,j¯2,…,j¯Nd0. Then, let xj1,j2,…,jNd0≔xj1,j2,…,jNd0 and wj¯1,j¯2,…,j¯Nd0≔wj¯1,j¯2,…,j¯Nd0.
Step  2. With scheme (26), move the grid from xj1,j2,…,jNdκ to xj1,j2,…,jNdκ+1 and the corresponding grid values from wj¯1,j¯2,…,j¯Ndκ to wj¯1,j¯2,…,j¯Ndκ+1 based on scheme (29) for κ=0,1,2,…. Repeat the updating procedure for a fixed number of iterations or until xj1,j2,…,jNdκ+1-xj1,j2,…,jNdκ≤ε. Then, let x~j1,j2,…,jNd=xj1,j2,…,jNdκ+1 and w~j¯1,j¯2,…,j¯Nd=wj¯1,j¯2,…,j¯Ndκ+1 at time level tn.
Step  3. Evolve the underlying PDEs using finite volume scheme (12) and time discretization scheme (19) on the new mesh x~j1,j2,…,jNd to obtain the numerical approximations wj¯1,j¯2,…,j¯Ndn+1 at the time level tn+1.
Step  4. If tn+1⩽T, then let xj1,j2,…,jNd0=x~j1,j2,…,jNd and wj¯1,j¯2,…,j¯Nd0=wj¯1,j¯2,…,j¯Ndn+1, and go to Step  2.Remark 2.
In practice, it is common to use some temporal or spatial smoothing on the weight functionω to obtain smoother meshes in order to avoid very singular mesh and/or large approximation errors near those regions where the solution has a large gradient. A low pass filter should be used to smooth the weight function 2~3 times for each κ: in one-dimensional space:(35)ωj=14ωj-1+2ωj+ωj+1;in two-dimensional space:(36)ωi,j⟵14ωi,j+18ωi+1,j+ωi,j+1+ωi-1,j+ωi,j-1+116ωi+1,j+1+ωi+1,j-1+ωi-1,j+1+ωi-1,j-1;in three-dimensional space:(37)ωi,j,k⟵18ωi,j,k+116ωi-1,j,k+ωi+1,j,k+ωi,j-1,k+ωi,j+1,k+ωi,j,k-1+ωi,j,k+1+132ωi-1,j-1,k+ωi-1,j+1,k+ωi+1,j-1,k+ωi+1,j+1,k+132ωi-1,j,k-1+ωi-1,j,k+1+ωi+1,j,k-1+ωi+1,j,k+1+132ωi,j-1,k-1+ωi,j-1,k+1+ωi,j+1,k-1+ωi,j+1,k+1+164ωi-1,j-1,k-1+ωi+1,j-1,k-1+ωi-1,j+1,k-1+ωi-1,j-1,k+1+ωi-1,j+1,k+1+ωi+1,j-1,k+1+ωi+1,j+1,k-1+ωi+1,j+1,k+1.Remark 3.
Because discontinuities may initially exist on boundaries or move to boundaries at a later time, redistribution of boundary points should be made in order to improve the quality of the solution near boundary. For convenience, our attention is restricted to the case in which the physical domainΩp is straight boundary. Assume that a new set of grid points xj1,j2,…,jNdκ+1 is obtained in Ωp by solving the moving mesh equation (26). Then the speeds of the internal grid points xj1,j2,…,jNdκ are given by(38)cj1,j2,…,jNdj1,cj1,j2,…,jNdj2,…,cj1,j2,…,jNdjNd=xj1,j2,…,jNdκ+1-xj1,j2,…,jNdκ,where 1≤jr≤Nxr, 1≤r≤Nd. We assume that the points of boundaries are moving with the same speed as the tangential component of the speed for the internal points adjacent to those boundary point. For an example, when the boundary surface is the boundary surface of x0,j2,…,jNd, we have (39)c0,j2,…,jNdj1,c0,j2,…,jNdj2,…,c0,j2,…,jNdjNd=0,cj1,j2,…,jNdj2,…,cj1,j2,…,jNdjNd,where 1≤jr≤Nxr, 1≤r≤Nd. Thus new boundary points x0,j2,…,jNdκ+1 are defined by (40)x0,j2,…,jNdκ+1=x0,j2,…,jNdκ+c0,j2,…,jNdj1,c0,j2,…,jNdj2,…,c0,j2,…,jNdjNd,where 1≤jr≤Nxr, 1≤r≤Nd. The redistribution for other boundaries can be carried out in a similar way. Numerical experiments show that the above procedure for moving the boundary points is useful in improving the solution resolution.


## 4. Parallelization Strategy

The perfect parallel strategy can reach the conflicting goals of portability and high performance. Here, we adopt the software architecture of MPI, which is a de facto standard for parallel programming on distributed memory systems [20]. The primary importance with the MPI model is its discrete memory view in programming, which makes it hard to write and often involves a very long development cycle [21, 22]. Carefully thought-out strategies for partitioning the data and managing the resultant communication are essential for a good MPI program. Secondly, due to the distributed nature of the model, global data may be duplicated for performance reasons, resulting in an increase in the overall memory requirement. Besides, because of machine instabilities and the lack of fault tolerance in the model, it is very challenging to get very large MPI jobs running on large parallel systems (although this is true in general for other programming models). Considering our numerical schemes, the schematic diagram for our parallel computation is given in Figure 6. Here, N corresponds to the number of processors. Next, the illustrations aiming at the critical problems for the schematic diagram will be discussed.Figure 6
Schematic diagram of the parallel computation.


### 4.1. Partition the Spatial Mesh Domain

In the parallel computation, the whole domain is divided into some blocks, and each block is distributed to different processors. The distribution is such that each processor gets allocated with one block. Because the computational complexity is related to the number of mesh points, the whole spatial domain is divided according to mesh points rather than spatial coordinates. In addition, the number of mesh points for each processor computation is the same as far as possible. The common partitions for whole spatial domain are given in Figure7. Figure 7(a) is line partition which fits the stripe domain. The surface partition is shown in Figure 7(b) which can be used for the surface domain and body domain. The body partition in Figure 7(c) is used for the three-dimensional spatial domain. Next, we will discuss relationship between the whole domain and processors, and pseudocodes will be given. It is assumed that CellNum is the number of cells in the whole spatial domain. The array Cell_to_P(CellNum) is created to store which processor is assigned to each cell, and array P_CellNum(N) is created to store the number of cells for each processor. The new number Cell_in_P(i) in processor for ith cell can be obtained by the following pseudocode:P_CellNum(1~N)=0doi=1,CellNumprocessor=Cell_to_P(i)P_CellNum(processor)=P_CellNum(processor)+1Cell_in_P(i)=P_CellNum(processor)enddoFigure 7
Partitions for spatial domain. (a) Line partition; (b) surface partition; (c) body partition.
(a)
(b)
(c)In addition, the arraySub_to_entire(max⁡(P_CellNum(N)),N) is created to store the relationship between the whole domain and processors. The pseudocode isdoi=1,CellNumprocessor=Cell_to_P(i)sub_i=Cell_in_P(i)Sub_to_entire(sub_i,processor)=ienddo


### 4.2. Reduction and Synchronization

The size of spatial mesh is different for each processor, which leads to the time step for each processor being different. Thus time step is needed to synchronize in the computation.reduction and broadcast are used to achieve this. Let ht′ be the time step of each processor computation with CFL condition. ht is the time step after synchronization. The codes for MPI Fortran areCALLMPI_Reduceht′,ht,1,MPI_DOUBLE_PRECISION,&MPI_MIN,0,COMM,IERRCALLMPI_Bcastht,1,MPI_DOUBLE_PRECISION,0,COMM,IERR


### 4.3. Processor Communication

For data communication, it is necessary to know neighbor processors for each processor. For line partition, there are two neighbors for current processor to communicate, and if the current processor is on the boundary or its neighbor does not exist, the neighbor processor will be signednull. For the adjacent two processors PA and PB (PA is left; PB is right), there are Ne_j adjacent cells for processors PA and PB. Thus, processor PA needs Ne_j cells to communicate with processor PB. In the same way, processor PB needs N_aj cells to communicate with processor PA. The following pseudocode is given to show the communication between PA and PB.C_A/C_B=1doi=1,P_CellNum(A)ifrighcellbelongtoP_BthenP_A_send(C_A)=iP_B_receive(C_A)=iC_A=C_A+1endifenddodoi=1,P_CellNum(B)ifleftcellbelongtoP_AthenP_B_send(C_B)=iP_A_receive(C_B)=iC_B=C_B+1endifenddoThe arraysP_A_send(N_aj) and P_A_receive(N_aj) are the cells which PA send and receive, while The arrays P_B_send(N_aj) and P_B_receive(N_aj) are the cells which PB send and receive. For the surface partition, there are eight neighbor processors for the current processor, which is shown in Figure 8(a) (here, Cen is the current processor). Because neighbors Nei, i=1,2,3,4, have the common communicating edge with the current processor Cen, we can use the similar way as the line partition to communication. The position of processor Nei, i=5,6,7,8, can be obtained with Nei, i=1,2,3,4. For example, the following pseudocode is to show the communication between Ne5 and Cen.C_Cen/C_Ne5=1doi=1,P_CellNum(Ne5)ifrighcellbelongtoNe1andabovecellbelongtoNe4thenNe5_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddodoi=1,P_CellNum(Cen)ifleftcellbelongtoNe4andbelowcellbelongtoNe1thenCen_send(C_Ne5)=iNe5_receive(C_Ne5)=iC_Ne5=C_Ne5+1endifenddoFigure 8
Processor communication for surface partition; (a) internal communication; (b) boundary communication.
(a)
(b)When the current processor is on the boundary, which is shown in Figure8(b), the ghost cells Cbc can be obtained with the current processor Cen, other ghost cells Cda and Cdb for Cen can be obtained by communicating with Nea and Neb. The following pseudocode is for Cen obtaining CdaC_Cen=1doi=1,P_CellNum(Nea)ifrighcellbelongtonullandbelowcellbelongtoCenthenNea_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoFor body partition, there are26 neighbor processors around the current processor Cen, which is shown in Figure 8(a). We have divided the neighbors into three kinds Nmi, Nlj, and Nvk,i=1,2,…,6, j=1,2,…,12, and k=1,2,…,8, where they are corresponding to three kinds of communication which is shown in Figure 9(b). Nm1, Nl1, and Nv1 are the examples, of which we will give the pseudocode to show the communication of body partition.Figure 9
Processor communication for body partition; (a) neighbor processors; (b) communication blocks.
(a)
(b)ForNm1,C_Cen=1doi=1,P_CellNum(A)ifbehindcellbelongtoCenthenNm1_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoForNl1,C_Cen=1doi=1,P_CellNum(A)ifbehindcellbelongtoNm3andrightcellbelongtoNm1thenNl1_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoForNv1,C_Cen=1doi=1,P_CellNum(A)if(belowcellbelongtoNl1andbehindcellbelongtoNl8&andrightcellbelongtoNl5thenNv1_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoFigure10 shows the communications when the current processor is the boundary processor. Figure 10(a) is the surface boundary while Figure 10(b) is the edge boundary. For Figure 10(a), the communication is the same as the situations in Figure 8, while Figure 10(b) is the similar situation as the line partition. Thus, we leave it out.Figure 10
Boundary processor communication; (a) boundary surface processors; (b) boundary edge processors.
(a)
(b)


## 4.1. Partition the Spatial Mesh Domain

In the parallel computation, the whole domain is divided into some blocks, and each block is distributed to different processors. The distribution is such that each processor gets allocated with one block. Because the computational complexity is related to the number of mesh points, the whole spatial domain is divided according to mesh points rather than spatial coordinates. In addition, the number of mesh points for each processor computation is the same as far as possible. The common partitions for whole spatial domain are given in Figure7. Figure 7(a) is line partition which fits the stripe domain. The surface partition is shown in Figure 7(b) which can be used for the surface domain and body domain. The body partition in Figure 7(c) is used for the three-dimensional spatial domain. Next, we will discuss relationship between the whole domain and processors, and pseudocodes will be given. It is assumed that CellNum is the number of cells in the whole spatial domain. The array Cell_to_P(CellNum) is created to store which processor is assigned to each cell, and array P_CellNum(N) is created to store the number of cells for each processor. The new number Cell_in_P(i) in processor for ith cell can be obtained by the following pseudocode:P_CellNum(1~N)=0doi=1,CellNumprocessor=Cell_to_P(i)P_CellNum(processor)=P_CellNum(processor)+1Cell_in_P(i)=P_CellNum(processor)enddoFigure 7
Partitions for spatial domain. (a) Line partition; (b) surface partition; (c) body partition.
(a)
(b)
(c)In addition, the arraySub_to_entire(max⁡(P_CellNum(N)),N) is created to store the relationship between the whole domain and processors. The pseudocode isdoi=1,CellNumprocessor=Cell_to_P(i)sub_i=Cell_in_P(i)Sub_to_entire(sub_i,processor)=ienddo


## 4.2. Reduction and Synchronization

The size of spatial mesh is different for each processor, which leads to the time step for each processor being different. Thus time step is needed to synchronize in the computation.reduction and broadcast are used to achieve this. Let ht′ be the time step of each processor computation with CFL condition. ht is the time step after synchronization. The codes for MPI Fortran areCALLMPI_Reduceht′,ht,1,MPI_DOUBLE_PRECISION,&MPI_MIN,0,COMM,IERRCALLMPI_Bcastht,1,MPI_DOUBLE_PRECISION,0,COMM,IERR


## 4.3. Processor Communication

For data communication, it is necessary to know neighbor processors for each processor. For line partition, there are two neighbors for current processor to communicate, and if the current processor is on the boundary or its neighbor does not exist, the neighbor processor will be signednull. For the adjacent two processors PA and PB (PA is left; PB is right), there are Ne_j adjacent cells for processors PA and PB. Thus, processor PA needs Ne_j cells to communicate with processor PB. In the same way, processor PB needs N_aj cells to communicate with processor PA. The following pseudocode is given to show the communication between PA and PB.C_A/C_B=1doi=1,P_CellNum(A)ifrighcellbelongtoP_BthenP_A_send(C_A)=iP_B_receive(C_A)=iC_A=C_A+1endifenddodoi=1,P_CellNum(B)ifleftcellbelongtoP_AthenP_B_send(C_B)=iP_A_receive(C_B)=iC_B=C_B+1endifenddoThe arraysP_A_send(N_aj) and P_A_receive(N_aj) are the cells which PA send and receive, while The arrays P_B_send(N_aj) and P_B_receive(N_aj) are the cells which PB send and receive. For the surface partition, there are eight neighbor processors for the current processor, which is shown in Figure 8(a) (here, Cen is the current processor). Because neighbors Nei, i=1,2,3,4, have the common communicating edge with the current processor Cen, we can use the similar way as the line partition to communication. The position of processor Nei, i=5,6,7,8, can be obtained with Nei, i=1,2,3,4. For example, the following pseudocode is to show the communication between Ne5 and Cen.C_Cen/C_Ne5=1doi=1,P_CellNum(Ne5)ifrighcellbelongtoNe1andabovecellbelongtoNe4thenNe5_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddodoi=1,P_CellNum(Cen)ifleftcellbelongtoNe4andbelowcellbelongtoNe1thenCen_send(C_Ne5)=iNe5_receive(C_Ne5)=iC_Ne5=C_Ne5+1endifenddoFigure 8
Processor communication for surface partition; (a) internal communication; (b) boundary communication.
(a)
(b)When the current processor is on the boundary, which is shown in Figure8(b), the ghost cells Cbc can be obtained with the current processor Cen, other ghost cells Cda and Cdb for Cen can be obtained by communicating with Nea and Neb. The following pseudocode is for Cen obtaining CdaC_Cen=1doi=1,P_CellNum(Nea)ifrighcellbelongtonullandbelowcellbelongtoCenthenNea_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoFor body partition, there are26 neighbor processors around the current processor Cen, which is shown in Figure 8(a). We have divided the neighbors into three kinds Nmi, Nlj, and Nvk,i=1,2,…,6, j=1,2,…,12, and k=1,2,…,8, where they are corresponding to three kinds of communication which is shown in Figure 9(b). Nm1, Nl1, and Nv1 are the examples, of which we will give the pseudocode to show the communication of body partition.Figure 9
Processor communication for body partition; (a) neighbor processors; (b) communication blocks.
(a)
(b)ForNm1,C_Cen=1doi=1,P_CellNum(A)ifbehindcellbelongtoCenthenNm1_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoForNl1,C_Cen=1doi=1,P_CellNum(A)ifbehindcellbelongtoNm3andrightcellbelongtoNm1thenNl1_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoForNv1,C_Cen=1doi=1,P_CellNum(A)if(belowcellbelongtoNl1andbehindcellbelongtoNl8&andrightcellbelongtoNl5thenNv1_send(C_Cen)=iCen_receive(C_Cen)=iC_Cen=C_Cen+1endifenddoFigure10 shows the communications when the current processor is the boundary processor. Figure 10(a) is the surface boundary while Figure 10(b) is the edge boundary. For Figure 10(a), the communication is the same as the situations in Figure 8, while Figure 10(b) is the similar situation as the line partition. Thus, we leave it out.Figure 10
Boundary processor communication; (a) boundary surface processors; (b) boundary edge processors.
(a)
(b)


## 5. Numerical Examples

In this section, the numerical convergence and the efficiency for our parallel schemes will be shown. In addition, the practical problem will be considered. In our simulation, the parameters are taken asγ=1.2, Eα=1.7, Eβ=0.35, Q=52.6, kα=287.0, and kβ=0.0359. The simulations are taken on the computer group of four nodes. There are two Intel E5620CPUs on each node, and each CPU is quad-core.


### 5.1. Example  1

Firstly, the numerical accuracy is tested for the one-dimensional Euler equations. A periodic boundary condition is used and the initial conditions are set to beρx,0=1.0+0.2sin⁡x, ux,0=0.5, px,0=1.0, αx,0=0.5+0.2sin⁡x, and βx,0=0.5+0.2sin⁡x. The exact solution for this problem is(41)ρx,t=1.0+0.2sin⁡x-0.5t,ux,t=0.5,px,t=1.0,αx,t=0.5+0.2sin⁡x-0.5t,βx,t=0.5+0.2sin⁡x-0.5t.The final time for this computation is T=2.0. The monitor function for this computation is(42)M=1+α1sξ,where s=ργ/p and α1=0.1. The parameter γ is 1.2. The computational domain is taken as -π,π. The errors and convergence orders in the density are obtained and compared in Table 2. From Table 2, we can see that our scheme is the second-order convergent scheme. In addition, Sod’s classical shock tube problem is studied to show that our schemes can capture the singularities and avoid the nonphysical phenomena. The Riemann initial values are(43)ρ,u,p=1.0,0.0,1.0if0≤x<0.5,0.125,0.0,0.1if0.5≤x<1.0and reflecting boundary conditions at x=0 and x=1. The terminal time is T=0.2. The solutions about density and the trajectory of mesh are shown in Figure 11. There are 100 cells used in the simulation. The monitor function employed for this computation is(44)M=1+α1ρξmaxξρξ+α2sξmaxξsξ.It is found that the mesh is adaptive with the solutions and our schemes can capture the discontinuous jump in numerical simulations. The computational complexity and errors analysis are compared between the fixed mesh and moving mesh schemes in Table 3. Comparing (F1) and (M1), we can use 50 mesh points to reach the effect of 100 fixed mesh points while costing the same CPU time. In addition, the different arc-length parameters with 100 mesh points are simulated to compare 200 and 250 fixed mesh points schemes. Our schemes can reach better results with less mesh points. Particularly, L2, L∞ errors are obviously reduced. Thus we can conclude that our schemes can capture the singularities and avoid the nonphenomena.Table 2
Errors and convergence orders in the density.
Cells
Fixed mesh
Moving mesh
Error
Order
Error
Order
20
3.4932E-02
—
3.5100E-02
—
40
1.0219E-02
1.7733
1.0244E-02
1.7767
80
2.6508E-03
1.9468
2.6380E-03
1.9572
160
6.3429E-04
2.0632
6.2496E-04
2.0776
320
1.4987E-04
2.0815
1.4569E-04
2.1009Table 3
Comparison of various fixed and moving mesh.
Description
CPU time
Time steps
L
1 error
L
2 error
L
∞ error
Fixed mesh (cells):
(F1)100
1.56
E
-
02
48
7.9754
E
-
03
1.5122
E
-
02
6.4432
E
-
02
(F2)200
6.24
E
-
02
97
4.2999
E
-
03
1.0988
E
-
02
7.1633
E
-
02
(F3)250
9.36
E
-
02
121
3.2221
E
-
03
8.5794
E
-
03
6.7484
E
-
02
Moving mesh
(
c
e
l
l
s
,
α
1
,
α
2
)
:
(M1)50,20,50
1.56
E
-
02
69
7.8921
E
-
03
1.4041
E
-
02
6.4284
E
-
02
(M2)100,20,100
4.68
E
-
02
129
4.2209
E
-
03
1.0347
E
-
02
7.1121
E
-
02
(M3)100,10,10
7.80
E
-
02
236
3.1650
E
-
03
7.9560
E
-
03
6.5914
E
-
02Figure 11
The simulation results; (a) solution about density; (b) trajectory of mesh.
(a)
(b)


### 5.2. Example  2

This is the two-dimensional example. Firstly, the numerical accuracy is tested with the periodic boundary condition and initial conditionsρx,0=1.0+0.2sin⁡x+y, ux,0=vx,0=0.5, px,0=1.0, αx,0=0.5+0.2sin⁡x+y, and βx,0=0.5+0.2sin⁡x+y. The source terms are adjusted so that the governing equations (2) can fit the exact solution(45)ρx,t=1.0+0.2sin⁡x+y-t,ux,t=vx,t=0.5,px,t=1.0,αx,t=0.5+0.2sin⁡x+y-t,βx,t=0.5+0.2sin⁡x+y-t.Here the monitor function is about entropy (42). The final time for the computation is T=2.0. The errors and convergence orders about density are shown in Table 4, which implies that our schemes are the second-order convergence. Next, we apply our schemes to a two-dimensional explosion problem. The computational domain is illustrated in Figure 12. The boundary condition is reflective, and the initial conditions in the unreacted region are ρ,u,v,w,α,β=1,0,0,0,1,1. The initial conditions in the reacted region are obtained by the ZND model [23]. There are 240×240 cells in our computation. The monitor function(46)M=1+α1ρ2+α2∇ρ2 is used in the computation. The terminal time is 0.5. The parallel cost with different processors is compared in Tables 5 and 6. Table 5 is the comparison of parallel efficiency with different processors for line partition, while Table 6 is for the surface partition. The Speedup and Efficiency for parallel computations are defined by(47)SpeedupwithNprocessors=CosttimewithoneprocessorCosttimewithNprocessors,EfficiencywithNprocessors=SpeedupwithNprocessorsN.Both tables illustrate that our parallel schemes are effective and the computational cost time is reduced. Figure 13 is the computational results at time 0.15, and Figure 14 is at time 0.3. The solutions figures show that our schemes can capture the shock waves and the mesh trajectory is adaptive.Table 4
Errors and convergence orders in the density.
Cells
Fixed mesh
Moving mesh
Error
Order
Error
Order
20
2.4050
E
-
01
—
2.4980
E
-
01
—
40
6.6037
E
-
02
1.8647
6.6661
E
-
02
1.9058
80
1.6605
E
-
02
1.9916
1.6590
E
-
02
2.0056
160
3.9552
E
-
03
2.0698
3.9543
E
-
03
2.0688
320
9.2949
E
-
04
2.0892
9.2871
E
-
04
2.0901Table 5
Parallel efficiency with different processors for line partition.
Number of processors
CPU time
Speedup
Efficiency
11×1
15778
1.0
1.0
22×1
7894
1.9987
0.9994
33×1
5281
2.9877
0.9959
44×1
3976
3.9683
0.9921
55×1
3354
4.7042
0.9408
66×1
2891
5.4576
0.9096Table 6
Parallel efficiency with different processors for surface partition.
Number of processors
CPU time
Speedup
Efficiency
11×1
15778
1.0
1.0
42×2
3965
3.9793
0.9948
93×3
2053
7.6853
0.8539
164×4
1620
9.7395
0.6087Figure 12
The computational domain.Figure 13
The computational results at time 0.15; (a) solution about density; (b) trajectory of mesh.
(a)
(b)Figure 14
The computational results at time 0.3; (a) solution about pressure; (b) trajectory of mesh.
(a)
(b)


### 5.3. Example  3

The last example is the three-dimensional problem. The numerical exact test is taken firstly. The periodic boundary condition is used and the initial conditions are set to beρx,0=1.0+0.2sin⁡x+y+z, ux,0=vx,0=wx,0=0.5, px,0=1.0, αx,0=0.5+0.2sin⁡x+y+z, and βx,0=0.5+0.2sin⁡x+y+z. The terminal time is T=2.0 and the monitor function is about entropy (42). The exact solution is(48)ρx,t=1.0+0.2sin⁡x+y+z-1.5t,ux,t=vx,t=wx,t=0.5,px,t=1.0,αx,t=0.5+0.2sin⁡x+y+z-1.5t,βx,t=0.5+0.2sin⁡x+y+z-1.5t. The errors and convergence orders about density are shown in Table 7, which implies that our schemes are the second-order convergence. The computational domain and boundary conditions are shown in Figure 15. The initial values for reacted region and unreacted region are the same as Example two. Table 8 is used to show the parallel efficiency for body partition, which illustrates that the computational time is reduced and our parallel scheme can deal with the three-dimensional problems. There are 120×120×120 cells in our simulations. The monitor function(49)M=1+α1ρ2+α2∇ρ2 is used in the computation. The terminal time is 0.5. Figures 16 and 17 show the adaptive solutions. Figure 16 shows the three-dimensional sliced effect of the adaptive mesh with the solutions at different time points. Figure 17(a) is the simulated result of the block in the center of the domain about density at time 0.15, while Figure 17(b) is a x-direction of the sliced simulation about density at time 0.3. The spatial cells are distorted to adapt to the shock wave front. They show us that our schemes can capture shock waves.Table 7
Errors and convergence orders in the density.
Cells
Fixed mesh
Moving mesh
Error
Order
Error
Order
30
1.0109
—
1.0138
—
40
5.8547E-01
1.8987
5.8571
E
-
01
1.9071
50
3.8052
E
-
01
1.9310
3.8027
E
-
01
1.9358
60
2.6448
E
-
01
1.9952
2.6421
E
-
01
1.9970
70
1.9246
E
-
01
2.0620
1.9223
E
-
01
2.0633Table 8
Parallel efficiency with different processors for body partition.
Number of processors
CPU time
Speedup
Efficiency
11×1×1
24298
1.0
1.0
82×2×2
3370
7.2101
0.9013
273×3×3
1145
21.2209
0.7859
324×4×2
1053
23.0750
0.7211Figure 15
The computational domain.Figure 16
The simulation solutions; (a) solution about density at time 0.15; (b) solution about pressure at time 0.3.
(a)
(b)Figure 17
The simulation solutions; (a) a block solutions about density at time 0.15; (b) a piece of solutions about density at time 0.3.
(a)
(b)


## 5.1. Example  1

Firstly, the numerical accuracy is tested for the one-dimensional Euler equations. A periodic boundary condition is used and the initial conditions are set to beρx,0=1.0+0.2sin⁡x, ux,0=0.5, px,0=1.0, αx,0=0.5+0.2sin⁡x, and βx,0=0.5+0.2sin⁡x. The exact solution for this problem is(41)ρx,t=1.0+0.2sin⁡x-0.5t,ux,t=0.5,px,t=1.0,αx,t=0.5+0.2sin⁡x-0.5t,βx,t=0.5+0.2sin⁡x-0.5t.The final time for this computation is T=2.0. The monitor function for this computation is(42)M=1+α1sξ,where s=ργ/p and α1=0.1. The parameter γ is 1.2. The computational domain is taken as -π,π. The errors and convergence orders in the density are obtained and compared in Table 2. From Table 2, we can see that our scheme is the second-order convergent scheme. In addition, Sod’s classical shock tube problem is studied to show that our schemes can capture the singularities and avoid the nonphysical phenomena. The Riemann initial values are(43)ρ,u,p=1.0,0.0,1.0if0≤x<0.5,0.125,0.0,0.1if0.5≤x<1.0and reflecting boundary conditions at x=0 and x=1. The terminal time is T=0.2. The solutions about density and the trajectory of mesh are shown in Figure 11. There are 100 cells used in the simulation. The monitor function employed for this computation is(44)M=1+α1ρξmaxξρξ+α2sξmaxξsξ.It is found that the mesh is adaptive with the solutions and our schemes can capture the discontinuous jump in numerical simulations. The computational complexity and errors analysis are compared between the fixed mesh and moving mesh schemes in Table 3. Comparing (F1) and (M1), we can use 50 mesh points to reach the effect of 100 fixed mesh points while costing the same CPU time. In addition, the different arc-length parameters with 100 mesh points are simulated to compare 200 and 250 fixed mesh points schemes. Our schemes can reach better results with less mesh points. Particularly, L2, L∞ errors are obviously reduced. Thus we can conclude that our schemes can capture the singularities and avoid the nonphenomena.Table 2
Errors and convergence orders in the density.
Cells
Fixed mesh
Moving mesh
Error
Order
Error
Order
20
3.4932E-02
—
3.5100E-02
—
40
1.0219E-02
1.7733
1.0244E-02
1.7767
80
2.6508E-03
1.9468
2.6380E-03
1.9572
160
6.3429E-04
2.0632
6.2496E-04
2.0776
320
1.4987E-04
2.0815
1.4569E-04
2.1009Table 3
Comparison of various fixed and moving mesh.
Description
CPU time
Time steps
L
1 error
L
2 error
L
∞ error
Fixed mesh (cells):
(F1)100
1.56
E
-
02
48
7.9754
E
-
03
1.5122
E
-
02
6.4432
E
-
02
(F2)200
6.24
E
-
02
97
4.2999
E
-
03
1.0988
E
-
02
7.1633
E
-
02
(F3)250
9.36
E
-
02
121
3.2221
E
-
03
8.5794
E
-
03
6.7484
E
-
02
Moving mesh
(
c
e
l
l
s
,
α
1
,
α
2
)
:
(M1)50,20,50
1.56
E
-
02
69
7.8921
E
-
03
1.4041
E
-
02
6.4284
E
-
02
(M2)100,20,100
4.68
E
-
02
129
4.2209
E
-
03
1.0347
E
-
02
7.1121
E
-
02
(M3)100,10,10
7.80
E
-
02
236
3.1650
E
-
03
7.9560
E
-
03
6.5914
E
-
02Figure 11
The simulation results; (a) solution about density; (b) trajectory of mesh.
(a)
(b)


## 5.2. Example  2

This is the two-dimensional example. Firstly, the numerical accuracy is tested with the periodic boundary condition and initial conditionsρx,0=1.0+0.2sin⁡x+y, ux,0=vx,0=0.5, px,0=1.0, αx,0=0.5+0.2sin⁡x+y, and βx,0=0.5+0.2sin⁡x+y. The source terms are adjusted so that the governing equations (2) can fit the exact solution(45)ρx,t=1.0+0.2sin⁡x+y-t,ux,t=vx,t=0.5,px,t=1.0,αx,t=0.5+0.2sin⁡x+y-t,βx,t=0.5+0.2sin⁡x+y-t.Here the monitor function is about entropy (42). The final time for the computation is T=2.0. The errors and convergence orders about density are shown in Table 4, which implies that our schemes are the second-order convergence. Next, we apply our schemes to a two-dimensional explosion problem. The computational domain is illustrated in Figure 12. The boundary condition is reflective, and the initial conditions in the unreacted region are ρ,u,v,w,α,β=1,0,0,0,1,1. The initial conditions in the reacted region are obtained by the ZND model [23]. There are 240×240 cells in our computation. The monitor function(46)M=1+α1ρ2+α2∇ρ2 is used in the computation. The terminal time is 0.5. The parallel cost with different processors is compared in Tables 5 and 6. Table 5 is the comparison of parallel efficiency with different processors for line partition, while Table 6 is for the surface partition. The Speedup and Efficiency for parallel computations are defined by(47)SpeedupwithNprocessors=CosttimewithoneprocessorCosttimewithNprocessors,EfficiencywithNprocessors=SpeedupwithNprocessorsN.Both tables illustrate that our parallel schemes are effective and the computational cost time is reduced. Figure 13 is the computational results at time 0.15, and Figure 14 is at time 0.3. The solutions figures show that our schemes can capture the shock waves and the mesh trajectory is adaptive.Table 4
Errors and convergence orders in the density.
Cells
Fixed mesh
Moving mesh
Error
Order
Error
Order
20
2.4050
E
-
01
—
2.4980
E
-
01
—
40
6.6037
E
-
02
1.8647
6.6661
E
-
02
1.9058
80
1.6605
E
-
02
1.9916
1.6590
E
-
02
2.0056
160
3.9552
E
-
03
2.0698
3.9543
E
-
03
2.0688
320
9.2949
E
-
04
2.0892
9.2871
E
-
04
2.0901Table 5
Parallel efficiency with different processors for line partition.
Number of processors
CPU time
Speedup
Efficiency
11×1
15778
1.0
1.0
22×1
7894
1.9987
0.9994
33×1
5281
2.9877
0.9959
44×1
3976
3.9683
0.9921
55×1
3354
4.7042
0.9408
66×1
2891
5.4576
0.9096Table 6
Parallel efficiency with different processors for surface partition.
Number of processors
CPU time
Speedup
Efficiency
11×1
15778
1.0
1.0
42×2
3965
3.9793
0.9948
93×3
2053
7.6853
0.8539
164×4
1620
9.7395
0.6087Figure 12
The computational domain.Figure 13
The computational results at time 0.15; (a) solution about density; (b) trajectory of mesh.
(a)
(b)Figure 14
The computational results at time 0.3; (a) solution about pressure; (b) trajectory of mesh.
(a)
(b)


## 5.3. Example  3

The last example is the three-dimensional problem. The numerical exact test is taken firstly. The periodic boundary condition is used and the initial conditions are set to beρx,0=1.0+0.2sin⁡x+y+z, ux,0=vx,0=wx,0=0.5, px,0=1.0, αx,0=0.5+0.2sin⁡x+y+z, and βx,0=0.5+0.2sin⁡x+y+z. The terminal time is T=2.0 and the monitor function is about entropy (42). The exact solution is(48)ρx,t=1.0+0.2sin⁡x+y+z-1.5t,ux,t=vx,t=wx,t=0.5,px,t=1.0,αx,t=0.5+0.2sin⁡x+y+z-1.5t,βx,t=0.5+0.2sin⁡x+y+z-1.5t. The errors and convergence orders about density are shown in Table 7, which implies that our schemes are the second-order convergence. The computational domain and boundary conditions are shown in Figure 15. The initial values for reacted region and unreacted region are the same as Example two. Table 8 is used to show the parallel efficiency for body partition, which illustrates that the computational time is reduced and our parallel scheme can deal with the three-dimensional problems. There are 120×120×120 cells in our simulations. The monitor function(49)M=1+α1ρ2+α2∇ρ2 is used in the computation. The terminal time is 0.5. Figures 16 and 17 show the adaptive solutions. Figure 16 shows the three-dimensional sliced effect of the adaptive mesh with the solutions at different time points. Figure 17(a) is the simulated result of the block in the center of the domain about density at time 0.15, while Figure 17(b) is a x-direction of the sliced simulation about density at time 0.3. The spatial cells are distorted to adapt to the shock wave front. They show us that our schemes can capture shock waves.Table 7
Errors and convergence orders in the density.
Cells
Fixed mesh
Moving mesh
Error
Order
Error
Order
30
1.0109
—
1.0138
—
40
5.8547E-01
1.8987
5.8571
E
-
01
1.9071
50
3.8052
E
-
01
1.9310
3.8027
E
-
01
1.9358
60
2.6448
E
-
01
1.9952
2.6421
E
-
01
1.9970
70
1.9246
E
-
01
2.0620
1.9223
E
-
01
2.0633Table 8
Parallel efficiency with different processors for body partition.
Number of processors
CPU time
Speedup
Efficiency
11×1×1
24298
1.0
1.0
82×2×2
3370
7.2101
0.9013
273×3×3
1145
21.2209
0.7859
324×4×2
1053
23.0750
0.7211Figure 15
The computational domain.Figure 16
The simulation solutions; (a) solution about density at time 0.15; (b) solution about pressure at time 0.3.
(a)
(b)Figure 17
The simulation solutions; (a) a block solutions about density at time 0.15; (b) a piece of solutions about density at time 0.3.
(a)
(b)


## 6. Conclusion and Discussion

In this paper, we have discussed the parallel computation for pseudo arc-length moving mesh schemes. Different from the traditional Euler numerical scheme, the communications between processors are more complex, which include adjacent processor and diagonal processor. The parallel schemes including line, surface, and body partition are all considered in this work and the pseudocodes and schematic diagram are given. Finally, the numerical examples are shown to illustrate that our parallel schemes are effective and the pseudo arc-length moving mesh schemes are convergent and can be used to capture the shock wave.



---

*Source: 5896940-2017-07-12.xml*