Scientific Programming
Volume 2017 (2017), Article ID 5896940, 17 pages
https://doi.org/10.1155/2017/5896940

Research Article

Parallel Pseudo Arc-Length Moving Mesh Schemes for Multidimensional Detonation

Jianguo Ning, Xinpeng Yuan, Tianbao Ma, and Jian Li

State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China

Correspondence should be addressed to Tianbao Ma; madabal@bit.edu.cn

Received 17 January 2017; Accepted 16 May 2017; Published 12 July 2017

Academic Editor: Piotr Luszczek

Copyright © 2017 Jianguo Ning et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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].

where is the mass density, the pressure, the temperature, the gas constant, the heat release parameter, and the constants of reaction rates, and and 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

where is the multidimensional vector, is the multidimensional matrix function, and is the chemical reaction source term. For the one-dimensional space,

For the two-dimensional space,

In the case of three-dimensional space,

where , , and are the Cartesian components of the fluid velocity in the , , and directions, respectively. Total energy density is defined as

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 . Given a partition of the mesh domain and a partition of the time interval . The average on the cell is

Here, the sign denotes the length for the region in one dimension, the area for the region in two dimensions, and the volumes for the region in three dimensions. Integrating (2) over , we have

Applying the Green-Gauss’s theorem and rewriting, we have

where is the outward unit normal vector of boundary external surface . The Lax-Friedrichs flux is defined by

where . It satisfies the conservation and consistency

Let be a partition of boundary external surface and be the outward unit normal vector of surface . Then we have a semidiscrete scheme of (2)

where or , , is the internal or external approximate value at .

For the one-dimensional space, , and is the left (right) point of line cell . Thus, (12) becomes

By the initial reconstruction technique [15] to reset , , , on the edge , we can obtain the second-order accurate spatial discretization:

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]

For the two-dimensional space, , and the diagram for edges , , of the quadrilateral element is shown in Figure 1.

Figure 1: Schematic diagram of the quadrilateral element .

The values , at edge , , can be defined by

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 is illustrated in Figure 2. The boundary surfaces , , are defined in Table 1.

Table 1: Definition of boundary surfaces , .

Surfaces
Vertexes

Figure 2: Schematic diagram of the quadrilateral element .

The second-order values , at edge , , can be obtained by

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

Here, the second-order TVD Runge-Kutta method for (18) in the time discretization is

3.2. Pseudo Arc-Length Moving Mesh Scheme

Let and denote the physical and computational coordinates, respectively. Here, denotes the number of spatial dimensions. A one-to-one coordinate transformation from the computational domain to the physical domain is denoted by

In the variational approach, the mesh map between the domains and is usually provided by

where and , , 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

In practice, 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

where . In addition, one of the choices of monitor functions is , where 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]

where and are some nonnegative constants. Thus, the corresponding Euler-Lagrange equations about (23) are

In our computations, we will use the Gauss-Seidel type iteration method to solve the mesh equation (25)

for , , and , where

where In addition, let us consider the solutions updating on new grids from the old grids . Here, the solution updating should preserve conservative property for in the sense of that

Let denote the region scanned by the boundary after one iterative step of (28). So we can have the following conservative interpolation scheme:

where and are the reconstructed left and right states on the boundary by (for more details see Section 3.2). , , denotes the integration of over the domain . Here we have omitted the subscripts , , of and for simplicity. Since only cell averages of over are known, we should give an approximation of instead of its exact value. Then can be approximately calculated as

Next, we will show the examples to calculate , , according to different dimensions.

In one-dimensional space, , can be obtained from the following:

In two-dimensional space, the changed area , , are illustrated in Figure 3. As an example, we compute . Noticing , we have

Obviously, is the area, which is positive if , , , and 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 areas , .

In the three-dimensional space, the control element is illustrated in Figure 4. Similarly, we compute as an example. The changed area is illustrated in Figure 5. is composed by three tetrahedrons, so we have

Each tetrahedron volume can be calculated by the volume formula

Obviously, is the volume, which is positive if , , and 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 element .

Figure 5: Three-dimensional changed area .

The solution procedure can be illustrated by the following flowchart.

Algorithm 1.
Step  1. Give a uniform mesh on the computational domain . Initialize the physical mesh on the physical domain , and compute the grid initial values . Then, let and .
Step  2. With scheme (26), move the grid from to and the corresponding grid values from to based on scheme (29) for . Repeat the updating procedure for a fixed number of iterations or until . Then, let and at time level .
Step  3. Evolve the underlying PDEs using finite volume scheme (12) and time discretization scheme (19) on the new mesh to obtain the numerical approximations at the time level .
Step  4. If , then let and , 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:

in two-dimensional space:

