High-order Spatial Discretization Research Articles

Spatially and temporally high-order dynamic nonlinear variational multiscale methods for generalized Newtonian flows

Non-Newtonian fluids are of interest in industrial sectors, biological problems and other natural phenomena. This work proposes rheologically-dependent, spatially and temporally high-order non-residual stabilized finite element formulations. The accuracy of the methods is assessed by tackling highly-convective time-dependent power-law flows. The spatial approximation uses Lagrangian finite elements up to fourth order. The temporal integration is done via backward differentiation formulas of order one, two and three. A key aspect of our work is using a non-residual orthogonal variational multiscale (VMS) formulation to stabilize dominant convection and to allow equal-order interpolation of velocity and pressure. Our VMS method uses dynamic nonlinear subscales, which have not been tested so far for generalized Newtonian fluids. In this work, the use of high-order temporal discretizations for the subscale components is systematically evaluated. Numerical experiments consider the flow over a confined cylinder for Reynolds numbers between 40 and 400 and power-law indices between 0.4 and 1.8. Numerical testing demonstrates the method to be stable in all combinations of spatial and temporal orders. Our results show that using high-order spatial discretizations more accurately approximates boundary layers and viscosity fields. Moreover, higher temporal orders allow using larger time steps while still capturing highly dynamic behaviors with better resolution in frequency spectra.

Implementation of a high-order spatial discretization into a finite volume solver: Applications to turbomachinery test cases using an eddy-viscosity turbulence closure

In this study, the implementation of a high-order spatial discretization method into a Finite Volume solver is presented. Specific emphasis is put on the analysis of the performance over selected turbomachinery test cases. High-order numerical discretization is achieved by adopting the cell-centered Least-Square reconstruction, which is implemented in the in-house solver HybFlow. The validation of the adopted methodology is performed by assessing the solution of a turbulent flat plate with zero pressure gradient, using a eddy-viscosity transitional model. The test case also evidences the effect of the discretization of gradient-based source terms when a high-order reconstruction methodology is used. In the second part of the paper, the solver is used for the solution of relevant two-dimensional turbomachinery test cases, assessing the impact of 2nd and 3rd order reconstruction on the prediction of the aerodynamics and the heat transfer for respectively a low-pressure blade and a high-pressure turbine vane. It is shown how a high-order reconstruction allows for obtaining a better prediction of turbomachinery aerodynamics, with lower number of elements. The benefits over heat transfer predictions in high Reynolds number conditions are instead limited to the reduction of heat transfer coefficient spikes in under-resolved regions of the blade. Eventually, the methodology is validated for a three-dimensional low-pressure turbine cascade with realistic boundary layer inflow conditions.

Convergence towards High-Speed Steady States Using High-Order Accurate Shock-Capturing Schemes

Creating time-marching unsteady governing equations for a steady state in high-speed flows is not a trivial task. Residue convergence in time cannot be achieved when using most low- and high-order spatial discretization schemes. Recently, high-order, weighted, essentially non-oscillatory schemes have been specially designed for steady-state simulations. They have been shown to be capable of achieving machine precision residues when simulating the Euler equations under canonical coordinates. In the present work, we review these schemes and show that they can also achieve machine residues when simulating the Navier–Stokes equations under generalized coordinates. This is carried out by considering three supersonic flows of perfect fluids, namely the flow upstream a cylinder, the flow over a blunt wedge, and the flow over a compression ramp.

Algebraically stable SDIRK methods with controllable numerical dissipation for first/second-order time-dependent problems

In this paper, a family of four-stage singly diagonally implicit Runge-Kutta methods are proposed to solve first-/second-order time-dependent problems, exhibiting the following numerical properties: fourth-order accuracy in time, unconditional stability, controllable numerical dissipation, and adaptive time step selection. The BN-stability condition is employed as a constraint to optimize parameters in the Butcher table, having significant benefits, and hence is recommended for nonlinear dynamics problems in contrast to existing methods. Numerical examples involving both first- and second-order linear/nonlinear dynamics problems validate the proposed method, and numerical results reveal that the proposed methods are free from the order reduction phenomenon when applied to nonlinear dynamics problems. The performance of adaptive time-stepping using the embedded scheme is further illustrated by the phase-field modeling problem. Additionally, the advantages and disadvantages of three-stage third-order accurate algebraically stable methods are discussed. The proposed high-order time integration can be readily integrated into high-order spatial discretization methods, such as the high-order spectral element method employed in this paper, to obtain high-order discretization in space and time dimensions.

