Next Article in Journal
Innovation in an Existing Backpressure Turbine for Ensure Better Sustainability and Flexible Operation
Next Article in Special Issue
Multi-Scale CFD Modeling of Plate Heat Exchangers Including Offset-Strip Fins and Dimple-Type Turbulators for Automotive Applications
Previous Article in Journal
Maintenance Optimization of Offshore Wind Turbines Based on an Opportunistic Maintenance Strategy
Previous Article in Special Issue
CFD Simulation of Defogging Effectivity in Automotive Headlamp
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Simple and Robust Shock-Capturing Approach for Discontinuous Galerkin Discretizations

1
Department of Aeronautics and Astronautics, Stanford University, Stanford, CA 94305, USA
2
Department of Mechanical Engineering, University of Twente, 7522 NB Enschede, The Netherlands
*
Author to whom correspondence should be addressed.
Energies 2019, 12(14), 2651; https://doi.org/10.3390/en12142651
Submission received: 1 June 2019 / Revised: 24 June 2019 / Accepted: 3 July 2019 / Published: 10 July 2019
(This article belongs to the Special Issue Computational Fluid Dynamics (CFD) 2018)

Abstract

:
The discontinuous Galerkin (DG) method has become popular in Computational Fluid Dynamics mainly due to its ability to achieve high-order solution accuracy on arbitrary grids, its high arithmetic intensity (measured as the ratio of the number of floating point operations to memory references), and the use of a local stencil that makes scalable parallel solutions possible. Despite its advantages, several difficulties hinder widespread use of the DG method, especially in industrial applications. One of the major challenges remaining is the capturing of discontinuities in a robust and accurate way. In our previous work, we have proposed a simple shock detector to identify discontinuities within a flow solution. The detector only utilizes local information to sense a shock/discontinuity ensuring that one of the key advantages of DG methods, their data locality, is not lost in transonic and supersonic flows. In this work, we reexamine the shock detector capabilities to distinguish between smooth and discontinuous solutions. Furthermore, we optimize the functional relationships between the shock detector and the filter strength, and present it in detail for others to use. By utilizing the shock detector and the corresponding filtering-strength relationships, one can robustly and accurately capture discontinuities ranging from very weak to strong shocks. Our method is demonstrated in a number of two-dimensional canonical examples.

1. Introduction

Ever since Computational Fluid Dynamics (CFD) began to play an important role in analyzing fluid motion and designing industrial products, engineers have sought techniques which can increase the accuracy of a flow simulation without increasing its associated computational cost. Many CFD researchers have developed a growing list of numerical methods to meet these demands. In particular, CFD developments over the past half century have focused on finite volume, finite difference, and spectral methods [1,2,3,4,5,6]. However, conventional CFD methods have faced difficulties in achieving higher-order (order of solution convergence larger than 2) simulations from both the implementation and high-performance computing perspectives. For example because of their low arithmetic intensity and large computational stencils, higher-order Weighted Essentially Non-Oscillatory (WENO) methods [7] on unstructured meshes and higher-order finite difference methods [8] can lose efficiency when running large-scale simulations.
To overcome these shortcomings in parallel performance and scalability, during the past two decades researchers have focused on finite-element methods such as the discontinuous Galerkin (DG) method. Reed and Hill first introduced the DG method in their seminal paper [9]. Later on, a solid theoretical foundation for the DG method was established by Cockburn and Shu in a series of papers [10,11,12,13,14]. For an overview of the DG method, the reader is referred to [15], which provides many more details than is practical to discuss in this paper. Because of its mathematical formulation, the DG method has high arithmetic intensity compared to the finite volume and the finite difference methods. Furthermore, higher-order accuracy in space can easily be obtained by adding more degrees of freedom in combination with higher-order polynomial bases within a given element. For these reasons, the DG method can obtain higher-order spatial accuracy relatively easily even on complex geometries and when compared to both finite volume and finite difference methods. Additionally, in the DG method only the nearest face-neighbors are needed, while many operations for a given element do not require information from neighboring elements at all, which is an attractive feature for modern computer hardware. Due to these advantages, the DG method has become popular in many fields.
Despite these advantages, the DG method has not yet been fully adopted in industry mainly because of some remaining challenges for practical use. One of these difficulties is the robust and accurate capturing of discontinuities within a flow solution, which is the main objective of the present work. Specifically, and like most other CFD methods, the polynomial nature of the basis functions used by the DG method lead to spurious oscillations in flow simulations if discontinuities emerge during a computation. Because of this fundamental issue, engineers and scientists have pursued methods that can help attain the two fundamental objectives for shock-capturing in the DG framework: the balancing of robustness and accuracy. In order to satisfy these goals, various approaches have been proposed in the literature [16,17,18,19,20,21,22,23,24,25,26]. Researchers have devised multiple methods to sense a discontinuity with information from a local stencil of surrounding neighbors. After the discontinuity is identified, a flow solution can be stabilized with any one several resolving methods including artificial viscosity, limiting, and filtering. A brief comparison of these approaches to shock detection and resolution can be found in [27].
In our previous work, we proposed an approach for discontinuity detection for the DG method and demonstrated its shock-capturing capabilities [27]. Our shock detector robustly distinguishes between smooth solutions and discontinuities even for polynomial orders as low as two and three, which are relatively low orders for such shock-capturing schemes in the DG literature. Furthermore, the shock detector parameter and the necessary filtering strengths were correlated through numerical experiments, enabling the user to avoid tedious parameter tuning processes. These correlations provide robust shock-capturing capabilities regardless of the strengths of shocks present in a flow simulation. However, in that early work, we did not provide a quantitative analysis of the mathematical relationships between the detector and the filtering strengths. At the time, it was not clear whether the preliminary relationships used provided excessive filtering which, despite its solution smoothness advantages, might lead to a severe loss in accuracy. For this reason, in this paper the relationships between the shock detector and the filtering strengths are optimized to obtain maximum accuracy while maintaining robustness of the shock-capturing capability. These relationships are presented for future use.
The remaining sections are organized as follows. The governing equations and the associated numerical discretization approach are introduced in Section 2. In Section 3, the shock detector proposed in [27] and a filtering method are briefly described. Additionally, functional relationships between the shock detector and the filtering method are optimized to obtain accurate solutions while maintaining robustness. Numerical simulation results are shown in Section 4. Finally, we close the paper with concluding remarks and ongoing work in Section 5.