in three-dimensional space:

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 is straight boundary. Assume that a new set of grid points is obtained in by solving the moving mesh equation (26). Then the speeds of the internal grid points are given by

where , . 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 , we have

where , . Thus new boundary points are defined by

where , . 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, 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 Figure 7. 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 is the number of cells in the whole spatial domain. The array is created to store which processor is assigned to each cell, and array is created to store the number of cells for each processor. The new number in processor for th cell can be obtained by the following pseudocode:

 1~N                        

(a)

(b)

(c)

Figure 7: Partitions for spatial domain. (a) Line partition; (b) surface partition; (c) body partition.

In addition, the array is created to store the relationship between the whole domain and processors. The pseudocode is

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. and are used to achieve this. Let be the time step of each processor computation with CFL condition. is the time step after synchronization. The codes for MPI Fortran are

 _  _ 

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 signed . For the adjacent two processors and ( is left; is right), there are adjacent cells for processors and . Thus, processor needs cells to communicate with processor . In the same way, processor needs cells to communicate with processor . The following pseudocode is given to show the communication between and .

The arrays and are the cells which send and receive, while The arrays and are the cells which send and receive. For the surface partition, there are eight neighbor processors for the current processor, which is shown in Figure 8(a) (here, is the current processor). Because neighbors , , have the common communicating edge with the current processor , we can use the similar way as the line partition to communication. The position of processor , , can be obtained with , . For example, the following pseudocode is to show the communication between and .

(a)

(b)

Figure 8: Processor communication for surface partition; (a) internal communication; (b) boundary communication.

When the current processor is on the boundary, which is shown in Figure 8(b), the ghost cells can be obtained with the current processor , other ghost cells and for can be obtained by communicating with and . The following pseudocode is for obtaining

For body partition, there are neighbor processors around the current processor , which is shown in Figure 8(a). We have divided the neighbors into three kinds , , and , , and , where they are corresponding to three kinds of communication which is shown in Figure 9(b). , , and are the examples, of which we will give the pseudocode to show the communication of body partition.

(a)

(b)

Figure 9: Processor communication for body partition; (a) neighbor processors; (b) communication blocks.

For ,

For ,

For ,

Figure 10 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.

(a)

(b)

Figure 10: Boundary processor communication; (a) boundary surface processors; (b) boundary edge processors.

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 , , , , , and . 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 , , , , and . The exact solution for this problem is

The final time for this computation is . The monitor function for this computation is

where and . The parameter is . 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

and reflecting boundary conditions at and . The terminal time is . The solutions about density and the trajectory of mesh are shown in Figure 11. There are cells used in the simulation. The monitor function employed for this computation is

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 mesh points to reach the effect of fixed mesh points while costing the same CPU time. In addition, the different arc-length parameters with mesh points are simulated to compare and fixed mesh points schemes. Our schemes can reach better results with less mesh points. Particularly, , 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

—

—

40

80

160

320

Table 3: Comparison of various fixed and moving mesh.

Description
CPU time
Time steps
error
error
error

Fixed mesh (cells):

(F1)

(F2)

(F3)

Moving mesh

(M1)

(M2)

(M3)

(a)

(b)

Figure 11: The simulation results; (a) solution about density; (b) trajectory of mesh.

5.2. Example  2

This is the two-dimensional example. Firstly, the numerical accuracy is tested with the periodic boundary condition and initial conditions , , , , and . The source terms are adjusted so that the governing equations (2) can fit the exact solution

Here the monitor function is about entropy (42). The final time for the computation is . 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 . The initial conditions in the reacted region are obtained by the model [23]. There are cells in our computation. The monitor function

is used in the computation. The terminal time is . 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 and for parallel computations are defined by

Both tables illustrate that our parallel schemes are effective and the computational cost time is reduced. Figure 13 is the computational results at time , and Figure 14 is at time . 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

—

—

40

80

160

320

Table 5: Parallel efficiency with different processors for line partition.

Number of processors
CPU time
Speedup
Efficiency

Table 6: Parallel efficiency with different processors for surface partition.

Number of processors
CPU time
Speedup
Efficiency

Figure 12: The computational domain.

(a)

(b)

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.

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 , , , , and . The terminal time is and the monitor function is about entropy (42). The exact solution is

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 cells in our simulations. The monitor function

is used in the computation. The terminal time is . 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 , while Figure 17(b) is a -direction of the sliced simulation about density at time . 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

—

—

40

50

60

70

Table 8: Parallel efficiency with different processors for body partition.

Number of processors
CPU time
Speedup
Efficiency

Figure 15: The computational domain.

(a)

(b)

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.

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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant nos. 11390363 and 11532012).

References

