A Riemann-Solver Free Spacetime Discontinuous Galerkin Method for General Conservation Laws

This paper summarizes a Riemann-solver-free spacetime discontinuous Galerkin method developed for general conservation laws. The method integrates the best features of the spacetime Conservation Element/Solution Element (CE/SE) method and the discontinuous Galerkin (DG) method. The core idea is to construct a staggered spacetime mesh through alternate cell-centered CEs and vertex-centered CEs within each time step. Inside each SE, the solution is approximated using high-order spacetime DG basis polynomials. The spacetime flux conservation is enforced inside each CE using the DG concept. The unknowns are stored at both vertices and cell centroids of the spatial mesh. However, the solutions at vertices and cell centroids are updated at different time levels within each time step in an alternate fashion. Thanks to the staggered spacetime formulation, there are no left and right states for the solution at the spacetime interface. Instead, the solution available to evaluate the flux is continuous across the interface. Therefore, no (approximate) Riemann solvers are needed to provide a unique numerical flux. The current method can be used to solve arbitrary conservation laws including the compressible Euler equations, shallow water equations and magnetohydrodynamics (MHD) equations without the need of any form of Riemann solvers. A set of benchmark problems of various conservation laws are presented to demonstrate the accuracy of the method.


Introduction
Traditionally, when solving conservation laws numerically, (approximate) Riemann solvers are used to provide unique interface fluxes in solvers based on the finite volume method, discontinuous Galerkin method and some other methods.Numerical methods based on Riemann solvers have achieved tremendous success in solving hyperbolic systems, e.g., compressible Euler equations, where the eigenstructure of the system is clearly known, in the past several decades.However, when physical processes get more complicated, for example, magnetohydrodynamics (MHD), the success of Riemann-solvers is far from satisfactory.For example, the Roe scheme [1] has been modified to solve the MHD system.Roe's scheme requires eigen-decomposition and becomes very complicated in MHD equations.In addition, the Roe averages of quantities are not clearly defined in MHD system.Moreover, due to the complexity and non-strict hyperbolicity of the MHD system, the validity of the eigensystems of the MHD system is not unanimously agreed upon among researchers.Therefore, Riemann-solver free approaches are highly desirable to avoid the uncertainties caused by imperfect Riemann solvers.
In recent years, the author has been developing a Riemann-solver-free high order spacetime method for arbitrary conservation laws [2]- [12].The method is inspired by the spacetime Conservation Element/Solution Element (CE/SE) [13] method and the discontinuous Galerkin (DG) [14] method.The method integrates the best features of the two methods, i.e. the Riemann-solver-free feature of the CE/SE method and the high-order accuracy of the DG method.The core idea of the method is to construct a staggered spacetime mesh through alternate cell-centered CEs and vertex-centered CEs (cf. Figure 1 (right)) within each time step.Note that CEs at the vertex level are defined via the dual mesh obtained by connecting the surrounding cell centers and face centers.The difference between SEs and CEs is that the SE includes the volumeless vertical spike.Inside each SE (cf. Figure 1 (left)), the solution is approximated using high-order spacetime DG basis polynomials.The spacetime flux conservation is enforced inside each CE using the DG discretization.The unknowns are stored at both vertices and cell centroids of the spatial mesh.However, the solutions at vertices and cell centroids are updated at different time levels within each time step in an alternate fashion.Due to the alternate cell-vertex solution updating strategy and its DG ingredient, the method has been termed as the discontinuous Galerkin cellvertex scheme (DG-CVS) [4].
The method offers benefits to both the spacetime CE/SE method and the spacetime DG method.The benefits for the spacetime CE/SE method is two fold.First, DG-CVS increases the CE/SE method's temporal and spatial accuracies by employing high-degree polynomial basis functions inside the SE.Second, DG-CVS presents universal definitions of CEs and SEs for arbitrary meshes.The CE and SE definitions in DG-CVS are independent of the underlying spatial mesh and the same definitions can be easily extended from 1-D (cf. Figure 1) to higher-dimensions (cf. Figure 2 and Figure 3) without any ambiguity.
The benefits for the spacetime DG method are also two-fold.First, DG-CVS provides a Riemann-solver-free approach for the DG methods.Second, DG-CVS does not need special treatment for the diffusion terms and thus is conceptually simpler than existing traditional DG methods for diffusion equations.
This paper summarizes the author's work on DG-CVS in recent years.This paper is organized as follows.Section 2 describes the current method in detail where the two ingredients of the method, spacetime DG method and the alternate cell-vertex solution updating strategy, will be explained.In Section 3, many numerical examples are presented to demonstrate the accuracy of the method in solving various conservation laws including scalar advection equations, compressible Euler equations, shallow water equations and MHD equations.Finally, conclusions will be summarized in Section 4.