2. Governing Equations and Numerical Discretization

This section briefly describes the governing equations and numerical discretization used in this paper. The majority of this section is restated, with some modifications to enhance clarity, from [27].

2.1. Governing Equations

In this work, we consider the Euler equations of gas dynamics in a differential conservative form
u t + · F = 0 ,
where u = u ( x , t ) is the vector of conserved variables and F = F ( u ) is the vector of convective fluxes. x is the position vector in a physical domain Ω , and t represents a partial derivative in the temporal dimension. Standard inviscid flux definitions apply.

2.2. Discontinuous Galerkin Discretization

The physical domain Ω and its boundary Ω are discretized with K non-overlapping elements
Ω = i = 1 K D i ,
where D i is a discrete element in the domain. The index i ranges from 1 to the number of elements, K, in the domain. In two dimensions, the domain is divided into triangular or quadrilateral elements. For each discrete element, a numerical solution u h of (1) is constructed as
u u h ( x , t ) = i = 1 N p u ^ i ( t ) ψ i ( x ) = j = 1 N p u ¯ j ( x j , t ) l j ( x ) .
u h can be described either in a modal form with modal basis functions, ψ i , or in a nodal form with Lagrange polynomials, l j . u ^ i ( t ) is the vector of modal coefficients for the corresponding ψ i ( x ) and u ¯ j ( x j , t ) is the vector of nodal coefficients for the corresponding l j ( x ) , where x j indicates the spatial location of u ¯ j . More details about the transformation between the two forms and the selection of the modal basis functions can be found in [15]. In this work, we use Legendre polynomials for ψ i to construct the numerical solution for the modal form. N p is the number of degrees of freedom in each discrete element. Similarly to the modal basis functions, various options exist to distribute the nodal points within an element. For simplicity, equidistant points are used to represent the degrees of freedom. Although distributing nodal points uniformly may lead to ill-conditioned interpolations (Runge’s phenomenon), equidistant nodal points up to polynomial order 4, the upper limit considered in this work, perform robustly.
With this discontinuous Galerkin discretization, the governing Equation (1) can be converted into the following weak form at time t (integration only takes place over the spatial element, D i ):
D i u t l j d V D i F k l j x k d V + D i n · F * l j d A = 0 , j = 1 , , N p ,
where n is the unit outward normal vector at the boundary of D i and F * is a numerical flux at this boundary. In order to stabilize the numerical discretization, the Roe flux [28] is used for the surface integral calculation. The Lax–Friedrichs flux [29] is only used in Section 4.1. For more information about the Roe and the Lax–Friedrichs flux formulations, the reader is referred to [30].

2.3. Space-Time Discretization: ADER-DG

The DG discretization enables simple p-refinement, which is one of its major advantages compared to different higher-order methods. Higher accuracy in space can be achieved by adding more degrees of freedom within an element while the communication remains limited to its face neighbors.
For the integration in time, it is possible to use standard time integration schemes for Equation (4). As the mass matrix for a DG scheme is local for an element, explicit schemes of the Runge–Kutta type can also be used. However, for practical applications, the time step for explicit time integration schemes is very much limited by the CFL condition and only when combined with time-accurate local time-stepping can it be an alternative to implicit methods. Unfortunately, standard time integration schemes do not allow for an order of accuracy beyond first order when applied in combination with time-accurate local time-stepping, which is clearly too low for practical applications.
In order to achieve high-order accuracy in both space and time while using time-accurate local time-stepping, Toro et al. suggested the Arbitrary high-order DERivatives (ADER) method in [31]. In that work, the authors obtain high-order accuracy in the spatial and temporal dimensions by using the Cauchy–Kovalevskaya procedure. Dumbser implemented the ADER method in a finite volume solver in [32]. Similarly, in order to apply the ADER method in combination with DG, we need to introduce a space-time discretization. Here, we summarize the major steps for the DG space-time discretization. Further details can be found in [33,34].

2.3.1. The Local Space-Time DG Predictor

In the predictor step, we look for a solution vector q ( x , t ) which takes the classical DG discretization at time step n as an input and returns its evolution in time up to time step n + 1 , the difference between these time levels being the time step Δ t :
u ( x , t n ) q ( x , t ) t n t t n + 1 = t n + Δ t .
q ( x , t ) can be expanded in space and time as follows:
q ( x , t ) = q ( ξ , τ ) = q ¯ i θ i ( ξ , τ ) ,
where θ i is a polynomial basis function in space-time, q ¯ i is a vector of space-time degrees of freedom, and summation over repeated indices is implied. ξ and τ are parametric coordinates in space and time, respectively. For simplicity, let us restrict ourselves to the one-dimensional case and use parametric coordinates ( τ in time and ξ in space). After multiplication by a space-time test function and integrating over a space-time control volume, Equation (1) becomes
1 1 1 1 θ j ( q τ + f ξ ) d ξ d τ = 0 ,
where f = Δ t 2 ( ξ / x ) · F . Additionally, we assume that the convective fluxes can also be expanded over the space-time basis functions
f = f ¯ i θ i .
Since we are using the nodal representation of a discrete solution, f ¯ i is just the evaluation of the convective fluxes at the space-time degrees of freedom. Therefore, f ¯ i = f ( q ¯ i ) . Note that the assumption (8) is in principle not correct for a nonlinear system of equations and aliasing errors may be introduced. However, this assumption is only used in the predictor step of the algorithm and, as such, the errors do not amplify. With all these equations, we can obtain a final system of algebraic equations for
K i j 1 q ¯ j + K i j ξ · f ¯ j = F i k 1 u ¯ k n ,
where
K i j 1 = 1 1 θ i ( ξ , 1 ) θ j ( ξ , 1 ) d ξ 1 1 1 1 θ i τ θ j d ξ d τ , K i j ξ = 1 1 1 1 θ i θ j ξ d ξ d τ , F i k 1 = 1 1 θ i ( ξ , 1 ) ψ k ( ξ ) d ξ .
The time integration integrals in (10) are approximated by using a quadrature rule with Gauss–Legendre points. After evaluating (10), the unknown q ¯ j can be solved with an iterative procedure up to a desired tolerance. Further details can be found in [33,34]. Note that integration by parts is not used in space to obtain (9). Also note that the space-time discretization described for q in the predictor step is completely local to an element and, therefore, a CFL limit is required to ensure that the procedure is stable.