F. Davis, Detonation, Univ. California Press, 1979.

G. Wang, D. Zhang, K. Liu, and J. Wang, “An improved CE/SE scheme for numerical simulation of gaseous and two-phase detonations,” Computers and Fluids. An International Journal, vol. 39, no. 1, pp. 168–177, 2010.

A. Harten and J. M. Hyman, “Self-adjusting grid methods for one-dimensional hyperbolic conservation laws,” Journal of Computational Physics, vol. 50, no. 2, pp. 235–269, 1983.

J. M. Stockie, J. A. Mackenzie, and R. . Russell, “A moving mesh method for one-dimensional hyperbolic conservation laws,” SIAM Journal on Scientific Computing, vol. 22, no. 5, pp. 1791–1813, 2000.

X. Yang, W. Huang, and J. Qiu, “A moving mesh WENO method for one-dimensional conservation laws,” SIAM Journal on Scientific Computing, vol. 34, no. 4, pp. A2317–A2343, 2012.

H. Tang and T. Tang, “Adaptive mesh methods for one- and two-dimensional hyperbolic conservation laws,” SIAM Journal on Numerical Analysis, vol. 41, no. 2, pp. 487–515, 2003.

J. Han and H. Tang, “An adaptive moving mesh method for two-dimensional ideal magnetohydrodynamics,” Journal of Computational Physics, vol. 220, no. 2, pp. 791–812, 2007.

P. He and H. Tang, “An adaptive moving mesh method for two-dimensional relativistic magnetohydrodynamics,” Computers and Fluids. An International Journal, vol. 60, pp. 1–20, 2012.

J. Ning, X. Yuan, T. Ma, and C. Wang, “Positivity-preserving moving mesh scheme for two-step reaction model in two dimensions,” Computers and Fluids. An International Journal, vol. 123, pp. 72–86, 2015.

M. G. Knepley and D. A. Karpeev, “Mesh algorithms for PDE with sieve I: mesh distribution,” Scientific Programming, vol. 17, no. 3, pp. 215–230, 2009.

K. Morris, D. W. I. Rouson, M. N. Lemaster, and S. Filippone, “Exploring capabilities within ForTrilinos by solving the 3D Burgers equation,” Scientific Programming, vol. 20, no. 3, pp. 275–292, 2012.

M. Mubarak, S. Seol, Q. Lu, and M. S. Shephard, “A parallel ghosting algorithm for the flexible distributed mesh database,” Scientific Programming, vol. 21, no. 1-2, pp. 17–42, 2013.

C. Wang, X. Dong, and C.-W. Shu, “Parallel adaptive mesh refinement method based on WENO finite difference scheme for the simulation of multi-dimensional detonation,” Journal of Computational Physics, vol. 298, pp. 161–175, 2015.

S. Taki and T. Fujiwara, “Numerical simulation on the establishment of gaseous detonation,” Prog Astronaut Aeronaut, vol. 9, pp. 186–200, 1984.

B. van Leer, “Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method,” Journal of Computational Physics, vol. 32, no. 1, pp. 101–136, 1979.

C.-W. Shu and S. Osher, “Efficient implementation of essentially nonoscillatory shock-capturing schemes,” Journal of Computational Physics, vol. 77, no. 2, pp. 439–471, 1988.

S. Gottlieb and C.-W. Shu, “Total variation diminishing Runge-Kutta schemes,” Mathematics of Computation, vol. 67, no. 221, pp. 73–85, 1998.

K. Chen, “Error equidistribution and mesh adaptation,” SIAM Journal on Scientific Computing, vol. 15, no. 4, pp. 798–818, 1994.

X. Wang, T.-B. Ma, and J.-G. Ning, “A pseudo arc-length method for numerical simulation of shock waves,” Chinese Physics Letters, vol. 31, no. 3, Article ID 030201, 2014.

H. Jin, D. Jespersen, P. Mehrotra, R. Biswas, L. Huang, and B. Chapman, “High performance computing using MPI and OpenMP on multi-core parallel systems,” Parallel Computing, vol. 37, no. 9, pp. 562–575, 2011.

L. Hochstein and V. R. Basili, “The ASC-alliance projects: a case study of large-scale parallel scientific code development,” Computer, vol. 41, no. 3, pp. 50–58, 2008.

L. Hochstein, F. Shull, and L. B. Reid, “The role of MPI in development time: A case study,” in Proceedings of the 2008 SC - International Conference for High Performance Computing, Networking, Storage and Analysis, SC 2008, usa, November 2008.

W. Fickett and W. W. Wood, “Flow calculations for pulsating one-dimensional detonations,” Physics of Fluids, vol. 9, no. 5, pp. 903–916, 1966.