Riemann-Solver Free Spacetime DG Method
To illustrate the Riemann-solver free discontinuous Galerkin (DG) formulation, we consider the following general one-dimensional time dependent partial differential equation system ( ) that is used to describe many physical processes.Here, u is the solution vector containing the physical quantities to be determined.f is the spatial flux vector.Without loss of generality, a source term r is also in- cluded in Equation (1).Note that f may contain both advective and diffusive fluxes.Both f and r can be functions of u or ∇u or both.Appropriate initial and boundary conditions must also be provided to solve Equation (1).

Spacetime Discontinuous Galerkin Formulation
Following the idea of the DG method, an approximate solution h u is sought within each spacetime solution element (SE), denoted as K .When restricted to the SE, h u belongs to the finite dimensional space ( ) . Assume all the components of the conservative variables inside each SE are approximated by polynomials of the same degree, i.e. where are some type of spacetime polynomial basis functions defined within the solution element, are the unknowns to be determined and N is the number of basis functions depending on the degree of the polynomial function.Using such spacetime DG solution space, the solution approximation can be arbitrarily high order accurate in both space and time.
Following the Galerkin orthogonality principle, multiply (1) with each of the basis functions i φ ( )  and integrate over a spacetime CE to obtain the following weak form d 0, 1,2, , where K Ω is the spacetime conservation element (CE) associated with the solution element K.Note that the conservation element is identical to the solution element except for the volumeless vertical spike in the solution element as shown in Figure 1.The spacetime flux conservation in weak form as in (3) is for each individual spacetime conservation element.Therefore, the current method can be considered as a spacetime discontinuous Galerkin method.

Integrating (3) by parts results in d d
where ( ) is the spacetime flux tensor and ( ) is the outward unit normal of the CE boundary, i.e.K Γ = ∂Ω , of the spacetime conservation element (CE) under consideration.Note that the partial integration is also performed on the time-dependent term, which is a salient difference between spacetime DG methods and semi-discrete DG methods.As can be seen, the formulation in (4) contains both the volume integral and the surface integral.

Cell-Vertex Solution Updating Strategy
The current method inherits the core idea of the CE/SE method using a staggered spacetime mesh to enforce the spacetime flux conservation.However, the construction of the staggered spacetime slabs in the current method deviates from the CE/SE method.In the current method, unknowns are stored at both vertices and cell centroids of the spatial mesh, and the solutions at vertices and cell centroids are updated at different time levels within each time step.At the beginning of each physical time step, the solution is assumed known at the cell centroids of the mesh, either given as the initial condition or obtained from the previous time step.Inside each new time step, the solution is updated in two successive steps.The first step updates the solution at vertices at the halftime level ( ) based on the known cell-centroid solutions at the previous time level ( ) n t .The second step updates the solution at cell centroids at the new time level ( ) t + based on the known vertex solutions at the previous half-time level ( ) 1 2 n t + .The same process is repeated for new time steps.The solution updating at the cell level or the vertex level is based on the key equation (4).To evaluate the integrals in (4), first divide the conservation element into the following portions: • The interior volume K Ω where the solution is associated with the spacetime node at the new time level.
where the left hand side contains the unknowns and the right hand side contains the known information.
Since the top and the bottom surface of the CE are horizontal, which leads to

⋅ = H n u on the top face and
u on the bottom face, Equation ( 6) can be simplified to Equation ( 7) is a nonlinear equation system for each spacetime node which can be solved by the Newton-Raphson iterative method.
To further illustrate the idea of enforcing the spacetime flux conservation, consider the conservation element shown in Figure 4. Suppose the solution at the spacetime node . The interior volume of the conservation element is also associated with where the solution is to be found.
Due to the cell-vertex solution updating strategy and its DG ingredient, the current method has been termed as the DG-CVS method [4].
Conventional DG methods use some type of Riemann solvers to provide numerical fluxes across element interfaces.Thanks to the above explained staggered space-time formulation, there is no "left" and "right" states for the solution at the spacetime interface.Instead, the solution available to evaluate the flux is continuous across the interface.Therefore, no (approximate) Riemann solvers are needed to provide a unique numerical flux.The Riemann-solver-free feature is highly desirable for systems where the eigenstructure of the underlying system is not readily known and for the purpose of avoiding some pathological behaviors (e.g., carbuncle problem) associated with some Riemann solvers.Note that in the CE/SE method [13], the time derivative of the solution is obtained from the spatial derivative of the solution using the original governing equation.By contrast, in the current method, the time derivative of the solution is treated as an independent unknown and handled in the same way as that for the spatial derivative of the solution.Therefore, even for the same second-order version, the current method deviates from the CE/SE method.Figure 5 compares the solutions obtained from the CE/SE method and the current method for the sinusoidal wave advection problem.It can be seen that the current method exhibits invisible phase error after 100 periods, which is an outstanding feature in contrast to the CE/SE solution where significant phase error can be observed.