2.3.2. The Corrector Step

After the space-time predictor solution q ( x , t ) is obtained, we can use it to calculate the solution at the next time level, n + 1 . By integrating the weak form (4) in time from t n to t n + 1 , the equation becomes
( D i l j l k d V ) ( u k n + 1 u k n ) t n t n + 1 D t F k ( u ) l j x k d V d t + t n t n + 1 D t n · F * ( u ) l j d A d t = 0 .
Since we obtained a time evolution of u until t n + 1 from the predictor step, the integrals over the time-volume and time-surface can be calculated with the predictor solution q
( D i l j l k d V ) ( u k n + 1 u k n ) t n t n + 1 D t F k ( q ) l j x k d V d t + t n t n + 1 D t n · F * ( q ) l j d A d t = 0 ,
to obtain the desired solution at the end of the time step.

3. Shock-Capturing in the Discontinuous Galerkin Method

Shock-capturing in the DG method comprises two fundamental steps: (a) detecting elements which may contain a flow discontinuity or may exhibit spurious oscillations, and (b) stabilizing the numerical solution by applying resolving methods in those elements. In [27], we proposed a simple shock detector for the DG method and introduced the functional relationships between the shock detector and filtering strengths. In this section, we revisit the methodology and refine the correlations between the sensor and the filtering method for use by others.

3.1. Shock Detector

Multiple researchers have proposed strategies for detecting discontinuities in the DG method [19,21,22,24]. All shock detectors in the literature can be categorized into three groups: detectors that only use local information, detectors that utilize local information and direct neighbors’ information, and detectors that exploit local information and the Voronoi neighbors’ information. To preserve one of the advantages of the DG method concerning high-performance computing, it is highly desirable to develop shock detectors that only use information from the local element.
Persson and Peraire proposed a novel shock detector that fulfills these demands in [24]. In this shock detector, a spectral decay rate is measured by computing the ratio between the highest modal coefficient and all modal coefficients. By using this ratio, the detector distinguishes between smooth regions and a shock/discontinuity, and it does this particularly well for elements with high polynomial order. However, the detector requires parameter tuning procedures. Inspired by Persson and Peraire’s work, we proposed a shock detector in [27] to circumvent calibrating steps. The basic implementation of the detector is revisited here.
Let ϕ n be the ratio between the norms of the modal coefficients of the highest and the lowest mode in the description of a flow quantity u at time level n:
ϕ n = | u ^ N n | 2 | u ^ 0 n | 2 ,
where u ^ N n is the modal coefficient of the N th degree polynomial at t n . For one-dimensional problems, u ^ N n is just one modal coefficient. In two-dimensional problems, u ^ N n are the modal coefficients for which the degree of the corresponding polynomial is greater than or equal to N. In addition, u ^ 0 n is the modal coefficient of the lowest polynomial at time level n: it is the mean value of u within the element. In this work, we use the density in an element to calculate ϕ n .
The value of ϕ n is usually small unless an element has oscillatory variations in its solution. Indeed, even for an oscillatory solution, ϕ n is not large if the domain is sufficiently refined with discrete elements. However, if a discontinuity is convected from an adjacent element or if it emerges within an element, then ϕ n becomes larger. Therefore, by comparing the values of ϕ at two different time levels, one can narrow down the candidates to flag as troubled elements: the elements that may contain a discontinuity that needs to be treated. Now let β be the ratio between ϕ n + 1 and ϕ n :
β = ϕ n + 1 ϕ n .
We can examine how the strength of an oscillation within an element behaves between two time levels by measuring the β value. With these definitions of ϕ and β , the shock detecting procedure can be stated as follows:
  • Calculate ϕ n from u ( x , t n ) ,
  • March solution forward one time-step and calculate a candidate solution u ( x , t n + 1 ) ,
  • Calculate ϕ n + 1 from u ( x , t n + 1 ) ,
  • If ϕ n + 1 ϕ 0 , an element is not considered to be a troubled element. Otherwise, the current element might be a troubled element. Here, ϕ 0 is a threshold value that depends on the polynomial order of the element.
  • For elements with ϕ n + 1 ϕ 0 , calculate β . If β β 0 , then the current element is a troubled element. Here, β 0 is an empirically determined value.
We set β 0 as 1.0 in order to capture a stationary shock in numerical simulations. In addition, the value of ϕ 0 is predetermined based on the type and polynomial order of an element. Details of determining appropriate ϕ 0 values are described in Section 3.3. Therefore, using the proposed sensor does not require tuning parameters from one problem to the next.

3.2. Filtering

After an element is flagged as a troubled element, additional treatment is required to stabilize or resolve a high-gradient solution near a shock or discontinuity. In this work, we use the filtering method. The filtering method is effectively equivalent to using artificial viscosity. However, it does not restrict the time step compared to directly adding artificial viscosity in the governing equations. In addition, the filtering method amounts to a simple matrix-vector multiplication that can be easily implemented. Therefore, it is not necessary to calculate a viscous flux to resolve a discontinuity, which is necessary when using artificial viscosity. Details of the theoretical background and implementation of the method can be found in [15]. The following second-order exponential filter is used in this work:
σ ( η ) = e x p ( α η 2 ) where η = min ( n / N , 1.0 ) .
Here, N is the polynomial order of a basis function and α is a parameter to change the strength of the filter. The σ values are determined based on the polynomial order and they are multiplicative factors for the modal coefficients, u ^ , of the element to obtain filtered modal coefficients, u ˜
u ˜ = σ u ^ .
For example, let us assume that a filter with α = 1.0 is applied to a one-dimensional n = 2 element. This element has three modal coefficients, u ^ 1 , u ^ 2 , u ^ 3 and their polynomial orders are 0, 1, and 2, respectively. Then, σ ( 0 / 2 ) = e x p ( 1.0 × 0 2 ) = 1.0 is the multiplicative factor for u ^ 1 , σ ( 1 / 2 ) = e x p ( 1.0 × 0 . 5 2 ) 0.7788 is the multiplicative factor for u ^ 2 , and σ ( 2 / 2 ) = e x p ( 1.0 × 1 2 ) 0.3679 is the multiplicative factor for u ^ 3 . Therefore, u ˜ 1 = 1.0 u ^ 1 , u ˜ 2 = 0.7788 u ^ 2 , and u ˜ 3 = 0.3679 u ^ 3 . Note that, for quadrilateral elements in two dimensions, some basis functions have a larger polynomial order value than N. Therefore, η is the minimum value of n / N and 1.0:

