Search

In this paper, we present a high-order arbitrary Lagrangian Eulerian (ALE) method for compressible single-material flow on adaptive moving unstructured meshes. The discretization of system is implemented by Runge-Kutta Discontinuous Galerkin (RKDG) method. The vertex velocity is given by the approach of mesh movement from G. Chen et al. (2008) [8] . The new mesh can be automatically redistributed and concentrated on the regions with large gradient value of the variables. A HWENO reconstruction is used to eliminate false oscillations which maintains compactness of DG methods. In addition, we prove the property that the material derivatives of the Lagrangian basis functions are equal to zero with which the scheme is shown to satisfy the geometric conservation law (GCL). Furthermore, the scheme is conservative for mass, momentum and total energy. Numerical examples are presented to illustrate the good performance and high-order accuracy of the scheme.

In the present paper, a multi-step reconstruction procedure is proposed for high order finite volume schemes on unstructured grids using compact stencil. The procedure is a recursive algorithm that can eventually provide sufficient relations for high order reconstruction in a multi-step procedure. Two key elements of this procedure are the partial inversion technique and the continuation technique. The partial inversion can be used not only to obtain lower order reconstruction based on existing reconstruction relations, but also to regularize the existing reconstruction relations to provide new relations for higher order reconstructions. The continuation technique is to extend the regularized relations on the face-neighboring cells to current cell as additional reconstruction relations. This multi-step procedure is operationally compact since in each step only the relations defined on a compact stencil are used. In the present paper, the third and fourth order finite volume schemes based on two-step quadratic and three-step cubic reconstructions are studied. (C) 2019 Elsevier Inc. All rights reserved.

In this paper, we propose a multi-symplectic quasi-interpolation method for solving multi-symplectic Hamiltonian partial differential equations. Based on the method of lines, we first discretize the multi-symplectic PDEs using quasi-interpolation method and then employ appropriate time integrators to obtain the full-discrete system. The local conservation properties including multi-symplectic conservation laws, energy conservation laws and momentum conservation laws are discussed in detail. For illustration, we provide two concrete examples: the nonlinear wave equation and the nonlinear Schrodinger equation. The salient feature of our multi-symplectic quasi-interpolation method is that it is valid both on uniform grids and nonuniform grids. The numerical results show the good accuracy and excellent conservation properties of the proposed method.

In this paper, we continue our work in [46] and propose a new type of high-order finite volume multi-resolution weighted essentially non-oscillatory (WENO) schemes to solve hyperbolic conservation laws on triangular meshes. Although termed “multi-resolution WENO schemes”, we only use the information defined on a hierarchy of nested central spatial stencils and do not introduce any equivalent multi-resolution representation. We construct new third-order, fourth-order, and fifth-order WENO schemes using three or four unequal-sized central spatial stencils, different from the classical WENO procedure using equal-sized biased/central spatial stencils for the spatial reconstruction. The new WENO schemes could obtain the optimal order of accuracy in smooth regions, and could degrade gradually to first-order of accuracy so as to suppress spurious oscillations near strong discontinuities. This is the first time that only a series of unequal-sized hierarchical central spatial stencils are used in designing arbitrary high-order finite volume WENO schemes on triangular meshes. The main advantages of these schemes are their compactness, robustness, and their ability to maintain good convergence property for steady-state computation. The linear weights of such WENO schemes can be any positive numbers on the condition that they sum to one. Extensive numerical results are provided to illustrate the good performance of these new finite volume WENO schemes.

Liu and Jiao proposed a WLS-ENO method (Weighted-Least-Squares based Essentially Non-Oscillatory). Comparing to the WENO method (Weighted ENO), the WLS-ENO method is more flexible for solving the hyperbolic conservation laws. Notice that, a WLS problem is generated for each WLS-ENO procedure, although the QR factorization with column pivoting is employed, the computational cost of solving WLS problem is still heavy. So in this paper, we will propose an improved WLS-ENO method for improving the efficiency. The basic idea of new method is that, let the reconstructed polynomial approximate the cell averages of some smooth cells in conservation fashion, and approximate other cell averages in WLS fashion, so that the resulted WLS problem only has one unknown. Because the resulted WLS problem can be solved by using the normal equation directly, the whole computational cost is small, thus the new method has better efficiency.

A numerical model based on decomposed vaporization mechanisms has been developed for the study of vaporization of liquid droplets with large deformation and high mass transfer rate. In this model, the vaporization-induced volume, energy, and species source terms in the gas phase are modeled using a body-fitted layer placed around the interface on the gas side. In this layer, continuous distributions of volume, energy, and species sources are specified by the conservation laws. A corresponding body-fitted sink layer is placed on the liquid side of the interface to account for mass, energy and species loss in the liquid phase due to vaporization. The formulation retains the advantages of interface localization methods and mitigates the error and complexity caused by the combined treatment of mass/heat transfer and interface localization. The new model is versatile and does not depend on the specific interface localization method, and can thus be implemented in any of the existing methods. As a specific example, the present approach is implemented in a volume-of-fraction (VOF)-based code “Gerris”. A number of test cases concerning liquid vaporization in both quiescent and convective environments are presented, with the constant-density assumption, to verify and validate this model in terms of interface localization. In particular, the vaporization of a single n-decane droplet in a uniform air flow of 1000 K at 1 atm is investigated. The predicted liquid volume evolution demonstrates the robustness and accuracy of interface localization for cases with large density ratio and high vaporization rate. Vaporization of n-decane droplets with large deformation and high mass transfer is also investigated systematically.