Quadrature-Free Implementation
If the underlying spatial mesh is unstructured (cf. Figure 5), the vertex-level CEs are general polygonal cylinders containing general polygonal bases where the Gaussian quadrature rule cannot be directly applied.Therefore, it is desired that the volume integral in Equation ( 7) is evaluated using a quadrature-free approach.In addition, quadrature-free approaches help to significantly improve the efficiency of the method.
To allow the quadrature-free implementation, the flux and source vectors must be expanded in terms of the basis polynomials similar to those used to expand the conservative variables, namely, for two-dimensional cases , , and .
where N  is the number of basis polynomials of one degree higher than those to expand the conservative variables as in Equation ( 2).The requirement of using basis polynomials of degree 1 p + for flux and source term expansions is necessary to ensure the accuracy of the scheme.
The method based on the Cauchy-Kovalewski (CK) procedure [15] [16] is especially suitable for this purpose.In the current method, the conservative variables are expanded using the Taylor polynomials.With the Taylor polynomials, the spatial and temporal derivatives of the solution are explicitly available.Using the CK procedure, the spacetime derivatives of the flux vectors can be obtained using the spacetime derivatives of the conservative variables.
With Taylor polynomials, the integrand of the volume integral has the following general form for twodimensional cases ( ) which results from the products between Taylor polynomials.In (9), , l m and n are integer exponents.The idea is to use the divergence theorem explained in [17] to reduce the volume integral to surface integrals and further to line integrals.The detailed procedure has been described in [4] and will not be repeated here.

Numerical Examples
In this section, various numerical cases are presented to demonstrate that current Riemann-solver free DG method is able to solve many different types of conservation laws without using any Riemann-solvers.

Linear Scalar Advection Equation
We first consider the following linear scalar advection equation: with the initial condition given as follows [18]: , , e and , , max 1 , 0 It is clear that this problem is about the advection of mixed Gaussian, square, triangle and elliptical waves.Periodic boundary conditions are assumed at the ends of the domain.
The computational domain is discretized by 200 evenly spaced cells which is a pretty coarse mesh for this case.The time step used in all simulations is 0.0025 t δ = corresponding to the Courant number 0.25 σ = . To resolve discontinuities, the limiting procedure described in [3] is used.Figure 6 shows the limited solutions at 8.0 t = using 1 p (second order) and 3 p (fourth order) basis polynomials, respectively.As can be seen, both the discontinuities and smooth regions are better captured with higher p .With increasing p , the transition of discontinuities is steeper.The second order run is not able to resolve the Gaussian wave which is steep and smooth.Another feature of the solution which is worthy of mentioning is that the wave profiles are symmetric, which indicates the phase error of the current method is very small.

Burgers Equation
Here, the following 2-D viscous Burgers equation is solved.
with the analytical solution given as x y = is a constant location.The 1-D version of this case was presented in [19].This case is constructed such that the original wave is propagated without changing shape under the effect of both nonlinear advection and linear diffusion.The initial In this test, a medium wave corresponding to 1 = /ν δx is simulated using the present method.p solution and the 1 p solution.As can be seen, the wave has been captured accurately in terms of both location and magnitude.Figure 7 also compares the computed solution gradients ( x u or y u ) with the exact solution.As can be seen, the 3 p solution is superior to the 1 p solution in resolving the local sharp extremum of the solution gradient.This case demonstrate that the current method is able to solve diffusion problems without any special tricks for diffusion terms as required by many other DG methods.Note that the treatment of diffusion terms in traditional semi-discrete DG methods are non-trivial because of the so-called "variational crime".More test cases on diffusion problems have been presented in [7].