3.3. Correlation between α and ϕ n + 1

Researchers seek two objectives in capturing a discontinuity: stability and accuracy. Of these two goals, stability has higher priority because measuring accuracy becomes relevant only if the solution is stable. Therefore, it is necessary to obtain a relationship between α and ϕ n + 1 that provides at least a stable solution when a shock exists within a flow, while adding the least amount of filtering to improve the accuracy of the computed solution.
The shock tube problem introduced in [35] is selected as our benchmark problem to correlate α and ϕ n + 1 . In the shock tube problem, one can change the strength of the right-running shock wave by varying the initial pressure ratio across the diaphragm. We fix the initial conditions just as was done in [35] except for the pressure on the left side of the diaphragm. Starting from an initial pressure ratio of 2, we gradually increased the pressure ratio up to a value of 15 in order to generate a series of shock waves with different strengths. For each initial pressure ratio, we performed numerical experiments with a set of α values, where α is constant in each experiment. During each simulation, we track the maximum ϕ n + 1 values as well as the undershoots that occur near a shock. The L error in density is then calculated with respect to the initial density value on the right side of the diaphragm. Near a discontinuity, it is critical to minimize the L error among other error measurements. This is because the L is directly related to undershoots near the right-running shock in this problem, and these undershoots must be suppressed to guarantee a stable solution throughout the simulations. By performing these numerical experiments, we can obtain a relationship between α and ϕ n + 1 to suppress undershoots as a function of the L percent error.
In Figure 1, the L percent errors in density are depicted as a function of α and ϕ n + 1 for each polynomial order of a triangular element. For a specific ϕ n + 1 value, the minimum α value to ensure that the L error does not exceed a certain value can be found. For example, for n = 1 triangles and ϕ n + 1 = 0.06 , at least α 0.15 is required to guarantee that undershoots do not surpass 2% error. If the L error from the undershoots is required to be lower than 0.5% for the same polynomial order and ϕ n + 1 value, the filtering strength must be larger than 0.25. To generate simple relationships between α and ϕ n + 1 , we constructed a piecewise linear function for each element and polynomial order, provided in [27].
The functional relationships given in [27] ensure that a solution will not diverge due to the presence of a shock. However, if a solution is filtered out excessively, its accuracy will be degraded, thereby losing the main advantage of higher-order methods. Therefore, the functional relationships between α and ϕ n + 1 need to be tailored to simulate a flow as accurately as possible while maintaining stability.
In designing the relationships between α and ϕ n + 1 , we must determine two unknown factors. The first is whether such relationships return excessive filtering strengths. The filtering strengths must be strong enough to stabilize a solution. However, we should not excessively filter a solution to maintain accuracy. The other one is the threshold value selected, ϕ 0 . Theoretically, one can lower the ϕ 0 value all the way to zero, which effectively makes the solver apply filtering everywhere in the domain. On the other hand, if one selects too large a value for ϕ 0 , the sensor cannot robustly distinguish the presence of a shock/discontinuity. Therefore, the ϕ 0 value should be sufficiently large up to a point where it does not harm the robustness of the solver.
To eliminate such ambiguities, we perform other numerical experiments with the shock-entropy interaction problem proposed by Shu and Osher in [36]. The objective of these numerical experiments is to minimize the L 2 error in density with respect to a reference solution. The reference solution is obtained by a one-dimensional simulation with 3,200 equidistant elements of polynomial order two.
In these experiments, we allow two parameters to change: α 0 , an offset value to shift piecewise linear functions, and ϕ 0 , a threshold value for the sensor described in Section 3.1. We optimize the α ( ϕ n + 1 ) relationships by shifting the functions along the α axis as well as by changing their ϕ 0 value. α 0 values range from –0.3 to 0.3, and a typical range for ϕ 0 depends on the element type and polynomial order to find a local minimum of the L 2 error.
Figure 2 shows the normalized L 2 error in density by solving the Shu–Osher problem with triangular elements and different polynomial orders. For a fixed ϕ 0 value, the L 2 error increases if we increase the filtering strength. This is because high-gradient information is lost by excessive filtering. On the other hand, if the filtering strength decreases too much, the error also increases, and eventually the simulation diverges. For a constant α 0 value, the L 2 error increases if we set a value of ϕ 0 that is too large. This is because spurious oscillations are introduced due to a failure of the shock detector. On the contrary, the L 2 error also increases if ϕ 0 becomes too small. In this case, the shock detector determines too many elements as troubled elements and filtering is applied to elements that might not require such resolving processes. In Figure 2, red dots are selected as design points for each polynomial order of a triangular element. After repeating the same procedures for quadrilateral elements, we propose the following relationships between α and ϕ n + 1 , which we use in all other test cases presented in this paper:
(1) Two dimensions, triangular elements:
α t r i , 1 ( ϕ n + 1 ) = ( 3.21 × ϕ n + 1 0.01994 ) 0.00350 ϕ n + 1 , α t r i , 2 ( ϕ n + 1 ) = ( 6.18 × ϕ n + 1 + 0.04383 ) 0.00150 ϕ n + 1 0.01980 , = ( 2.93 × ϕ n + 1 + 0.10818 ) 0.01980 ϕ n + 1 , α t r i , 3 ( ϕ n + 1 ) = ( 35.7 × ϕ n + 1 + 0.03362 ) 0.00020 ϕ n + 1 , α t r i , 4 ( ϕ n + 1 ) = ( 38.0 × ϕ n + 1 + 0.05520 ) 0.00012 ϕ n + 1 .
(2) Two dimensions, quadrilateral elements:
α q u a d , 1 ( ϕ n + 1 ) = ( 4.00 × ϕ n + 1 0.19400 ) 0.06000 ϕ n + 1 , α q u a d , 2 ( ϕ n + 1 ) = ( 3.67 × ϕ n + 1 + 0.01004 ) 0.00300 ϕ n + 1 0.019160 , = ( 2.46 × ϕ n + 1 + 0.03323 ) 0.019160 ϕ n + 1 , α q u a d , 3 ( ϕ n + 1 ) = ( 20.54 × ϕ n + 1 0.014566 ) 0.00200 ϕ n + 1 , α q u a d , 4 ( ϕ n + 1 ) = ( 30.15 × ϕ n + 1 + 0.013116 ) 0.00090 ϕ n + 1 .