Validation of wall boundary conditions for simulating complex fluid flows via the Boltzmann equation: Momentum transport and skin friction

The influence and validity of wall boundary conditions for non-equilibrium fluid flows described by the Boltzmann equation remains an open problem. The substantial computational cost of directly solving the Boltzmann equation has limited the extent of numerical validation studies to simple, often two-dimensional, flow problems. Recent algorithmic advancements for the Boltzmann–Bhatnagar–Gross–Krook equation introduced by the authors [Dzanic et al., J. Comput. Phys. 486, 112146 (2023)], consisting of a highly efficient high-order spatial discretization augmented with a discretely conservative velocity model, have made it feasible to accurately simulate unsteady three-dimensional flow problems across both the rarefied and continuum regimes. This work presents a comprehensive evaluation and validation of wall boundary conditions across a variety of flow regimes, primarily for the purpose of exploring their effects on momentum transfer in the low Mach limit. Results are presented for a range of steady and unsteady wall-bounded flow problems across both the rarefied and continuum regimes, from canonical two-dimensional laminar flows to unsteady three-dimensional transitional and turbulent flows, the latter of which are the first instances of wall-bounded turbulent flows computed by directly solving the Boltzmann equation. We show that approximations of the molecular gas dynamics equations can accurately predict both non-equilibrium phenomena and complex hydrodynamic flow instabilities and show how spatial and velocity domain resolution affect the accuracy. The results indicate that an accurate approximation of particle transport (i.e., high spatial resolution) is significantly more important than particle collision (i.e., high velocity domain resolution) for predicting flow instabilities and momentum transfer consistent with that predicted by the hydrodynamic equations and that these effects can be computed accurately even with very few degrees of freedom in the velocity domain. These findings suggest that highly accurate spatial schemes (e.g., high-order schemes) are a promising approach for solving molecular gas dynamics for complex flows and that the direct solution of the Boltzmann equation can be performed at a reasonable cost when compared to hydrodynamic simulations at the same level of resolution.

TARS: A parallel tetrahedral discontinuous finite element code for the solution of the discrete ordinates neutron transport equation

In this study, we introduce TARS, a neutron transport code based on the discrete-ordinates discontinuous finite element (SN-DFEM) method, which tackles existing challenges in the field by providing accurate and efficient solutions to the neutron transport equation. The use of DFEM enables higher-order spatial discretization, improving numerical stability, while parallelizing transport sweeping reduces computational time and resource requirements. TARS' performance was evaluated using the Takeda benchmark problems. The results showed that TARS generates highly accurate solutions, with eigenvalue errors no more than 90 pcm compared to MC solutions. The relative differences of region-averaged fluxes were below 1% in most regions, except those with severe flux gradients. Parallel calculation efficiency analysis revealed a PCE of over 70% even at 1024 CPU cores. However, hybrid MPI and OpenMP parallel programming did not outperform pure MPI programming, suggesting that the latter might be more suitable for TARS.

High order spatial discretization for variational time implicit schemes: Wasserstein gradient flows and reaction-diffusion systems

We design and compute first-order implicit-in-time variational schemes with high-order spatial discretization for initial value gradient flows in generalized optimal transport metric spaces. We first review some examples of gradient flows in generalized optimal transport spaces from the Onsager principle. We then use a one-step time relaxation optimization problem for time-implicit schemes, namely generalized Jordan-Kinderlehrer-Otto schemes. Their minimizing systems satisfy implicit-in-time schemes for initial value gradient flows with first-order time accuracy. We adopt the first-order optimization scheme ALG2 (Augmented Lagrangian method) and high-order finite element methods in spatial discretization to compute the one-step optimization problem. This allows us to derive the implicit-in-time update of initial value gradient flows iteratively. We remark that the iteration in ALG2 has a simple-to-implement point-wise update based on optimal transport and Onsager's activation functions. The proposed method is unconditionally stable for convex cases. Numerical examples are presented to demonstrate the effectiveness of the methods in two-dimensional PDEs, including Wasserstein gradient flows, Fisher–Kolmogorov-Petrovskii-Piskunov equation, and two and four species reversible reaction-diffusion systems.