Level Set Equation
The level set method introduced by Osher and Sethian [20] has been widely used to solve interface problems emerging in applications such as the two-phase fluid dynamics, image segmentation, computer vision and computational geometry etc.The biggest advantages of the level set method is that it naturally allows the interface to have topological changes due to the interface breaking and merging.
In the level set method, the interface is represented by the zero level set ( ) of a level set function ( ) ,t ψ x .ψ is often taken to be the signed distance function to the interface.The level set equation can be expressed as where ( ) is the velocity field which evolves the interface.Equation ( 13) can be rewritten in conservative form ( ) Here, the DG-CVS method is used to solve Equation ( 14).In the two examples presented here, the com- The interface is represented by the zero level set of the signed distance field.Since the interface evolution often occurs within a narrow band, the initial signed distance field is narrowed by the following , if ; where  is chosen to be 0.1. .800 time steps were finish one revolution.Figure 9 shows the fourth order ( ) Compared to the exact solution, the smearing is very small which indicates the small numerical dissipation of this high order scheme.
Another example is to simulate the deformation of a circular bubble of radius 015 centered at ( ) 0.5, 0.75 .The circle is deformed due to the following velocity field [22]: sin 2π sin π cos π , and sin 2π sin π cos π which are obtained from ( ) where Ψ is the following stream function 1 sin π sin π π x y Ψ = proposed in [22].The time dependent term in (15) was proposed in [23] to control the magnitude and direction of the velocity.At 2 t T = , the velocity starts to reverse direction and at time t T = , the original circle is ex- pected to be restored.The velocity field given by ( 15 Figure 10 shows the fourth order ( ) evolution of the signed distance field and the zero level set at various times.At 0.5 2 t T = = , the bubble starts to restore its shape.At 4.0 t = , the bubble recovers its original circular shape.The comparison with the exact solution shows the accuracy of the current method.Note that reinitialization is not needed with the current high order method in this case.

Compressible Euler Equations
The two-dimensional compressible Euler equations can be written in the following vector form: where 2 2 , , and are the conservative state vector, flux vectors in x -and y -directions, respectively.In (17), ρ , u , v , p , E , and H are density, x -and y -components of the velocity, pressure and the total energy, and the total enthalpy, respectively.Here   ( ) where e is the specific internal energy and γ is the ratio of specific heats.
This first Euler equation case is the classic Sod's shocktube problem.The initial conditions are a highpressure region on the left and a low-pressure region on the right separated by a diaphragm located at 0.5 x = .At time 0 t = , the diaphragm is broken, a shock wave will propagate to the right and an expansion wave pro- pagate to the left.There is a contact discontinuity in between across which the density and entropy are discontinuous while the pressure and velocity are continuous.The spatial domain is discretized into 100 equally spaced elements.Figure 11 shows the comparison of 1 p numerical solution with exact solution at time

t =
. The good agreement in both wave location and sharpness demonstrates the accuracy of curren Riemann-solver free method.
The second case is a supersonic flow with 3 = M flowing through a channel of length 3 and height 1.A step of height 0.2 stands at location 0.6 downstream of the inlet.The step acts as an obstacle to the supersonic inflow.A transient detached shock wave will be formed and hit the top wall.The shock wave will reflect from the top wall and further hit the bottom wall and reflect again.The flow expands rapidly at the singularity corner.In the present simulation, a rectangular mesh of 1 40 is used.The linear basis polynomials ( ) p is employed in this test.Figure 12 shows the density field at 4.0 t = .As can be seen, the main features including the detached bow shock, the Mach stem below the top wall, slip line, reflected shock waves are well captured.Here, no any type of boundary conditions is applied at the singularity corner and no numerical difficulty was encountered.No expansion shocks are seen around the corner.In addition, the present method does not exhibit the so-called carbuncle problem in the subsonic region behind the bow shock.
To illustrate the mesh flexibility of the current method, a simulation about a nearly incompressible subsonic flow with Mach number 0.1 M = passing around a circular cylinder with radius 0.5 R = on an overset Cartesian/quadrilateral mesh is presented.Figure 13 shows the near-body pressure field at the steady state with the background overset mesh also displayed.Indeed, the current method works for arbitrarily unstructured meshes.

Shallow Water Equations
The shallow water equations (SWE) are widely used as the mathematical model to numerically simulate the dam break, river inundation, failure of levees and tide of ocean in coastal and civil engineering.The frictionless shallow water equation can be expressed in the following conservative form ( ) in which the conservative vector, flux vectors and source vector are Here, the total water depth x where ( ) is the free-surface elevation and ( ) d x is the still water depth.u and v are the x -and y -component of the depth-averaged velocity.g is the acceleration due to gravity.0 x S and 0 y S are the bottom slopes in x -and y -directions, respec- tively.If the bottom has zero slope, the source terms in Equation ( 18) would be absent.
The first shallow water case is a two-dimensional shocktube problem which is analogous to the shocktube problem used to verify a compressible Euler solver.The computational domain is a [ ] [ ]   symmetry boundaries are identical to those in the 1-D case.Therefore, the available analytical solutions [24] ca be used to verify the present numerical solution.As can be seen, the shock location and strength have been captured very well.This is in contrast to the results reported in [25] showing that the least-square finite element method is not able to capture the strong shock in the right location.
The second case is another classical benchmark problem widely used to verify shallow water equation solvers [26].Figure 15 (left) depicts the domain geometry.Initially, two reservoirs are separated by a lock gate which is of 75 meters wide.The water levels are 10 meters and 5 meters, respectively.The reservoir on the left has higher  water level.The resolution of the computational mesh is 2.5 meters.At 0 t = , the lock gate is opened.Figure 15 shows the carpet view of the water level at 7.2 t = seconds.The 1 p solution qualitatively compares well with that in [26] where a Riemann-solver-based method is used.

3 3 × MHD Model Equations
The present Riemann-solver-free method is first applied to solve the following 3 3 × MHD model system in the phase space which exactly preserves the MHD singularities [27]: where the longitudinal viscosity µ and the magnetic resistivity η are the dissipative coefficients, and the Hall coefficient χ is a dispersive coefficient.The above quantities are related to the MHD system via is the Alfven wave speed, p is the pressure, ρ is the density and γ is the ratio of specific heats.The model system (19) generates three wave families which are the Alfven wave, the slow and the fast MHD waves.This model system has been investigated by several authors [27] [28] to provide insights to improve Godunov-typed numerical schemes for MHD equations.
The test case presented here is taken from Ref. [28] where conventional and modified Roe's scheme were used to solve the system.In this case, 3 c = and the computational domain [ ] 0,1 is evenly divided into 150 cells.The dissipative term on the right hand side of ( 19) is assumed to be zero.
The initial discontinuity is an entropy-violating shock defined as ( ) . The final time is 0.0376506024 t = . 150 time steps have been run to reach the final time.Figure 16 shows u vs. x and v vs.
x solutions from the second order simulation along with the analytical solutions.The agreement is very well and the overshoots/undershoots around the discontinuity region have been suppressed using a slope limiter.

Ideal Compressible MHD Equations
The ideal compressible MHD equations combines the Euler equations of gas dynamics with Maxwell equations of electromagnetics for problems in which viscous, resistive and relativistic effects can be neglected.The ideal MHD equations contain a set of nonlinear hyperbolic equations as shown in the following conservative form: ( ) represent the mass density, the velocity vector, the total energy, the hydrodynamic pressure, permeability of vacuum and the magnetic field.In addition to satisfying the induction equation (Equation ( 23) above), the magnetic field is also divergence free, i.e.

∇ ⋅ = B
which can be considered as an extra constraint.Note that in Equations ( 20)-( 23), the so-called Powell approach [29] is followed to deal with 0 ∇ ⋅ = B by adding source terms proportional to ∇ ⋅B to the equations.
The first MHD case presented here is known as the Brio-Wu case [30] which is a one-dimensional MHD  .As can be seen, the complex MHD waves including the compound wave, the shock wave, contact discontinuity and rarefaction waves can be captured by the present Riemann-solver free method.However, the unlimited solutions exhibit some non-physical overshoots and undershoots near discontinuities.

Conclusions
This paper summarizes a Riemann-solver free spacetime discontinuous Galerkin method for general hyperbolic conservation laws.Numerical examples demonstrate that the method is able to accurately solve various conservation laws including the scalar advection equations, compressible Euler equations, shallow water equations and magnetohydrodynamics equations without using any Riemann solvers.To summarize, the method has the following features: • Based on space-time formulation.Therefore, the current method automatically satisfies the Geometric Conservation Law (GCL) for moving mesh problems [10] without special care as needed by ALE-based methods.
• High-order accuracy in both space and time for smooth problems.Space and time are handled in a unified way based on space-time flux conservation and high-order space-time discontinuous basis functions.This is in contrast to semi-discrete methods where the temporal order of accuracy is limited by the order of Backward Difference Formula (BDF) or the multi-step Runge-Kutta method.
• Riemann-solver free.DG-CVS does not need (approximate) Riemann solvers to provide numerical fluxes as needed in finite volume or traditional DG methods.The Riemann-solver-free feature is highly desirable for systems where the eigenstructure of the underlying system is not readily known and for the purpose of avoiding some pathological behaviors associated with Riemann solvers.
• Reconstruction free.DG-CVS solves for the solution and its all spatial and temporal derivatives simultaneously at each spacetime node, thus eliminating the need of reconstruction.
• Suitable for arbitrary spatial meshes.The CE and SE definitions in DG-CVS are independent of the underlying spatial mesh and the same definitions can be easily extended from 1-D to 2-D arbitrarily unstructured meshes without any ambiguity.
• Highly compact regardless of order of accuracy.Only information at the immediate neighboring nodes will be needed to update the solution at the new time level.Compactness eases the parallelization of the flow solver.
• Also accurately solves time-dependent diffusion equations.In the current method, the inclusion of diffusion terms is straightforward and simple in exactly the same way how advection terms are handled by simply incorporating the diffusive flux into the space-time flux.No extra reconstruction or recovery or ad hoc penalty and coupling terms are needed to ensure the consistency of the variational form for diffusive terms.
The future work remains to further improve the efficiency of the method and the robustness of the method when high order polynomials are used for nonsmooth problems.

Figure 1 .
Figure 1.Solution elements (SEs) and conservation elements (CEs) in the x t − domain.Left: solution elements and right: conservation elements.

Figure 4 .
Figure 4. Illustration of spacetime flux conservation on a 1-D cell-level CE.

Figure 5 .
Figure 5.Comparison of the solutions of sinusoidal wave advection after 100 periods.Left: CE/SE.right: current method (second-order version).

.
G and E are Gaussian and ellip- tical function, respectively, where

Figure 6 .
Figure 6.Linear advection of mixed waves.Limited solution at 8.0 t = .The computational domain consists of 200 uniform elements.Thick solid line: exact solution; circles and thin colorful solid lines: current numerical solutions; and stars: oscillation indicators.solution at 0 t = and the boundary conditions on all four boundaries are provided by the analytical solution.Therefore the boundary conditions are time dependent.In this test, a medium wave corresponding to 1 = /ν

Figure 7 plots the 1 p
solution along the diagonal line 0 x y − = together with the exact solution.The 3 p solution is not shown since no much visual difference can be seen between the 3

Figure 7 .
Figure 7. Solution of 2-D viscous Burgers equation along the diagonal line 0 x y − = at 20 T = .Left: u , middle: 1 p ) represents a vortex which stretch the circular fluid body into a thin filament toward the vortex center.The simulation is performed on a relatively coarse 40 40 ×