4. Numerical Simulation Results

4.1. Isentropic Vortex Problem

We first solve the isentropic vortex problem introduced in [7] for measuring the performance of our shock detector in the absence of discontinuities. This problem is well suited for evaluating whether a numerical scheme can preserve the initial shape of a vortex as well any initial disturbances in flow quantities. To obtain the optimum accuracy from the DG discretization, the shock detector should not flag any element in the domain as a troubled element, as the solution to this problem is smooth. Otherwise, the accuracy of a simulation will degrade due to the misuse of the filtering process. The initial non-dimensional conditions for this problem are as follows:
( ρ , u , v , p ) = ( [ 1 ( γ 1 ) ϵ 2 8 γ π 2 e ( 1 r 2 ) ] 1 γ 1 , 1 ϵ 2 π e 0.5 ( 1 r 2 ) y ¯ , 1 + ϵ 2 π e 0.5 ( 1 r 2 ) x ¯ , ρ γ ) ,
where ( x ¯ , y ¯ ) = ( x 5 , y 5 ) , r 2 = x ¯ 2 + y ¯ 2 , and ϵ = 5.0 . The physical domain becomes Ω = [ 0 , 10 ] × [ 0 , 10 ] and it is discretized with two different types of meshes: a fully unstructured mesh with triangular elements and a structured mesh with quadrilateral elements.The CFL value is 0.5 for calculating a global time step. As the time step for stability does not vary much over the domain, it is not necessary to use the time-accurate local time-stepping capabilities of ADER-DG.
The L 2 errors in density are listed in Table 1 and Table 2. Except for grid levels 1 to 3 in the n = 2 quadrilateral element grid, the L 2 errors are identical regardless of the filtering. Even for the quadrilateral elements with polynomial order two, the optimum accuracy is recovered as grids are refined. From this numerical test, the sensor and the filtering method proposed in Section 3 can identify a smooth solution in a flow and does not degrade the accuracy of the solution for both triangular and quadrilateral elements.

4.2. Shock Tube Problem

This numerical test case is a fundamental problem to prove whether a numerical method can capture discontinuities [35]. After the diaphragm is removed at t = 0 , a left-running rarefaction wave, a right-running contact discontinuity, and a right-running shock wave are created. Since C 0 -continuous and discontinuous physical waves both exist in the problem, this test case is well suited to examine the performance of our shock-capturing capability. The initial conditions at t = 0 are given by
( ρ , u , v , p ) = ( 1.0 , 0.0 , 0.0 , p L ) , i f x 0.5 , = ( 0.125 , 0.0 , 0.0 , p R ) , otherwise ,
where p L and p R are pressure values on the left and right sides of the diaphragm, respectively. The CFL value used for calculating the explicit time step is 0.45 and each simulation runs until the right-running shock wave arrives at x = 0.95 . For comparison, an exact solution is calculated using the method in [30] based on the selected initial pressure ratio at t = 0 . In Sod’s original work in [35], p L = 1.0 and p R = 0.1 . However, in our numerical experiments, we only fix the value of p R to 0.1 and set p L values to 0.2, 0.55, and 1.0. By selecting these p L values, we can generate normal shock waves with normal Mach numbers of approximately 1.1, 1.38, and 1.66, respectively.
We choose the numerical domain as Ω = [ 0 , 1 ] × [ 0 , 0.1 ] . The domain is discretized with triangles or quadrilaterals of characteristic length 0.01. Figure 3 shows the three different meshes used in the simulations: a structured mesh with quadrilateral elements, a structured mesh decomposed into triangular elements, and an unstructured mesh with triangular elements. Periodic boundary conditions are imposed on the surface elements located at y = 0.0 and y = 0.1 .
Numerical simulation results are shown in Figure 4, Figure 5 and Figure 6. In each figure, the dashed black line represents the exact solution and the solid blue line is the numerical simulation result with the proposed sensor and filtering approach. Numerical solutions are interpolated along the line at y = 0.054 with 501 equidistant sampling points. The shock detector identifies the shock and the contact discontinuity both in the two-dimensional structured and unstructured meshes. Furthermore, the filtering is appropriately applied to elements flagged as troubled elements to resolve both weak and strong discontinuities.

4.3. Shu–Osher Problem

This problem was suggested by Shu and Osher in [36]. The physical domain is divided into a high-pressure part on the left side and a sinusoidal density wave on the right side of a diaphragm. Then, a Mach 3.0 shock wave propagates to the right and interacts with the sinusoidal wave after the diaphragm is removed. Robust and accurate shock-capturing is necessary to preserve the initial high-frequency behavior as well as local extrema after the shock is propagated. The initial non-dimensional conditions at t = 0 are:
( ρ , u , v , p ) = ( 27 7 , 2.629369 , 0.0 , 31 3 ) i f x 0.1 , = ( 1 + 0.2 sin ( 50 x ) , 0.0 , 0.0 , 1.0 ) otherwise .
The numerical domain for this problem is Ω = [ 0 , 1 ] × [ 0 , 0.1 ] and the domain is discretized similarly as in Figure 3: a structured mesh with quadrilateral elements, a structured mesh decomposed into triangular elements, and an unstructured mesh with triangular elements. However, and unlike in the shock tube cases, we intentionally fix the total number of degrees of freedom along the x-direction to 1200. Consequently, elements with n = 2, 3, and 4 have characteristic lengths of 0.0025 , 0.0033 , and 0.004 , respectively, regardless of the element type. By restricting the total number of degrees of freedom, we can compare the efficiency of the shock-capturing qualitatively as a function of the polynomial order of the discretized elements. Each simulation runs until t = 0.18 and the CFL value is 0.3 . A reference solution is obtained with 3,200 second-order elements in one dimension with the proposed sensor and filtering method in [27]. Periodic boundary conditions are imposed on the surface elements located at y = 0.0 and y = 0.1 .
Simulation results with different element types and polynomial orders are shown in Figure 7. The results are qualitatively the same for quadrilaterals regardless of the polynomial order. For triangular grids, polynomial order two elements preserve local extrema and high-frequency behavior almost the same as quadrilaterals, but the accuracy of the solutions for polynomial order three and four degrades significantly. However, in [27], we found that triangular elements with polynomial order three and four can also preserve accuracy after the shock if the characteristic length of the triangles is 0.0025 . Therefore, we conclude that polynomial order two is cost effective in capturing shocks and pressure oscillations in these situations, compared to polynomial orders three and four for triangular elements with the sensor and the filtering methods.