In this paper a new three-dimensional Eulerian interface-capturing model is proposed for the simulation of compressible solid-fluid problems. The model is a development of a well-established five equation multi-fluid model to include material strength, by incorporating a new kinematic evolution equation for the elastic stretch tensor, and augmenting the Mie-Gruneisen equation-of-state to include a contribution from elastic strain. Principal to the development of the model is the fact that a great number of solid materials can be described by the particular general equation-of-state framework, which reduces to the Mie-Gruneisen equation-of-state in the limit of zero strength and is thus equally applicable to a wide class of fluids. The constitutive models are founded on hyperelastic theory and are therefore thermodynamically compatible. Only one kinematic equation is required for arbitrary numbers of components by assuming that mixtures are described by a common deviatoric strain tensor. The system of evolution equations can be written in conservation law form, and are therefore suitable for application of Godunov's method; specifically an HLLD Riemann solver is formulated for the calculation of numerical fluxes. Interfaces are sharpened by using the THINC method in conjunction with MUSCL reconstruction. A time operator splitting strategy is used to divide the time integration into an explicit elastic update, which uses the third order TVD Runge-Kutta method, followed by an implicit plastic update. The method is verified using a number of challenging solid-fluid problems, including a high-velocity impact of a viscoplastic strain-hardening cylinder in air. The potential of the model for practical problems is demonstrated through three-dimensional simulation of an explosive buried in solid ground material.

We present an efficient nodal discontinuous Galerkin method for approximating nearly incompressible flows using the Boltzmann equations. The equations are discretized with Hermite polynomials in velocity space yielding a first order conservation law. A stabilized unsplit perfectly matching layer (PML) formulation is introduced for the resulting nonlinear flow equations. The proposed PML equations exponentially absorb the difference between the nonlinear fluctuation and the prescribed mean flow. We introduce semi-analytic time discretization methods to improve the time step restrictions in small relaxation times. We also introduce a multirate semi-analytic Adams–Bashforth method which preserves efficiency in stiff regimes. Accuracy and performance of the method are tested using distinct cases including isothermal vortex, flow around square cylinder, and wall mounted square cylinder test cases.

We propose a general class of genuinely two-dimensional incomplete Riemann solvers for systems of conservation laws. In particular, extensions of Balsara's multidimensional HLL scheme [J. Comput. Phys. 231 (2012) 7476-7503] to two-dimensional PVM/RVM (Polynomial/Rational Viscosity Matrix) finite volume methods are considered. The numerical flux is constructed by assembling, at each edge of the computational mesh, a one-dimensional PVM/RVM flux with two purely two-dimensional PVM/RVM fluxes at vertices, which take into account transversal features of the flow. The proposed methods are applicable to general hyperbolic systems, although in this paper we focus on applications to magnetohydrodynamics. In particular, we propose an efficient technique for divergence cleaning of the magnetic field that provides good results when combined with our two-dimensional solvers. Several numerical tests including genuinely two-dimensional effects are presented to test the performances of the proposed schemes.

International audience; In this work, we consider the discretization of nonlinear hyperbolic systems in nonconservative form with the high-order discontinuous Galerkin spectral element method (DGSEM) based on collocation of quadrature and interpolation points (Kopriva and Gassner, J. Sci. Comput., 44 (2010), pp.136-155). We present a general framework for the design of such schemes that satisfy a semi-discrete entropy inequality for a given convex entropy function at any approximation order. The framework is closely related to the one introduced for conservation laws by Fisher and Carpenter (J. Comput. Phys., 252 (2013), pp. 518-557) and Gassner (SIAM J. Sci. Comput., 35 (2013), pp. A1233-A1253) and relies on the modification of the integral over discretization elements where we replace the physical fluxes by entropy conservative numerical fluxes from Castro et al. (SIAM J. Numer. Anal., 51 (2013), pp. 1371-1391), while entropy stable numerical fluxes are used at element interfaces. Time discretization is performed with strong-stability preserving Runge-Kutta schemes. The present method is first validated on the Burgers and Euler equations in nonconservative form. We then use this framework for the discretization of two-phase flow systems in one space-dimension: a 2 × 2 system with a nonconservative product associated to a linearly-degenerate field for which the DGSEM fails to capture the physically relevant solution, and the isentropic Baer-Nunziato model. For the latter, we derive conditions on the numerical parameters of the discrete scheme to further keep positivity of the partial densities and a maximum principle on the void fractions. Numerical experiments support the conclusions of the present analysis and highlight stability and robustness of the present schemes.

In the article [Yang et al. (2012) [37] ], the authors have developed a high order moving mesh WENO method for one-dimensional (1D) hyperbolic conservation laws, which is shown to be effective in resolving shocks and other complex solution structures. In this paper, in the light of the similar moving mesh technique, we develop a novel WENO scheme with non-polynomial bases, in particular, the exponential bases to further improve the performance of WENO schemes for solving 1D conservation laws. Furthermore, we modify the original moving mesh technique by developing a new monitor function as well as a different mesh smoothing strategy. A collection of numerical examples is presented to demonstrate high order accuracy and robustness of the method in capturing smooth and non-smooth solutions including the strong δ shock arising from the weakly hyperbolic pressureless Euler equations.

We construct entropy conservative and entropy stable discontinuous Galerkin (DG) discretizations for time-dependent nonlinear hyperbolic conservation laws on curvilinear meshes. The resulting schemes preserve a semi-discrete quadrature approximation of a continuous global entropy inequality. The proof requires the satisfaction of a discrete geometric conservation law, which we enforce through an appropriate polynomial approximation. We extend the construction of entropy conservative and entropy stable DG schemes to the case when high order curvilinear mass matrices are approximated using low-storage weight-adjusted approximations, and describe how to retain global conservation properties under such an approximation. For certain types of curvilinear meshes, these weight-adjusted approximations deliver optimal rates of convergence. Finally, the high order accuracy, local conservation, and discrete conservation or dissipation of entropy for under-resolved solutions are verified through numerical experiments for the compressible Euler equations on triangular and tetrahedral meshes.

A high-order numerical method for interface-preserving level-set reinitialization is presented in this paper. In the interface cells, the gradient of the level-set function is determined by a weighted local projection scheme and the missing additive constant is determined such that the position of the zero level set is preserved. In the non-interface cells, we compute the gradient of the level-set function by solving a Hamilton–Jacobi equation as a conservation law system using the discontinuous Galerkin method, following the work by Hu and Shu [SIAM J. Sci. Comput. 21 (1999) 660–690]; the missing constant is then recovered by the continuity of the level-set function while taking into account the characteristics. To handle highly distorted initial conditions, we develop a hybrid numerical flux that combines the Lax–Friedrichs flux and the penalty flux. Our method is stable for non-trivial test cases and handles singularities away from the interface very well. When derivative singularities are present on the interface, a second-derivative limiter is designed to suppress the oscillations. At least ( N + 1 ) th order accuracy in the interface cells and Nth order in the whole domain are observed for smooth solutions when Nth degree polynomials are used. Two dimensional test cases are presented to demonstrate superior properties such as accuracy, long-term stability, interface-preserving capability, and easy treatment of contact lines. We also show some preliminary results on the pinch-off process of a pendant drop, where topological changes of the fluid interface are involved. Our method is readily extendable to three dimensions and adaptive meshes.

We introduce a new level set method for representing evolving interfaces. In the case of divergence-free velocity fields, the new method satisfies a conservation principle. Conservation is important for many applications such as modeling two-phase incompressible flow. In the present implementation, the conserved quantity is defined as the integral of a smoothed characteristic function. The new approach embeds level sets into a volume of fluid formulation. The evolution of an approximate signed distance function is governed by a conservation law for its (smoothed) sign. The non-linear level set transport equation is regularized by adding a flux correction term that assures a non-singular Jacobian and penalizes deviations from a distance function. The result is a locally conservative level set method with built-in elliptic redistancing. The continuous model is monolithic in the sense that the level set transport model, the volume of fluid law of mass conservation, and the minimization problem that preserves the approximate distance function property are incorporated into a single equation. There is no need for any extra stabilization, artificial compression, flux limiting, redistancing, mass correction, and other numerical fixes which are commonly used in level set or volume of fluid methods. In addition, there is just one free parameter that controls the strength of regularization and penalization in the model. The accuracy and conservation properties of the monolithic finite element / level set method are illustrated by the results of numerical studies for passive advection of free interfaces., Ergebnisberichte des Instituts für Angewandte Mathematik;586

The current paper presents a large time step wave adding scheme (hereafter LTS-WA) originated in LeVeque's large time Godunov scheme, i.e. constructing numerical schemes via adding strong waves of discontinuity decomposition. Compared with the original scheme, we use a different strategy for wave adding and extend it to multi-dimensional cases. For rarefaction waves, we propose a grid cell decomposition method which can automatically satisfy the entropy condition and avoid nonphysical solutions. We give detailed formulae of the scheme, and make numerical experiments using scalar equations and Euler equations in one and multi-dimensional cases. Numerical results show that, besides the advantage of large time step, the new scheme has a lower numerical dissipation and a higher resolution of shocks and contact discontinuities with the increase of CFL number in certain range.

In this paper, a new type of high-order finite difference and finite volume multi-resolution weighted essentially non-oscillatory (WENO) schemes is presented for solving hyperbolic conservation laws. We only use the information defined on a hierarchy of nested central spatial stencils and do not introduce any equivalent multi-resolution representation. These new WENO schemes use the same large stencils as the classical WENO schemes in [25] , [45] , could obtain the optimal order of accuracy in smooth regions, and could simultaneously suppress spurious oscillations near discontinuities. The linear weights of such WENO schemes can be any positive numbers on the condition that their sum equals one. This is the first time that a series of unequal-sized hierarchical central spatial stencils are used in designing high-order finite difference and finite volume WENO schemes. These new WENO schemes are simple to construct and can be easily implemented to arbitrary high order of accuracy and in higher dimensions. Benchmark examples are given to demonstrate the robustness and good performance of these new WENO schemes.

We develop two local energy-preserving integrators and a global energy-preserving integrator for the general multisymplectic Hamiltonian system. When applied to the 1D and multi-dimensional N-coupled nonlinear Schrodinger equations, the given schemes have the exact preservation of the local/global conservation law and are decoupled in the components ψ n , n = 1 , 2 , … , N , i.e., each of the components can be solved independently. The decoupled feature is significant and helpful for overcoming the computational difficulty of the N-coupled ( N ≥ 3 ) nonlinear Schrodinger equations, especially of the multi-dimensional case. The composition method is employed to improve the accuracy of the schemes in time and the discrete fast Fourier transform is used to reduce the computational complexity. Several numerical experiments are carried out to exhibit the behaviors of the wave solutions. Numerical results confirm the theoretical results.

An alternative formulation of conservative weighted essentially non-oscillatory (WENO) finite difference scheme with the classical WENO-JS weights (Jiang et al. (2013) [6] ) has been successfully used for solving hyperbolic conservation laws. However, it fails to achieve the optimal order of accuracy at the critical points of a smooth function. Here, we demonstrate that the WENO-Z weights (Borges et al. (2008) [1] ) should be employed to recover the optimal order of accuracy at the critical points. Several one- and two-dimensional benchmark problems show the improved performance in terms of accuracy, resolution and shock capturing.

In Fu et al. [23] , a family of high-order targeted ENO (TENO) schemes with a new nonlinear weighting strategy has been proposed. Building upon this strategy, in the current paper, a new class of adaptive TENO schemes for hyperbolic conservation laws is proposed based on three new concepts: (1) a hierarchical voting strategy is proposed to improve the ENO-like stencil selection; (2) the TENO weighting strategy is extended to function as a built-in discontinuity-location detector. Since the reconstruction scheme must not cross any discontinuity, corresponding target schemes are selected from a set of predefined linear schemes which are optimized towards maximum accuracy order or spectral resolution; (3) based on the observation that the cut-off parameter C T in the TENO weighting strategy determines nonlinear dissipation, a C T adaptation strategy is developed to minimize numerical dissipation for high-wavenumber fluctuations while maintaining robustness for shock-capturing. Six-point and eight-point TENO schemes are constructed, and their spectral properties are analyzed by the ADR analysis. A set of benchmark cases is considered to demonstrate the performance of proposed adaptive TENO schemes. Numerical results suggest that the proposed TENO schemes preserve the accuracy order at first and second order critical points and show less numerical dissipation compared with typical WENO schemes. Moreover, the TENO8-NA scheme exhibits good results for both highly compressible flows and nearly incompressible turbulence.

We present a fourth-order accurate finite volume method for the solution of ideal magnetohydrodynamics (MHD). The numerical method combines high-order quadrature rules in the solution of semi-discrete formulations of hyperbolic conservation laws with the upwind constrained transport (UCT) framework to ensure that the divergence-free constraint of the magnetic field is satisfied. A novel implementation of UCT that uses the piecewise parabolic method (PPM) for the reconstruction of magnetic fields at cell corners in 2D is introduced. The resulting scheme can be expressed as the extension of the second-order accurate constrained transport (CT) Godunov-type scheme that is currently used in the Athena astrophysics code. After validating the base algorithm on a series of hydrodynamics test problems, we present the results of multidimensional MHD test problems which demonstrate formal fourth-order convergence for smooth problems, robustness for discontinuous problems, and improved accuracy relative to the second-order scheme., 42 pages, 22 figures, accepted for publication in J. Comp. Phys

A structure-preserving discretization of the shallow-water equations on unstructured spherical grids is introduced. The unstructured grids that we consider have triangular cells with a C-type staggering of variables, where scalar variables are located at centres of grid cells and normal components of velocity are placed at cell boundaries. The staggering necessitates reconstructions and these reconstructions are build into the algorithm such that the resulting discrete equations obey a weighted weak form. This approach, combined with a mimetic discretization of the differential operators of the shallow-water equations, provides a conservative discretization that preserves important aspects of the mathematical structure of the continuous equations, most notably the simultaneous conservation of quadratic invariants such as energy and enstrophy. The structure-preserving nature of our discretization is confirmed through theoretical analysis and through numerical experiments on two different triangular grids, a symmetrized icosahedral grid of nearly uniform resolution and a non-uniform triangular grid whose resolution increases towards the poles.

Structure-preserving discretization of the Rosenbluth–Fokker–Planck equation is still an open question especially for unlike-particle collision. In this paper, a mass-energy-conserving isotropic Rosenbluth–Fokker–Planck scheme is introduced. The structure related to the energy conservation is skew-symmetry in mathematical sense, and the action–reaction law in physical sense. A thermal relaxation term is obtained by using integration-by-parts on a volume integral in the energy moment equation, so the discontinuous Galerkin method is selected to preserve the skew-symmetry. The discontinuous Galerkin method enables ones to introduce the nonlinear upwind flux without violating the conservation laws. Some experiments show that the conservative scheme maintains the mass-energy-conservation only with round-off errors, and analytic equilibria are reproduced only with truncation errors of its formal accuracy.

In this work we present a novel exact Riemann solver for an artificial Equation of State (EoS) based modification of the incompressible Euler equations. Differently from the well known artificial compressibility method, this new approach overcomes the lack of the pressure-velocity coupling by using a suitably designed artificial EoS. The modified set of equations fits into the framework of first order hyperbolic conservation laws. Accordingly, an exact Riemann solver is derived. This new solver has the advantage to avoid a wave pattern violation issue which can affect the exact Riemann problem solution obtained by the standard artificial compressibility approach. The new artificial EoS based Riemann solver can be used as a tool for the definition of the advective Godunov fluxes in a Finite Volume or a discontinuous Galerkin discretization of the incompressible Navier–Stokes (INS) equations. We assess and analyse the new exact Riemann solver on 1D Riemann problems. Its capability and effectiveness are then shown in the context of a high-order accurate discontinuous Galerkin discretization of the INS equations. In particular, we verify the convergence properties, both in space and time, considering two test cases in 2D, namely, the Kovasznay test case and a damped travelling waves test case. Finally, we perform an implicit large eddy simulation of the incompressible turbulent flow over periodic hills where the classic forcing term formulation is modified to deal with variable time steps. Solution comparisons with respect to numerical and experimental results from the literature are given.

We present new Lagrangian discontinuous Galerkin (DG) hydrodynamic methods for compressible flows on unstructured meshes in axisymmetric coordinates. The physical evolution equations for the specific volume, velocity, and specific total energy are discretized using a modal DG method with linear Taylor series polynomials. Two different approaches are used to discretize the evolution equations – the first one is the true volume approach and the second one is the area-weighted approach. For the true volume approach, the DG equations are derived using the true 3D volume that is consistent with the geometry conservation law (GCL). The Riemann velocity at the nodes on the surface of the element, and the corresponding surface forces, are calculated by solving a multidirectional approximate Riemann problem using surfaces areas for axisymmetric coordinates. This true volume approach conserves mass, momentum, and total energy and satisfies the GCL. However, it can not preserve spherical symmetry on an equal-angle polar grid with 1D radial flows. For the area-weighted approach, the DG equations are based on the 2D Cartesian geometry that is rotated about the axis of symmetry using a single, element average radius. With this approach, the Riemann velocity at the nodes on the surface of the element, and the corresponding surface forces, are calculated by solving a multidirectional approximate Riemann problem in 2D Cartesian geometry. This area-weighted approach, in the limit of an infinitesimal mesh size, conserves physical momentum, and physical total energy. The area-weighted approach preserves spherical symmetry on an equal-angle polar grid for 1D radial flows, but it does not satisfy the GCL. A suite of test problems are calculated to demonstrate stable mesh motion, the expected second order accuracy of these methods, and that the new area-weighted DG method preserves spherical symmetry on 1D radial flow problems with equal-angle polar meshes.

Large eddy simulation (LES) has been increasingly used to tackle vortex-dominated turbulent flows. In LES, the quality of the simulation results hinges upon the quality of the numerical discretizations in both space and time. It is in this context we perform a Fourier analysis of several popular methods in LES including the discontinuous Galerkin (DG), finite difference (FD), and compact difference (CD) methods. We begin by reviewing the semi-discrete versions of all methods under consideration, followed by a fully-discrete analysis with explicit Runge–Kutta (RK) time integration schemes. In this regard, we are able to unravel the true dispersion/dissipation behavior of DG and Runge–Kutta DG (RKDG) schemes for the entire wavenumber range using a combined-mode analysis. In this approach, we take into account all eigenmodes in DG and RKDG schemes. The physical-mode is verified to be a good approximation for the asymptotic behavior of these DG schemes in the low wavenumber range. After that, we proceed to compare the DG, FD, and CD methods in dispersion and dissipation properties. Numerical tests are conducted using the linear advection equation to verify the analysis. In comparing different methods, it is found that the overall numerical dissipation strongly depends on the time step. Compact difference (CD) and central FD schemes, in some particular settings, can have more numerical dissipation than the DG scheme with an upwind flux. This claim is then verified through a numerical test using the Burgers' equation.

A new computational model is developed to solve three phase immiscible compressible flow in strongly heterogeneous poroelastic media. Within the context of the iterative coupled formulation based on the fixed stress split algorithm the governing equations are decomposed into three subsystems associated with the geomechanics, hydrodynamics and transport problems. The hydrodynamics involves an overall compressibility which plays essential role in the magnitude of additional source in the pressure equation associated the time derivative of the total mean stress. We show that the additional coefficients in the source can be decomposed into a rock/grain and fluid components involving the formation volume factors of each phase. The flow equations are written in a proper conservative form for application of mixed-hybrid formulation whereas the nonlinear coupling between flow and transport is handled by a proper sequential iterative scheme. A new post-processing approach is proposed to update the Lagrangian porosity based on the freezing of the total mean stress which reproduces the one-way formulation in a straightforward fashion under the stationarity of the total mean stress. The geomechanical subsystem is also discretized by a mixed formulation based on the decomposition of the effective stress into spherical and deviatoric components in order to handle the nearly-incompressible undrained limit. The system of conservation laws for the liquid saturations is rephrased in an alternative form for an appropriate operator splitting method, leading to the appearance of an extra source term arising from the transient behavior of the Lagrangian porosity induced by the geomechanical coupling. Within the framework of a predictor–corrector scheme, the predictor step is discretized by an extended version of a higher order central-upwind scheme whereas the corrector, which captures the influence of rock skeleton deformation upon transport, is captured by a Godunov splitting. The potential of the new computational model is illustrated in numerical simulations of water-flooding and of WAG (Water Alternating Gas) injection problems in poroelastic media. Numerical results illustrate precisely the role of the transient nature of the total stress upon subsidence, rock compaction, oil and gas production and finger grow in primary and secondary recovery regimes.

We present a new efficient method to compute Uehling–Uhlenbeck collision integral for all two-particle interactions in relativistic plasma with drastic improvement in computation time with respect to existing methods. Plasma is assumed isotropic in momentum space. The set of reactions consists of: Moeller and Bhabha scattering, Compton scattering, two-photon pair annihilation, and two-photon pair production, which are described by QED matrix elements. In our method exact energy and particle number conservation laws are fulfilled. Reaction rates are compared, where possible, with the corresponding analytical expressions and convergence of numerical rates is demonstrated.

This paper studies finite element discretizations for three types of time-dependent PDEs, namely heat equation, scalar conservation law and wave equation, which we reformulate as first order systems in a least-squares setting, subject to a space-time conservation constraint (coming from the original PDE). Available piecewise polynomial finite element spaces in ( n + 1 ) -dimensions for functional spaces from the ( n + 1 ) -dimensional de Rham sequence for n = 2 , 3 are used for the implementation of the method. Computational results illustrating the error behavior, iteration counts and performance of block-diagonal and monolithic geometric multigrid preconditioners are presented for the discrete CFOSLS system. The results are obtained from a parallel implementation of the methods for which we report reasonable scalability.

We present in this work a new reconstruction scheme, so-called MUSCL-THINC-BVD scheme, to solve the five-equation model for interfacial two phase flows. This scheme employs the traditional shock capturing MUSCL (Monotone Upstream-centered Schemes for Conservation Law) scheme as well as the interface sharpening THINC (Tangent of Hyperbola for INterface Capturing) scheme as two building-blocks of spatial reconstruction on the BVD (boundary variation diminishing) principle that minimizes the variations (jumps) of the reconstructed variables at cell boundaries, and thus effectively reduces the dissipation error in numerical solutions. The MUSCL-THINC-BVD scheme is implemented to the volume fraction and other state variables under the same finite volume framework, which realizes the consistency among volume fraction and other physical variables. Numerical results of benchmark tests show that the present method is able to capture the material interface as a well-defined sharp jump in volume fraction, and obtain numerical solutions of superior quality in comparison to other existing methods. The proposed scheme is a simple and effective method of practical significance for simulating compressible interfacial multiphase flows.