Figure 8 .
Figure 8. Zalesak disk.Initial level set is the signed distance to the disk.

Figure 9 .
Figure 9. Rotation of Zalesak disk.The locations of the disk at several instants within one revolution.The image order is from left to right and from top to bottom.

t
= , the lower left corner [ ] [ ] of the domain is filled with water of height 4 3.0 h = , and the rest of the domain is filled with water of height 1 0.03 h = .Symmetry conditions are assumed on the left and the bottom boundaries of the domain.The domain is discretized into 99 99 × evenly distributed rectangular cells.The simulation is run till 2.0 t = at which the shock is still within the domain.The solutions along the

Figure 12 .
Figure 12.Density field of a supersonic ( ) 3 M = flow through a channel with a forward-facing step at

Figure 13 .
Figure 13.Steady pressure field of the inviscid flow ( ) 0.1 M = passing around a circular cylinder.The

Figure 14
Figure 14 shows the 1 p solutions along the bottom symmetry boundary.The analytical solutions are alsoshown for comparison purpose.As can be seen, the shock location and strength have been captured very well.This is in contrast to the results reported in[25] showing that the least-square finite element method is not able to capture the strong shock in the right location.The second case is another classical benchmark problem widely used to verify shallow water equation solvers[26].Figure15(left) depicts the domain geometry.Initially, two reservoirs are separated by a lock gate which is of 75 meters wide.The water levels are 10 meters and 5 meters, respectively.The reservoir on the left has higher

Figure 15 .
Figure 15.2-D dam break problem.Left: geometry, right: carpet view of the water elevation at 7.2 s t = . with

Figure 16 .
Figure 16.Solutions of the Myong case.

Figure 17
Figure 17 shows the simulated solution at 0.2 t =

Figure 17 .Finally
Figure 17.Solutions of the Brio-Wu case.

Figure 18 .
shows the magnetic field and velocity magnitude field at 0It can seen that the solution fields have been preserved well after one period's advection.