4.4. Shock–Vortex Interaction

This problem was introduced by Jiang and Shu in [37]. A stationary shock resides within the domain and a right-running vortex interacts with a shock as it goes through it. To capture the highly-dynamic flow after the vortex traverses the shock, the detector and the filter must capture shocks robustly and precisely.
The spatial domain used is Ω = [ 0 , 2 ] × [ 0 , 1 ] , and it is discretized using a structured grid with quadrilateral elements. The grid spacing is uniform in the y-direction, but a finer spacing is used in the x-direction around the stationary shock location, x = 0.5 . Similarly to the Shu–Osher problem, we fix the total number of degrees of freedom to 600 along the x-direction to evaluate the efficiency of capturing shocks as a function of the polynomial order of the quadrilateral elements. The non-dimensional initial conditions of the left side of the shock at t = 0 are as follows:
( ρ , u , v , p ) = ( ( 1 ( γ 1 ) ϵ 2 e 2 α ( 1 τ 2 ) 4 α γ ) 1 γ 1 , M γ + ϵ τ e α ( 1 τ 2 ) s i n θ , ϵ τ e α ( 1 τ 2 ) c o s θ , ρ γ ) ,
where M = 1.1 , tan θ = ( y y c ) / ( x x c ) , τ = r / r c , r = ( x x c ) 2 + ( y y c ) 2 , ϵ = 0.3 , r c = 0.05 , ( x c , y c ) = ( 0.25 , 0.5 ) , and α = 0.204 . Initial conditions on the right side of the shock are calculated from the Rankine–Hugoniot conditions. Inviscid wall boundary conditions are imposed on the surface elements located at y = 0.0 and y = 1.0 . The CFL value used is 0.1 and each simulation runs until t = 0.8 .
Figure 8 shows isolines of pressure at t = 0.05, 0.2, and 0.35. The pressure contours are qualitatively identical for all polynomial orders at any time. Note that, in these figures, the isolines are not shown in some areas around the shock. This is especially the case for regions where the vortex has not disturbed the flow yet. As the value of the solution at the interfaces between elements is multiply defined in the DG discretization, i.e., allowing a jump in the solution, and the shock is located exactly at these interfaces, the postprocessing software does not show the isolines in these regions.
Pressure contours at t = 0.6 are shown in Figure 9. For comparison purposes, we imported the result by [37] in Figure 9a. Pressure contours are qualitatively the same. Among three different polynomial orders, our n = 2 solution exhibits less numerical noise compared to the n = 3 and n = 4 solutions.

5. Conclusions

In this work, we further examined whether the shock detector in [27] can properly distinguish between shocks/discontinuities and regions of the flow with smooth solutions. In addition, we proposed problem-independent, optimized relationships between the shock detector and the filtering strength as a function of polynomial order and element type in two dimensions. By using the suggested relationships, weak and strong discontinuities can be captured robustly with no additional parameter tuning for two-dimensional elements with polynomial orders 2, 3, and 4. Moreover, the proposed method is based only on operations which utilize local information. The shock sensor and the filtering method have been implemented in the discontinuous Galerkin finite element method (DG-FEM) solver in the SU2 framework [38] for users to easily utilize and reproduce these results.
The whole method is extensible to three-dimensional elements as well as viscous problems in a straightforward way. In the future, we expect our implementation of the SU2 DG-FEM solver to be applied to problems where interactions between shock waves and turbulence are present. Additionally, we would like to exploit the proposed method in realistic industrial cases that require highly-scalable methods.

Author Contributions

Conceptualization, formal analysis, investigation, methodology, software, validation, visualization, and writing—original draft: J.H.C.; Supervision, project administration, writing—review and editing: J.J.A. and E.v.d.W.

Funding