High-order methods for the simulation of unsteady counterflow flames subject to stochastic forcing of large amplitude

Unsteady counterflows are employed to understand and model the effect of turbulence on flames. We present a novel numerical approach for the simulation of one-dimensional unsteady counterflow flames with fourth order spatial discretization and up to fourth order time discretization. The approach couples a three-stage Lobatto IIIa formula for boundary value problems and variable-order, variable time step size Backward Differentiation Formulas for time integration. The framework is explained in detail, its computational performance is analysed, and its use is demonstrated for the case of stochastic forcing of premixed counterflow flames, whereby the imposed rate of strain is a multi-scale lognormal discrete random process with exponential autocorrelation. High-order spatial and temporal discretization make the approach well-suited for the accurate and computationally efficient simulation of the effect of turbulence on flames, which are characterised by large amplitude stochastic fluctuations of the local rate of strain.

A filtering monotonization approach for DG discretizations of hyperbolic problems

We introduce a filtering technique for Discontinuous Galerkin approximations of hyperbolic problems. Following an approach already proposed for the Hamilton-Jacobi equations by other authors, we aim at reducing the spurious oscillations that arise in presence of discontinuities when high order spatial discretizations are employed. This goal is achieved using a filter function that keeps the high order scheme when the solution is regular and switches to a monotone low order approximation if it is not. The method has been implemented in the framework of the $deal.II$ numerical library, whose mesh adaptation capabilities are also used to reduce the region in which the low order approximation is used. A number of numerical experiments demonstrate the potential of the proposed filtering technique.

A high-order convergence analysis for semi-Lagrangian scheme of the Burgers' equation

<abstract><p>In this article, we provide a comprehensive convergence and stability analysis of a semi-Lagrangian scheme for solving nonlinear Burgers' equations with a high-order spatial discretization. The analysis is for the iteration-free semi-Lagrangian scheme comprising the second-order backward finite difference formula (BDF2) for total derivative and the fourth-order central finite difference for diffusion term along the trajectory. The main highlight of the study is to thoroughly analyze the order of convergence of the discrete $ \ell_2 $-norm error $ \mathcal{O}(h^2+\triangle x^4+ \triangle x^{p+1}/h) $ by managing the relationship between the local truncation errors from each discretization procedure and the interpolation properties with a symmetric high-order discretization of the diffusion term. Furthermore, stability is established by the uniform boundedness of the numerical solution using the discrete Grönwall's Lemma. We provide numerical examples to support the validity of the theoretical convergence and stability analysis for the propounded backward semi-Lagrangian scheme.</p></abstract>

Mathematical Modeling of Real-Time Systems Using Heun and Piecewise Methods

It is often said that mathematical modeling is an implementation of mathematics in real-world problems with the aim of better understanding them, so we can say that mathematical modeling is linked to the solution of problems. Some of the essential principles and procedures of mathematical modeling are discussed using formulas and equations. We investigate the stability and convergence characteristics and demonstrate the suitability of different mathematical methods in a set of numerical examples. The described methods in our paper are the best choices for the simulation of linear phenomena and are more efficient for use with high-order spatial discretization. We emphasized the importance of mathematical modeling technologies used in computational tools. Our study shows that these new methods are more stable with lower errors.

Positivity-preserving and energy-dissipative finite difference schemes for the Fokker–Planck and Keller–Segel equations

Abstract In this work we introduce semi-implicit or implicit finite difference schemes for the continuity equation with a gradient flow structure. Examples of such equations include the linear Fokker–Planck equation and the Keller–Segel equations. The two proposed schemes are first-order accurate in time, explicitly solvable, and second-order and fourth-order accurate in space, which are obtained via finite difference implementation of the classical continuous finite element method. The fully discrete schemes are proved to be positivity preserving and energy dissipative: the second-order scheme can achieve so unconditionally while the fourth-order scheme only requires a mild time step and mesh size constraint. In particular, the fourth-order scheme is the first high order spatial discretization that can achieve both positivity and energy decay properties, which is suitable for long time simulation and to obtain accurate steady state solutions.