Shen et al. [1] , [2] proposed an upwind space-time conservation element and solution element (CESE) scheme for 1D and 2D hydrodynamics (HD) in rectangular coordinates, which combined the advantages of CESE and upwind scheme, namely, guaranteed strictly the space–time conservation law as well as captured discontinuities very efficiently. All kinds of upwind schemes can be combined very flexibly for different problems to achieve the perfect combination of CESE and finite volume method (FVM). However, in many physical applications, we need to consider geometries that are more sophisticated. Hence, the main objective of this paper is to extend the upwind CESE scheme to multidimensional magneto-hydrodynamics (MHD) in general curvilinear coordinates by transforming the MHD equations from the physical domain (general curvilinear coordinates) to the computational domain (rectangular coordinates) and the new equations in the computational domain can be still written in the conservation form. For the 3D case, the derivations of some formulas are much more abstract and complex in a 4D Euclidean hyperspace, and some technical problems need to be solved in the debugging process. Unlike in HD, keeping the magnetic field divergence-free for MHD problems is also a challenge especially in general curvilinear coordinates. These are the main obstacles we have overcome in this study. The test results of benchmarks demonstrate that we have successfully extended the upwind CESE scheme to general curvilinear coordinates for both 2D and 3D MHD problems.

High-resolution schemes for conservation laws need to suitably limit the numerical solution near discontinuities, in order to avoid Gibbs oscillations. The solution quality and the computational cost of such schemes strongly depend on their ability to correctly identify troubled-cells, namely, cells where the solution loses regularity. Motivated by the objective to construct a universal troubled-cell indicator that can be used for general conservation laws, we propose a new approach to detect discontinuities using artificial neural networks (ANNs). In particular, we construct a multilayer perceptron (MLP), which is trained offline using a supervised learning strategy, and thereafter used as a black-box to identify troubled-cells. The proposed MLP indicator can accurately identify smooth extrema and is independent of problem-dependent parameters, which gives it an advantage over traditional limiter-based indicators. Several numerical results are presented to demonstrate the robustness of the MLP indicator in the framework of Runge–Kutta discontinuous Galerkin schemes, and its performance is compared with the minmod limiter and the minmod-based TVB limiter.

The mathematical model of shear shallow water flows of constant density is studied. This is a 2D hyperbolic non-conservative system of equations that is mathematically equivalent to the Reynolds-averaged model of barotropic turbulent flows. The model has three families of characteristics corresponding to the propagation of surface waves, shear waves and average flow (contact characteristics). The system is non-conservative: for six unknowns (the fluid depth, two components of the depth averaged horizontal velocity, and three independent components of the symmetric Reynolds stress tensor) one has only five conservation laws (conservation of mass, momentum, energy and mathematical ‘entropy’). A splitting procedure for solving such a system is proposed allowing us to define a weak solution. Each split subsystem contains only one family of waves (either surface or shear waves) and contact characteristics. The accuracy of such an approach is tested on 2D analytical solutions describing the flow with linear with respect to the space variables velocity, and on the solutions describing 1D roll waves. The capacity of the model to describe the full transition scenario as commonly seen in the formation of roll waves: from uniform flow to 1D roll waves, and, finally, to 2D transverse ‘fingering’ of the wave profiles, is shown.

Radiation hydrodynamics simulations based on a single fluid two-temperature model may violate the law of energy conservation, because the governing equations are expressed in a nonconservative formulation. In this study, we maintain the important physical requirements by employing a strategy based on the key concept that mathematical structures associated with conservative and nonconservative equations are preserved, even at the discrete level. To this end, we discretize the conservation laws and transform them using exact algebraic operations. The proposed scheme maintains global conservation errors within the round-off level. In addition, a numerical experiment concerning the shock tube problem suggests that the proposed scheme agrees well with the jump conditions at the discontinuities regulated by the Rankine–Hugoniot relationship. The generalized derivation allows us to employ arbitrary central difference, artificial dissipation, and Runge–Kutta methods.