This research was funded by the United States Department of Energy under the Predictive Science Academic Alliance Program 2 (PSAAP2) at Stanford University.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Wang, Q.; Ren, Y.X.; Li, W. Compact high order finite volume method on unstructured grids I: Basic formulations and one-dimensional schemes. J. Comput. Phys. 2016, 314, 863–882. [Google Scholar] [CrossRef]
  2. Gao, X.; Owen, L.D.; Guzik, S.M. A high-order finite-volume method for combustion. In Proceedings of the AIAA Scitech 2016 Forum, San Diego, CA, USA, 4–8 January 2016; pp. 1–14. [Google Scholar]
  3. Lele, S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 1992, 103, 16–42. [Google Scholar] [CrossRef]
  4. Deng, X.; Mao, M.; Tu, G.; Zhang, H.; Zhang, Y. High-order and high accurate CFD methods and their applications for complex grid problems. Commun. Comput. Phys. 2012, 11, 1081–1102. [Google Scholar] [CrossRef]
  5. Allahyari, M.; Mohseni, K. Numerical simulation of flows with shocks and turbulence using observable methodology. In Proceedings of the AIAA Scitech 2018 Forum, Kissimmee, FL, USA, 8–12 January 2018; pp. 1–10. [Google Scholar]
  6. Jause-Labert, C.; Godeferd, F.; Favier, B. Numerical validation of the volume penalization method in three-dimensional pseudo-spectral simulations. Comput. Fluids 2012, 67, 41–56. [Google Scholar] [CrossRef]
  7. Shu, C.W. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. In Advanced Numerical Approximation of Nonlinear Hyperbolic Equations; Springer: Berlin/Heidelberg, Germany, 1998; pp. 325–442. [Google Scholar]
  8. van der Weide, E.; Giangaspero, G.; Svärd, M. Efficiency benchmarking of an energy stable high-order finite difference discretization. AIAA J. 2015, 53, 1845–1860. [Google Scholar] [CrossRef]
  9. Reed, W.; Hill, T. Triangular Mesh Methods for the Neutron Transport Equation; Technical Report LA-UR-73-479; Los Alamos Scientific Laboratory: Los Alamos, NM, USA, 1973. [Google Scholar]
  10. Cockburn, B.; Shu, C.W. The Runge–Kutta local projection P1-discontinuous Galerkin finite element method for scalar conservation laws. Math. Model. Numer. Anal. 1991, 25, 337–361. [Google Scholar] [CrossRef]
  11. Cockburn, B.; Shu, C.W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework. Math. Comput. 1989, 52, 411–435. [Google Scholar]
  12. Cockburn, B.; Lin, S.Y.; Shu, C.W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems. J. Comput. Phys. 1989, 84, 90–113. [Google Scholar] [CrossRef] [Green Version]
  13. Cockburn, B.; Hou, S.; Shu, C.W. The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV: The multidimensional case. Math. Comput. 1990, 54, 545–581. [Google Scholar]
  14. Cockburn, B.; Shu, C.W. The Runge–Kutta discontinuous Galerkin method for conservation laws V: Multidimensional systems. J. Comput. Phys. 1998, 141, 199–224. [Google Scholar] [CrossRef]
  15. Hesthaven, J.S.; Warburton, T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, 1st ed.; Springer: New York, NY, USA, 2008. [Google Scholar]
  16. Atkins, H.L.; Pampell, A. Robust and accurate shock capturing method for high-order discontinuous Galerkin methods. In Proceedings of the 20th AIAA Computational Fluid Dynamics Conference, Honolulu, HI, USA, 27–30 June 2011; pp. 1–22. [Google Scholar]
  17. Barter, G.; Darmofal, D. Shock capturing with higher-order, PDE-based artificial viscosity. In Proceedings of the 18th AIAA Computational Fluid Dynamics Conference, Miami, FL, USA, 25–28 June 2007; pp. 1–14. [Google Scholar]
  18. Burbeau, A.; Sagaut, P.; Bruneau, C.H. A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods. J. Comput. Phys. 2001, 169, 111–150. [Google Scholar] [CrossRef]
  19. Clain, S.; Diot, S.; Loubère, R. A high-order finite volume method for systems of conservation laws—Multi -dimensional optimal order detection (MOOD). J. Comput. Phys. 2011, 230, 4028–4050. [Google Scholar] [CrossRef]
  20. Dumbser, M.; Loubère, R. A simple robust and accurate a posteriori sub-cell finite volume limiter for the discontinuous Galerkin method on unstructured meshes. J. Comput. Phys. 2016, 319, 163–199. [Google Scholar] [CrossRef] [Green Version]
  21. Klöckner, A.; Warburton, T.; Hesthaven, J.S. Viscous shock capturing in a time-explicit discontinuous Galerkin method. Math. Model. Nat. Phenom. 2011, 6, 57–83. [Google Scholar] [CrossRef]
  22. Lv, Y.; See, Y.C.; Ihme, M. An entropy-residual shock detector for solving conservation laws using high-order discontinuous Galerkin methods. J. Comput. Phys. 2016, 322, 448–472. [Google Scholar] [CrossRef]
  23. Nguyen, C.; Peraire, J. An adaptive shock-capturing HDG method for compressible flows. In Proceedings of the 20th AIAA Computational Fluid Dynamics Conference, Honolulu, HI, USA, 27–30 June 2011; pp. 1–20. [Google Scholar]
  24. Persson, P.O.; Peraire, J. Sub-cell shock capturing for discontinuous Galerkin methods. In Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 9–12 January 2006; pp. 1–13. [Google Scholar]
  25. Sheshadri, A. An analysis of stability of the flux reconstruction formulation with applications to shock capturing. Ph.D. Thesis, Stanford University, Stanford, CA, USA, 2016. [Google Scholar]
  26. Vuik, M.J.; Ryan, J.K. Multiwavelet troubled-cell indicator for discontinuity detection of discontinuous Galerkin schemes. J. Comput. Phys. 2014, 270, 138–160. [Google Scholar] [CrossRef] [Green Version]
  27. Choi, J.H.; Alonso, J.J.; van der Weide, E. Simple shock detector for discontinuous Galerkin method. In Proceedings of the AIAA Scitech 2019 Forum, San Diego, CA, USA, 7–11 January 2019. [Google Scholar]
  28. Roe, P.L. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 1981, 43, 357–372. [Google Scholar] [CrossRef]
  29. Lax, P.D. Weak solutions of nonlinear hyperbolic equations and their numerical computation. Commun. Pure Appl. Math. 1954, 7, 159–193. [Google Scholar] [CrossRef]
  30. Toro, E.F. Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction; Springer: Berlin/Heidelberg, Germany, 2009; Chapter 5,6,7,11. [Google Scholar]
  31. Toro, E.F.; Millington, R.C.; Nejad, L. Towards very high order Godunov schemes. In Godunov Methods: Theory and Applications; Springer: Boston, MA, USA, 2001; pp. 907–940. [Google Scholar]
  32. Dumbser, M.; Enaux, C.; Toro, E.F. Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws. J. Comput. Phys. 2008, 227, 3971–4001. [Google Scholar] [CrossRef]
  33. Rezzolla, L.; Zanotti, O. Relativistic Hydrodynamics; Oxford University Press: Oxford, UK, 2013. [Google Scholar]
  34. Zanotti, O. Discontinuous Galerkin Methods for Hyperbolic PDEs. Lecture 2. In Proceedings of the Prospects in Theoretical Physics 2016, Princeton, NJ, USA, 18–29 July 2016; pp. 1–59. Available online: https://static.ias.edu/pitp/2016/node/1026.html (accessed on 10 July 2019).
  35. Sod, G. Survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J. Comput. Phys. 1978, 27, 1–31. [Google Scholar] [CrossRef]
  36. Shu, C.W.; Osher, S. Efficient implementation of essentially non-oscillatory shock-capturing schemes II. J. Comput. Phys. 1989, 83, 32–78. [Google Scholar] [CrossRef]
  37. Jiang, G.S.; Shu, C.W. Efficient Implementation of Weighted ENO Schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef] [Green Version]
  38. SU2, The Open-Source CFD Code. Available online: https://su2code.github.io (accessed on 1 June 2019).