On the properties of high-order least-squares finite-volume schemes

High-order finite-volume schemes based on polynomial least-squares methods are an active research topic for the discretization of hyperbolic equations as they allow to obtain high-order spatial discretization schemes in arbitrary meshes. However, few studies have analyzed their performance in good-quality/near-to-uniform meshes, which are commonly used as a meshing strategy in zones where turbulent effects are important. In this paper, the theoretical numerical properties of commonly used least-squares (LSQ) k-exact high-order finite volume schemes are studied in one-dimensional and in several two-dimensional meshes (with some remarks regarding their properties in three-dimensional meshes). These results are compared to those obtained using fully-constrained polynomial reconstructions only compatible with structured meshes. The numerical properties of the schemes are investigated through the von Neumann analysis methodology applied to the one-dimensional and two-dimensional finite volume formulation, including temporal discretization errors. This analysis is also extended to non-uniform and unstructured two-dimensional meshes. At last, the schemes are tested with several numerical experiments using the linear advection, the Euler equations and the Navier-Stokes equations. The analysis of both studies yields similar conclusions regarding the numerical errors and stability of the different studied schemes showing that the high-order least-squares finite volume schemes yield stable and robust results across different uniform and non-uniform unstructured meshes. However, their performance is heavily degraded in simulations with low mesh resolution compared to schemes specially catered to structured meshes. On the other hand, the latter schemes lack stability and robustness in general structured meshes and its formulation cannot be straightforwardly extended to unstructured meshes. Moreover, this work shows that the use of weighted-LSQ can drastically improve the results of LSQ schemes in under-resolved simulations.

Development of a Balanced Adaptive Time-Stepping Strategy Based on an Implicit JFNK-DG Compressible Flow Solver

A balanced adaptive time-stepping strategy is implemented in an implicit discontinuous Galerkin solver to guarantee the temporal accuracy of unsteady simulations. A proper relation between the spatial, temporal and iterative errors generated within one time step is constructed. With an estimate of temporal and spatial error using an embedded Runge-Kutta scheme and a higher order spatial discretization, an adaptive time-stepping strategy is proposed based on the idea that the time step should be the maximum without obviously influencing the total error of the discretization. The designed adaptive time-stepping strategy is then tested in various types of problems including isentropic vortex convection, steady-state flow past a flat plate, Taylor-Green vortex and turbulent flow over a circular cylinder at {Re}=3,900. The results indicate that the adaptive time-stepping strategy can maintain that the discretization error is dominated by the spatial error and relatively high efficiency is obtained for unsteady and steady, well-resolved and under-resolved simulations.

Generalized Petrov-Galerkin time finite element weighted residual methodology for designing high-order unconditionally stable algorithms with controllable numerical dissipation

In this paper, a novel stabilized time-weighted residual methodology under the umbrella of Petrov-Galerkin time finite element formulation is developed to design a generalized computational framework, which permits unconditionally stable, high-order time accuracy, and features with controllable numerical dissipation for solving transient first-order systems. Various unconditionally stable (A/L-stable) algorithms can be readily obtained in the proposed framework having not only high-order accuracy but also controllable numerical dissipation in the high frequency. Quadratic and cubic basis functions are utilized to illustrate the specific design process of the proposed method, which consequently ends up with numerous third-/fifth-order time accurate algorithms, QUAD3(γ,ρ∞) and CUBE5(γ,ρ∞), with controllable numerical dissipation. Comparing with the well-known implicit Runge-Kutta (RK) family of algorithms, these newly developed QUAD3(γ,ρ∞) and CUBE5(γ,ρ∞) can (a) recover the Radau IIA3/RK3 and Radau IIA5/RK5 schemes by certain selection of algorithmic parameters (γ,ρ∞); (b) obtain numerous new algorithms with improved solution accuracy that is superior to the RK family of algorithms, and (c) has similar computational efficiency as that of the implicit RK family of algorithms. Several single/multi-degree of freedom (SDOF and MDOF) problems are investigated to validate the proposed developments. In addition, it is worth noting that the proposed methodology can also use high-order (not limited to the quadratic and cubic) basis functions to design more advanced schemes and can be integrated with high-order spatial discretization methods in the solution of space-time PDEs, such as isogeometric methods, SEM, Discontinuous Galerkin (DGM), p-versions FEM, etc.

Delta-form-based method of solving high order spatial discretization schemes for neutron transport

Delta-form-based methods for solving high order spatial discretization schemes are introduced into the reactor SN transport equation. Due to the nature of the delta-form, the final numerical accuracy only depends on the residuals on the right side of the discrete equations and have nothing to do with the parts on the left side. Therefore, various high order spatial discretization methods can be easily adopted for only the transport term on the right side of the discrete equations. Then the simplest step or other robust schemes can be adopted to discretize the increment on the left hand side to ensure the good iterative convergence. The delta-form framework makes the sweeping and iterative strategies of various high order spatial discretization methods be completely the same with those of the traditional SN codes, only by adding the residuals into the source terms. In this paper, the flux limiter method and weighted essentially non-oscillatory scheme are used for the verification purpose to only show the advantages of the introduction of delta-form-based solving methods and other high order spatial discretization methods can be also easily extended to solve the SN transport equations. Numerical solutions indicate the correctness and effectiveness of delta-form-based solving method.

HTR-1.2 solver: Hypersonic Task-based Research solver version 1.2

We present an updated version of the open-source Hypersonics Task-based Research (HTR) solver for hypersonic aerothermodynamics. The solver, whose first version was presented in Di Renzo et al. (2020), is designed for direct numerical simulation (DNS) of canonical hypersonic flows at high Reynolds numbers in which thermo-chemical effects induced by high temperatures are relevant. The solver relies on high-order spatial discretization on structured meshes and efficient time integrators for stiff systems within the Regent/Legion software stack, which makes the code highly portable and scalable in CPU and GPU-based supercomputers. The new version herein presented includes several optimizations and new tools for data analysis, along with novel user option for hybrid skew-symmetric/targeted essentially non-oscillatory numerics, to offer higher computational efficiency and lower numerical dissipation at moderate Mach numbers, inclusion of a new combustion mechanism for methane and oxygen, and new recycling–rescaling inlet boundary conditions targeted to the simulation of fully developed turbulent boundary layers. New version program summaryProgram Title: Hypersonics Task-based Research solverCPC Library link to program files:https://doi.org/10.17632/9zsxjtzfr7.2Developer’s repository link:https://github.com/stanfordhpccenter/HTR-solver.gitLicensing provisions: BSD 2-clauseProgramming language: Regent, C++, and CUDAJournal Reference of previous version: Di Renzo, M., Fu, L., & Urzay, J. (2020). HTR solver: An open-source exascale-oriented task-based multi-GPU high-order code for hypersonic aerothermodynamics. Computer Physics Communications, 107262.Does the new version supersede the previous version?: YesReasons for the new version: Release of new featuresSummary of revisions:•New optional sixth-order hybrid scheme has been implemented (activated by the flag “hybridScheme” in the input file). The scheme combines the energy-preserving properties of a sixth-order skew-symmetric central difference scheme [1] in smooth flow regions with the shock-capturing properties of a sixth-order targeted essentially non-oscillatory (TENO) scheme at points where shocks are involved. The switch between the two schemes is controlled by a TENO sensor whose cutoff value is adapted based on the maximum value of a modified Ducros sensor [2] across the reconstruction stencil. If the flag “hybridScheme” is set to false, the numerical scheme will revert to the TENO6-A scheme released in the previous version of the solver [3];•New recycling–rescaling inflow boundary conditions [4,5] for the simulation of turbulent compressible boundary layers are now available to the user;•Support for Legion tracing, which significantly improves the strong scalability of the solver, has been implemented;•A diagnostic tool to monitor the time evolution of the flow variables in a subvolume of the computational domain is now available;•A single-step chemistry mechanism for methane/oxygen combustion has been added to the mixtures handled by the HTR solver;•Sample scripts for strong scaling have been added to the “testcases” directory;•Unit test and regression test suites have been added to the repository;•The input file scheme has been modified in order to reduce verbosity and increase flexibility in specifying the boundary conditions and type of gas mixture;•Hyperbolic sine stretching functions have been made available to users during the grid generation process;•Computationally intensive tasks have been ported to C++ and CUDA in order to achieve higher efficiency on all hardware;•Several optimizations of the tasks body and mapper have been implemented in order to increase the computational efficiency and reduce the memory footprint.Nature of problem: This code solves the Navier–Stokes equations at hypersonic Mach numbers including finite-rate chemistry for dissociating air and multicomponent transport. The solver is designed for direct numerical simulations (DNS) of transitional and turbulent hypersonic turbulent flows under high-enthalpy conditions, and it accounts for thermochemical effects (vibrational excitation and chemical dissociation).Solution method: This code uses a low-dissipation sixth-order schemes for the spatial discretization of the conservation equations on Cartesian stretched meshes. Time advancement is carried out by either an explicit method if chemistry is slow, hence not introducing additional stiffness, or by an operator-splitting algorithm whereby chemical production rates are handled implicitly.Additional comments including restrictions and unusual features: The HTR solver builds on the runtime Legion [6,7] and is written in the programming language Regent [8,9] developed at Stanford University. Instructions for installation of the components are provided in the README file enclosed with the HTR solver and in the Legion repository [6].

Structure aware Runge\u2013Kutta time stepping for spacetime tents

We introduce a new class of Runge–Kutta type methods suitable for time stepping to propagate hyperbolic solutions within tent-shaped spacetime regions. Unlike standard Runge–Kutta methods, the new methods yield expected convergence properties when standard high order spatial (discontinuous Galerkin) discretizations are used. After presenting a derivation of nonstandard order conditions for these methods, we show numerical examples of nonlinear hyperbolic systems to demonstrate the optimal convergence rates. We also report on the discrete stability properties of these methods applied to linear hyperbolic equations.

Fully-Discrete Analysis of High-Order Spatial Discretizations with Optimal Explicit Runge–Kutta Methods

High-order unstructured methods have become a popular choice for the simulation of complex unsteady flows. Flux reconstruction (FR) is a high-order spatial discretization method, which has been found to be particularly accurate for scale-resolving simulations of complex phenomena. In addition, it has been shown to provide sufficient dissipation for implicit large-eddy simulation (ILES). In conjunction with an FR discretization, an appropriate temporal scheme must be chosen. A common choice is explicit schemes due to their efficiency and ease of implementation. However, these methods usually require a small time-step size to remain stable. Recently, the development of optimal explicit Runge–Kutta (OERK) schemes has enabled stable simulations with larger time-step sizes. Hence, we analyze the fully-discrete properties of the FR method with OERK temporal schemes. We show results for first, second, third, fourth and eighth-order OERK schemes. We observe that OERK schemes modify the spectral behaviour of the semidiscretization. In particular, dissipation decreases in the region of high wavenumbers. We observe that higher-order OERK schemes require a smaller time step than the low-order schemes. However, they follow the dispersion relations of the FR scheme for a larger range of wavenumbers. We validate our analysis with simple advection test cases. It was observed that first and second-degree temporal schemes introduce a relatively large amount of error in the solutions. A one-dimensional ILES test case showed that, as long as the time-step size is not in the vicinity of the stability limit, results are generally similar to classical RK schemes.

A directional ghost-cell immersed boundary method for incompressible flows

This paper describes an efficient ghost-cell immersed boundary method for incompressible flow simulations. The method employs a locally directional extrapolation scheme along all discretization directions for the ghost values. Additionally, it involves fictitious discrete boundary forcing terms instead of the ghost values in the governing equations. In this way, the boundary is represented more accurately than in prior IBMs and it is possible to fulfill the divergence-free condition. When combined with high-order spatial discretization schemes, the IBM order is reduced locally near the immersed boundary in a step-wise manner. In this way, the method delivers more compact stencils and is able to deal with sharp interfaces. By handling the distance from the ghost point to the boundary carefully, the proposed method delivers lower truncation errors than standard IBM, with clean (persistent) convergence rate and enhanced stability. The parallel implementation of this approach is straightforward. Its accuracy has been checked by considering a variety of test cases, including irregular, three-dimensional, and moving boundaries. The local accuracy for the proposed method is formally second order, and is measured to be close to this value using numerical tests.