High-order methods for conservation laws can be highly efficient if their stability is ensured. A suitable means mimicking estimates of the continuous level is provided by summation-by-parts (SBP) operators and the weak enforcement of boundary conditions. Recently, there has been an increasing interest in generalised SBP operators both in the finite difference and the discontinuous Galerkin spectral element framework. However, if generalised SBP operators are used, the treatment of the boundaries becomes more difficult since some properties of the continuous level are no longer mimicked discretely --- interpolating the product of two functions will in general result in a value different from the product of the interpolations. Thus, desired properties such as conservation and stability are more difficult to obtain. Here, new formulations are proposed, allowing the creation of discretisations using general SBP operators that are both conservative and stable. Thus, several shortcomings that might be attributed to generalised SBP operators are overcome (cf. J.~Nordstr\"om and A.~A.~Ruggiu, "On Conservation and Stability Properties for Summation-By-Parts Schemes", \emph{Journal of Computational Physics} 344 (2017), pp. 451--464, and J. Manzanero, G. Rubio, E. Ferrer, E. Valero and D. A. Kopriva, "Insights on aliasing driven instabilities for advection equations with application to Gauss-Lobatto discontinuous Galerkin methods", Journal of Scientific Computing (2017), https://doi.org/10.1007/s10915-017-0585-6).

The objective of this paper is to present a local bounds preserving stabilized finite element scheme for hyperbolic systems on unstructured meshes based on continuous Galerkin (CG) discretization in space. A CG semi-discrete scheme with low order artificial dissipation that satisfies the local extremum diminishing (LED) condition for systems is used to discretize a system of conservation equations in space. The low order artificial diffusion is based on approximate Riemann solvers for hyperbolic conservation laws. In this case we consider both Rusanov and Roe artificial diffusion operators. In the Rusanov case, two designs are considered, a nodal based diffusion operator and a local projection stabilization operator. The result is a discretization that is LED and has first order convergence behavior. To achieve high resolution, limited antidiffusion is added back to the semi-discrete form where the limiter is constructed from a linearity preserving local projection stabilization operator. The procedure follows the algebraic flux correction procedure usually used in flux corrected transport algorithms. To further deal with phase errors (or terracing) common in FCT type methods, high order background dissipation is added to the antidiffusive correction. The resulting stabilized semi-discrete scheme can be discretized in time using a wide variety of time integrators. Numerical examples involving nonlinear scalar Burgers equation, and several shock hydrodynamics simulations for the Euler system are considered to demonstrate the performance of the method. For time discretization, Crank–Nicolson scheme and backward Euler scheme are utilized.

In this paper, we propose a new sixth-order WENO scheme for solving one dimensional hyperbolic conservation laws. The new WENO reconstruction has three properties: (1) it is central in smooth region for low dissipation, and is upwind near discontinuities for numerical stability; (2) it is a convex combination of four linear reconstructions, in which one linear reconstruction is sixth order, and the others are third order; (3) its linear weights can be any positive numbers with requirement that their sum equals one. Furthermore, we propose a simple smoothness indicator for the sixth-order linear reconstruction, this smooth indicator not only can distinguish the smooth region and discontinuities exactly, but also can reduce the computational cost, thus it is more efficient than the classical one.

We present a new predictor-corrector approach to enforcing local maximum principles in piecewise-linear finite element schemes for the compressible Euler equations. The new element-based limiting strategy is suitable for continuous and discontinuous Galerkin methods alike. In contrast to synchronized limiting techniques for systems of conservation laws, we constrain the density, momentum, and total energy in a sequential manner which guarantees positivity preservation for the pressure and internal energy. After the density limiting step, the total energy and momentum gradients are adjusted to incorporate the irreversible effect of density changes. Antidiffusive corrections to bounds-compatible low-order approximations are limited to satisfy inequality constraints for the specific total and kinetic energy. An accuracy-preserving smoothness indicator is introduced to gradually adjust lower bounds for the element-based correction factors. The employed smoothness criterion is based on a Hessian determinant test for the density. A numerical study is performed for test problems with smooth and discontinuous solutions.

A level-set method is presented for solving two-phase flows with soluble surfactant. The Navier–Stokes equations are solved along with the bulk surfactant and the interfacial surfactant equations. In particular, the convection–diffusion equation for the bulk surfactant on the irregular moving domain is solved by using a level-set based diffusive-domain method. A conservation law for the total surfactant mass is derived, and a re-scaling procedure for the surfactant concentrations is proposed to compensate for the surfactant mass loss due to numerical diffusion. The whole numerical algorithm is easy for implementation. Several numerical simulations in 2D and 3D show the effects of surfactant solubility on drop dynamics under shear flow.

We propose a new high-order accurate numerical boundary treatment for solving hyperbolic systems of conservation laws and Euler equations using a Lagrange-remap approach on Cartesian grids in cases of physical boundaries not aligned with the mesh. The method is an adaptation of the Inverse Lax–Wendroff procedure [34] , [35] , [36] , [37] , [38] to the Lagrange-remap approach, which considerably alleviates the algebra. High-order accurate ghost values of conservative variables are imposed using Taylor expansions whose coefficients are found by inverting a (linear or non-linear) system which is well posed in all our examples. For 2D problems, a least-square procedure is added to prevent extrapolation instabilities. The Lagrange-remap formalism also provides a simpler fluid–structure coupling which is also described. Numerical examples are given for the linear case and Euler equations in 1D and 2D.

The fifth order WENO scheme with adaptive order is competent for solving hyperbolic conservation laws, its reconstruction is a convex combination of a fifth order linear reconstruction and three third order linear reconstructions. Note that, on uniform mesh, the computational cost of smoothness indicator for fifth order linear reconstruction is comparable with the sum of ones for three third order linear reconstructions, thus it is too heavy; on non-uniform mesh, the explicit form of smoothness indicator for fifth order linear reconstruction is difficult to be obtained, and its computational cost is much heavier than the one on uniform mesh. In order to overcome these problems, a simple smoothness indicator for fifth order linear reconstruction is proposed in this paper.

This paper discusses a new class of optimized implicit-explicit Runge-Kutta methods for the numerical solution of the dispersive and non-dispersive hyperbolic systems. Optimized implicit-explicit methods are formulated for the better stability and dispersion properties. Moreover, for present methods inversion of the coefficient matrices is not necessary, which makes these methods very attractive in terms of computational cost and complexity. To validate the efficiency of the developed methods, we have solved the one- and two-dimensional dispersive rotating shallow water equations and benchmark problems from acoustics. Computed solutions are also compared with the exact and experimental results available in the literature. The present methods compete well with the existing multi-stage time-integration methods in terms of accurately resolving the physical characteristics for the chosen problems. Furthermore, the computational costs of the proposed methods are significantly lower as compared to the four-stage, fourth-order explicit Runge-Kutta ( RK 4 ) method.

In this paper, a new type of high-order finite volume and finite difference multi-resolution Hermite weighted essentially non-oscillatory (HWENO) schemes are designed for solving hyperbolic conservation laws on structured meshes. Here we only use the information defined on a hierarchy of nested central spatial stencils but do not introduce any equivalent multi-resolution representation, the terminology of multi-resolution HWENO follows that of the multi-resolution WENO schemes (Zhu and Shu, 2018) [29] . The main idea of our spatial reconstruction is derived from the original HWENO schemes (Qiu and Shu, 2004) [19] , in which both the function and its first-order derivative values are evolved in time and used in the reconstruction. Our HWENO schemes use the same large stencils as the classical HWENO schemes which are narrower than the stencils of the classical WENO schemes for the same order of accuracy. Only the function values need to be reconstructed by our HWENO schemes, the first-order derivative values are obtained from the high-order linear polynomials directly. Furthermore, the linear weights of such HWENO schemes can be any positive numbers as long as their sum equals one, and there is no need to do any modification or positivity-preserving flux limiting in our numerical experiments. Extensive benchmark examples are performed to illustrate the robustness and good performance of such finite volume and finite difference HWENO schemes.

This article describes the extension of the conformal symplectic method to solve the damped acoustic wave equation and the elastic wave equations in the framework of the spectral element method. The conformal symplectic method is a variation of conventional symplectic methods to treat non-conservative time evolution problems, which has superior behaviors in long-time stability and dissipation preservation. To reveal the intrinsic dissipative properties of the model equations, we first reformulate the original systems in their equivalent conformal multi-symplectic structures and derive the corresponding conformal symplectic conservation laws. We thereafter separate each system into a conservative Hamiltonian system and a purely dissipative ordinary differential equation system. Based on the splitting methodology, we solve the two subsystems respectively. The dissipative one is cheaply solved by its analytic solution. While for the conservative system, we combine a fourth-order symplectic Nystrm method in time and the spectral element method in space to cover the circumstances in realistic geological structures involving complex free-surface topography. The Strang composition method is adopted thereby to concatenate the corresponding two parts of solutions and generate the completed conformal symplectic method. A relative larger Courant number than that of the traditional Newmark scheme is found in the numerical experiments in conjunction with a spatial sampling of approximately 5 points per wavelength. A benchmark test for the damped acoustic wave equation validates the effectiveness of our proposed method in precisely capturing dissipation rate. The classical Lamb problem is used to demonstrate the ability of modeling Rayleigh wave in elastic wave propagation. More comprehensive numerical experiments are presented to investigate the long-time simulation, low dispersion and energy conservation properties of the conformal symplectic methods in both the attenuating homogeneous and heterogeneous media.

This paper presents a new approach to constraining the eigenvalue range of symmetric tensors in numerical advection schemes based on the flux-corrected transport (FCT) algorithm and a continuous finite element discretization. In the context of element-based FEM-FCT schemes for scalar conservation laws, the numerical solution is evolved using local extremum diminishing (LED) antidi usive corrections of a low order approximation which is assumed to satisfy the relevant inequality constraints. The application of a limiter to antidi usive element contributions guarantees that the corrected solution remains bounded by the local maxima and minima of the low order predictor. The FCT algorithm to be presented in this paper guarantees the LED property for the largest and smallest eigenvalues of the transported tensor at the low order evolution step. At the antidi usive correction step, this property is preserved by limiting the antidi usive element contributions to all components of the tensor in a synchronized manner. The definition of the element-based correction factors for FCT is based on perturbation bounds for auxiliary tensors which are constrained to be positive semidefinite to enforce the generalized LED condition. The derivation of sharp bounds involves calculating the roots of polynomials of degree up to 3. As inexpensive and numerically stable alternatives, limiting techniques based on appropriate approximations are considered. The ability of the new limiters to enforce local bounds for the eigenvalue range is confirmed by numerical results for 2D advection problems., Ergebnisberichte des Instituts für Angewandte Mathematik;564

In this paper, it is shown that three-dimensional stochastic Maxwell equations with multiplicative noise are stochastic Hamiltonian partial differential equations possessing a geometric structure (i.e. stochastic multi-symplectic conservation law), and the energy of system is a conservative quantity almost surely. We propose a stochastic multi-symplectic energy-conserving method for the equations by using the wavelet collocation method in space and stochastic symplectic method in time. Numerical experiments are performed to verify the excellent abilities of the proposed method in providing accurate solution and preserving energy. The mean square convergence result of the method in temporal direction is tested numerically, and numerical comparisons with finite difference method are also investigated.

In this paper, we aim to propose a novel mathematical model for traffic flow and apply a newly developed characteristic particle method to solve the associate governing equations. As compared with the existing non-equilibrium higher-order traffic flow models, the present one is put forward to satisfy the following three conditions: 1. Preserve the equilibrium state in the smooth region. 2. Yield an anisotropic propagation of traffic flow information. 3. Expressed with a conservation law form for traffic momentum. These conditions will ensure a more practical simulation in traffic flow physics: The current traffic will not be influenced by the condition in the behind and result in unambiguous condition across a traffic shock. Through analyses of characteristics, stability condition and steady-state solution adherent to the equation system, it is shown that the proposed model actually conform to these conditions. Furthermore, this model can be cast into its characteristic form which, incorporated with the Rankine-Hugoniot relation, is appropriate to be simulated by the characteristic particle method to obtain accurate computational results.

A moving discontinuous Galerkin finite element method with interface conservation enforcement (MDG+ICE) is developed for solving the compressible Euler equations. The MDG+ICE method is based on the space-time DG formulation, where both flow field and grid geometry are considered as independent variables and the conservation laws are enforced both on discrete elements and element interfaces. The element conservation laws are solved in the standard discontinuous solution space to determine conservative quantities, while the interface conservation is enforced using a variational formulation in a continuous space to determine discrete grid geometry. The resulting over-determined system of nonlinear equations arising from the MDG+ICE formulation can then be solved in a least-squares sense, leading to an unconstrained nonlinear least-squares problem that is regularized and solved by Levenberg-Marquardt method. A number of numerical experiments for compressible flows are conducted to assess the accuracy and robustness of the MDG+ICE method. Numerical results obtained indicate that the MDG+ICE method is able to implicitly detect and track all types of discontinuities via interface conservation enforcement and satisfy the conservation law on both elements and interfaces via grid movement and grid management, demonstrating that an exponential rate of convergence for Sod and Lax-Harden shock tube problems can be achieved and highly accurate solutions without overheating to both double-rarefaction wave and Noh problems can be obtained., 29 pages

ALMA (Anti-symmetric, Large-Moment, Accelerated) is a fast, flexible, and scalable toolkit designed to solve hyperbolic conservation law systems in hybrid supercomputers. This manuscript describes the theoretical background and implementation of ALMA , which uses the anti-symmetric formulation of fluids to obtain simple, robust, and easily paralellizable code. Practical GPU acceleration is realized on entire applications with an overall gain factor of 2 to 4. ALMA also provides a parallel, GPU accelerated sparse solver based on geometric multigrid, capable of diagonalizing linear systems with 239 unknowns. We demonstrate ALMA 's scaling and performance in petascale supercomputers and use standard fluid models to verify the overall approach with canonical benchmark problems.