Figure 1. L percent errors in density from a series of shock tube problems with various initial pressure ratios. The errors are plotted as a function of α and ϕ n + 1 for each polynomial order of triangular elements.
Figure 1. L percent errors in density from a series of shock tube problems with various initial pressure ratios. The errors are plotted as a function of α and ϕ n + 1 for each polynomial order of triangular elements.
Energies 12 02651 g001
Figure 2. Normalized L 2 errors in density for the Shu–Osher problem for each polynomial order of triangular elements. Red dots are selected design points for α 0 and ϕ 0 . For diverged simulation cases, the normalized L 2 error values are replaced with the maximum L 2 error value among stable solutions.
Figure 2. Normalized L 2 errors in density for the Shu–Osher problem for each polynomial order of triangular elements. Red dots are selected design points for α 0 and ϕ 0 . For diverged simulation cases, the normalized L 2 error values are replaced with the maximum L 2 error value among stable solutions.
Energies 12 02651 g002
Figure 3. Three different meshes for shock tube problem: (a) structured mesh with quadrilateral elements; (b) structured mesh with triangular elements; (c) unstructured mesh with triangular elements.
Figure 3. Three different meshes for shock tube problem: (a) structured mesh with quadrilateral elements; (b) structured mesh with triangular elements; (c) unstructured mesh with triangular elements.
Energies 12 02651 g003
Figure 4. Two-dimensional shock tube problem with initial pressure ratio 2, various polynomial orders, and different types of grids. In each plot, the top-right inset figure is the result near the shock wave and the bottom-left inset figure is the result near the contact discontinuity.
Figure 4. Two-dimensional shock tube problem with initial pressure ratio 2, various polynomial orders, and different types of grids. In each plot, the top-right inset figure is the result near the shock wave and the bottom-left inset figure is the result near the contact discontinuity.
Energies 12 02651 g004
Figure 5. Two-dimensional shock tube problem with initial pressure ratio 5.5, various polynomial orders, and different types of grids. In each plot, the top-right inset figure is the result near the shock wave and the bottom-left inset figure is the result near the contact discontinuity.
Figure 5. Two-dimensional shock tube problem with initial pressure ratio 5.5, various polynomial orders, and different types of grids. In each plot, the top-right inset figure is the result near the shock wave and the bottom-left inset figure is the result near the contact discontinuity.
Energies 12 02651 g005aEnergies 12 02651 g005b
Figure 6. Two-dimensional shock tube problem with initial pressure ratio 10, various polynomial orders, and different types of grids. In each plot, the top-right inset figure is the result near the shock wave and the bottom-left inset figure is the result near the contact discontinuity.
Figure 6. Two-dimensional shock tube problem with initial pressure ratio 10, various polynomial orders, and different types of grids. In each plot, the top-right inset figure is the result near the shock wave and the bottom-left inset figure is the result near the contact discontinuity.
Energies 12 02651 g006aEnergies 12 02651 g006b
Figure 7. Two-dimensional Shu–Osher problem with β 0 = 1.00 . The dashed black line is the reference solution and the solid red, blue, black lines are interpolated solutions with 1,201 equidistant points for n = 2 , 3 , and 4 elements, respectively.
Figure 7. Two-dimensional Shu–Osher problem with β 0 = 1.00 . The dashed black line is the reference solution and the solid red, blue, black lines are interpolated solutions with 1,201 equidistant points for n = 2 , 3 , and 4 elements, respectively.
Energies 12 02651 g007
Figure 8. Two-dimensional shock-vortex interaction problem at various times. Forty-one pressure contours are used from 0.84 to 1.50.
Figure 8. Two-dimensional shock-vortex interaction problem at various times. Forty-one pressure contours are used from 0.84 to 1.50.
Energies 12 02651 g008
Figure 9. Two-dimensional shock-vortex interaction problem at t = 0.6 . 90 pressure contours are used from 1.19 to 1.37. (a) WENO-LF-5 results by Jiang et al.
Figure 9. Two-dimensional shock-vortex interaction problem at t = 0.6 . 90 pressure contours are used from 1.19 to 1.37. (a) WENO-LF-5 results by Jiang et al.
Energies 12 02651 g009
Table 1. L 2 error of density on unstructured grids with triangular elements.
Table 1. L 2 error of density on unstructured grids with triangular elements.
Grid LevelElement Lengthn = 2n = 3n = 4
FilteringW/O FilteringFilteringW/O FilteringFilteringW/O Filtering
110/42.112E-022.112E-022.708E-022.708E-021.332E-021.332E-02
210/61.598E-021.598E-021.422E-021.422E-024.114E-034.114E-03
310/89.339E-039.339E-036.474E-036.474E-031.135E-031.135E-03
410/123.376E-033.376E-031.179E-031.179E-031.434E-041.434E-04
510/161.843E-031.843E-033.415E-043.415E-043.417E-053.417E-05
Table 2. L 2 error of density on structured grids with quadrilateral elements.
Table 2. L 2 error of density on structured grids with quadrilateral elements.
Grid LevelElement Lengthn = 2n = 3n = 4
FilteringW/O FilteringFilteringW/O FilteringFilteringW/O Filtering
110/43.535E-023.931E-022.580E-022.580E-021.815E-021.815E-02
210/61.643E-021.648E-021.130E-021.130E-026.296E-036.296E-03
310/81.300E-021.251E-026.668E-036.668E-031.033E-031.033E-03
410/125.167E-035.167E-032.611E-032.611E-031.440E-041.440E-04
510/162.099E-032.099E-038.461E-048.461E-043.923E-053.923E-05

Share and Cite

MDPI and ACS Style

Choi, J.H.; Alonso, J.J.; van der Weide, E. A Simple and Robust Shock-Capturing Approach for Discontinuous Galerkin Discretizations. Energies 2019, 12, 2651. https://doi.org/10.3390/en12142651

AMA Style

Choi JH, Alonso JJ, van der Weide E. A Simple and Robust Shock-Capturing Approach for Discontinuous Galerkin Discretizations. Energies. 2019; 12(14):2651. https://doi.org/10.3390/en12142651

Chicago/Turabian Style

Choi, Jae Hwan, Juan J. Alonso, and Edwin van der Weide. 2019. "A Simple and Robust Shock-Capturing Approach for Discontinuous Galerkin Discretizations" Energies 12, no. 14: 2651. https://doi.org/10.3390/en12142651

APA Style

Choi, J. H., Alonso, J. J., & van der Weide, E. (2019). A Simple and Robust Shock-Capturing Approach for Discontinuous Galerkin Discretizations. Energies, 12(14), 2651. https://doi.org/10.3390/en12142651

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop