Revision as of 16:38, 4 August 2020

Code Description

This page is dedicated to a short presentation of the three low-Mach number codes used for the rest of this benchmarl: YALES2, DINO and Nek5000.

Since these codes have already been the subject of many publications and are not completely new, only the most relevant features are discussed in what follows, with suitable references for those readers needing more details.

YALES2

YALES2 is a massively parallel multiphysics platform^[1][2] developed since 2009 by Moureau, Lartigue and co-workers at CORIA (Rouen, France). It is dedicated to the high-fidelity simulation of low-Mach number flows in complex geometries. It is based on the Finite Volumes formulation of the Navier-Stokes equations and it can solve both non-reacting and reacting flows. It can actually solve various physical problems thanks to its original structure which is composed of a main numerical library accompanied with tens of dedicated solvers (for acoustics, multiphase flows, heat transfer, radiation\ldots) which can be coupled with one another. YALES2 relies on unstructured meshes and a fully parallel dynamic mesh adaptation technique to improve the resolution in physically-relevant zones to mitigate the computational cost^[3]. As a result it can easily handle meshes composed of billions of tetrahedra, thus enabling the Direct Numerical Simulation of laboratory and semi-industrial configurations. It is now composed of nearly 500,000 lines of object-oriented Fortran and the parallelism is currently ensured by a pure MPI paradigm, although a hybrid OpenMP/MPI as well as a GPU version are under development.

As most low-Mach number codes, the time-advancement is based on a projection-correction method following the pioneering work of^[4]. The prediction step uses a method which is a blend between a 4th-order Runge-Kutta method and a 4th-order Lax-Wendroff-like method^[5], combined with a 4th-order node-based centered finite-volume discretization of the convective and diffusive terms. Moreover, to improve the performance of the correction step, the pressure of the previous iteration is included in the prediction step to limit the splitting errors. The correction step is required to ensure mass conservation, using a pressure that arises from a Poisson equation. This is performed numerically by solving a linear system on the pressure at each node of the mesh thanks to a dedicated in-house version of the deflated conjugate gradient algorithm, which has been optimized for solving elliptic equations on massively parallel machines^[6].

When considering multispecies and non-isothermal flows, the extension of the classical projection method proposed by Pierce et al.^[7] is used to account for variable-density flows. All the thermodynamic and transport properties are provided by the Cantera software^[8], which has been fully re-implemented in Fortran to avoid any performance issues. Each species thermodynamic property is specified by 5th-order polynomials on two temperature ranges (below and above 1000 K). The mixture is supposed to be both thermally and mechanically perfect. Regarding transport properties, several type of models are implemented in YALES2: 1) the default approach is based on the computation of transport properties for each species (using tabulated molecular potentials), then combining those to obtain mixture-averaged coefficients for viscosity, conductivity and diffusion velocities (Hirschfelder and Curtiss approximation^[9]; 2) alternatively, simplified laws (for example the Sutherland law^[10] can be used for the viscosity, while fixed values for Prandtl and Schmidt numbers allow the computation of conductivity and diffusion coefficients; and 3) a mix of both approaches, for example computing viscosity and conductivity with a mixture-averaged approximation and then imposing a Lewis number for each species. This last approach has been retained in the present benchmark.

Finally, the source terms used in the reacting simulations are modeled with an Arrhenius law with the necessary modifications needed to take into account three-body or pressure-dependent reactions.

From a numerical point of view, it must be noticed that both the diffusion and reaction processes occur at time scales which can be orders of magnitude smaller than the convective time scale. Solving these phenomena with explicit methods would thus drastically limit the global timestep of the whole simulation and induce an overwhelming CPU cost. To mitigate these effects when dealing with multi-species reacting flows, the classical operator splitting technique is used. The diffusion process is solved with a fractional timestep method inside each convective iteration, each substep being limited by a Fourier condition to ensure stability. This method gets activated only when the mesh is very fine (typically when performing DNS with very diffusive species like ${\displaystyle H_{2}}$ or ${\displaystyle H}$, as in the present project); otherwise an explicit treatment is sufficient. The chemical source terms are then integrated with a dedicated stiff solver, namely the CVODE library from SUNDIALS^[11][12][13]. To that purpose, each control volume is considered as an isolated reactor with both constant pressure and enthalpy.

Using the stiff integration technique results in a very strong load imbalance between the various regions of the flow: the fresh gases are solved in a very small number of integration steps while the inner flame region may require tens or even hundreds of integration steps. To overcome this difficulty, a dedicated MPI dynamic scheduler based on a work-sharing algorithm ensures global load-balancing.

DINO

DINO is a 3-D DNS code used for incompressible or low-Mach number flows, the latter approach being used in this project. The development of DINO started in the group of D. Thevénin (Univ. of Magdeburg) in 2013. DINO is a Fortran-90 code, written on top of a 2-D pencil decomposition to enable efficient large-scale parallel simulations on distributed-memory supercomputers by coupling with the open-source library 2-DECOMP&FFT^[14]. The code offers different features and algorithms in order to investigate different physicochemical processes. Spatial derivative are computed by default using sixth-order central finite differences. Time integration relies on several Runge-Kutta solvers. In what follows, an explicit 4th-order Runge-Kutta approach has been used. A 3rd-order semi-implicit Runge-Kutta integration can be activated as needed, when considering stiff chemistry. In this case, non-stiff terms are still computed with explicit Runge-Kutta, while the PyJac package^[15] is used to integrate in an implicit manner all chemistry terms with an analytical Jacobian computation. All thermodynamic, chemical and transport properties are computed using the open-source library Cantera 2.4.0^[16]. The transport properties can be computed based on three different models: 1) Constant Lewis numbers, 2) mixture-averaged, 3) full multicomponent diffusion, by coupling either again with Cantera or with the open-source library EGlib^[17]. The Poisson equation is solved using Fast Fourier Transform (FFT) for periodic as well as for non-periodic boundary conditions, relying in the latter case on an in-house pre- and post-processing technique. The I/O operations are implemented using two different approaches: (1) binary MPI-I/O using 2-DECOMP&FFT for check-points and restart files; (2) HDF5 files used for analysis and visualization.

Multi-phase flows can be simulated in DINO using resolved or non-resolved (point) particles and droplets using a Lagrangian approach^[18][19]. Complex boundaries are represented by a novel second-order immersed boundary method implementation (IBM) based on a directional extrapolation scheme^[20]. More details about the implemented algorithms can be found in particular in^[21]. Since DINO has been developed as a multi-purpose code for analyzing many different reacting and non-reacting flows^[22][23][24], a detailed verification and validation is obviously essential.

Nek5000

This spectral element low-Mach number reacting flow solver is based on the highly-efficient open-source solver Nek5000^[25] extended by a plugin developed at ETH implementing a high-order splitting scheme for low-Mach number reacting flows^[26]. The spectral element method (SEM) is a high-order weighted residual technique for spatial discretization that combines the accuracy of spectral methods with the geometric flexibility of the finite element method allowing for accurately representation of complex geometries^[27]. The computational domain is decomposed into ${\displaystyle E}$ conforming elements, which are quadrilaterals (in 2-D) or hexahedra (in 3-D) that conform to the domain boundaries. Within each element, functions are expanded as ${\displaystyle N}$th-order polynomials so that resolution can be increased either by decreasing the element size (${\displaystyle h}$-type refinement) or by increasing the polynomial order (${\displaystyle p}$-type refinement; typically ${\displaystyle N=7-15}$). The grids can be unstructured and allow for static local refinement, while adaptive mesh refinement has been recently developed^[28]. By casting the polynomial approximation in tensor-product form, the differential operators on ${\displaystyle N^{3}}$ gridpoints per element can be evaluated with only ${\displaystyle O(N^{4})}$ work and ${\displaystyle O(N^{3})}$ storage.

The principal advantage of the SEM is that convergence is exponential in ${\displaystyle N}$, yielding minimal numerical dispersion and dissipation, so that significantly fewer grid points per wavelength are required in order to accurately propagate a turbulent structure over the extended time required in high Reynolds number simulations. Nek5000 uses locally structured basis coefficients (${\displaystyle N\times N\times N}$ arrays), which allow direct addressing and tensor-product-based derivative evaluation that can be cast as efficient matrix-matrix products involving ${\displaystyle N^{2}}$ operators applied to ${\displaystyle N^{3}}$ data values for each element. As a result, data movement per grid point is the same as for low-order methods. It uses scalable domain-decomposition-based iterative solvers with efficient preconditioners. Communication is based on the Message Passing Interface (MPI) standard, and the code has proven scalability to over one million ranks. Nek5000 provides balanced I/O latency among all processors and reduces the overhead or even completely hides the I/O latency by using dedicated I/O communicators in the optimal case.

Time advancement is performed using the splitting scheme proposed in^[26] to decouple the highly non-linear and stiff thermochemistry (species and energy governing equations) from the hydrodynamic system (continuity and momentum). Species and energy equations are integrated without further splitting using the implicit stiff integrator solver CVODE from the SUNDIALS package^[11][12][13] that uses backward differentiation formulas (BDF). The continuity and momentum equations are integrated using a second- or third-order semi-implicit formulation (EXT/BDF) treating the non-linear advection term explicitly^[27]. The thermodynamic properties, detailed chemistry, and transport properties are provided by optimized subroutines compatible with Chemkin^[29].

The reacting flow solver can handle complex time-varying geometries and has been used for instance to simulate laboratory-scale internal combustion engines^[30][31]. It can account for conjugate fluid-solid heat transfer and detailed gas phase as well as surface kinetics^[32].

General comments on the codes

The three aforementioned solvers are all unsteady, high-fidelity codes based on the low-Mach number formulation of the Navier-Stokes equations. They can perform both Direct Numerical Simulations or Large Eddy Simulations of reacting flows, only DNS being considered here. However, they differ in a certain number of points, mainly from the numerical point of view. The aim of this section is to emphasize those major differences. First, YALES2 is an unstructured code, designed to handle any type of elements; its main application field pertains to LES of industrially relevant flows, though DNS is possible as well. On the other hand, both DINO and Nek5000 are mostly dedicated to DNS of configurations found in fundamental research. Both DINO and Nek5000 can only deal with quads or hexas, with the major difference that DINO is based on a structured connectivity while Nek5000 can use unstructured meshes (pavings). As a consequence, both DINO and Nek5000 employ higher-order numerical schemes compared to YALES2, limited at best to a 4th-order scheme. All codes rely on dedicated libraries to compute the thermo-chemical properties of the flow, either Chemkin, Cantera, or in-house versions of those. Moreover, they also rely at least to some extent on external software to perform the temporal integration of the stiff chemical source terms. Regarding the Poisson equation for pressure which must be solved by all codes, DINO relies on a spectral formulation by performing direct and inverse Fourier transforms, which is possible thanks to its structured mesh. On the other hand, both YALES2 and Nek5000 use an iterative solver with an efficient preconditioning technique; this method is more versatile and should be computationally more efficient for large and complex geometrical configurations. Nek5000 employs CVODE to integrate the thermochemical equations without further splitting of the different terms accounting for convection, diffusion and chemistry. The main differences between the three codes are summarized in Table~\ref{Tab:codes}. Please note that the presented values are those used for the benchmark, even though some other options are available in each codes.

1. ↑ V. Moureau. Yales2 public website www.coria-cfd.fr/index.php/YALES2.
2. ↑ V. Moureau, P. Domingo, and L. Vervisch. From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: Filtered laminar flame-pdf modeling. Combust. Flame, 158(7):1340–1357, 2011.
3. ↑ P. Benard, G. Balarac, V. Moureau, C. Dobrzynski, G. Lartigue, and Y. D'Angelo. Mesh adaptation for large eddy simulations in complex geometries. Int. J. Numer Methods Fluids, 81:719-740, 2015.
4. ↑ A.J. Chorin. Numerical solution of the Navier-Stokes equations. Math. Computation, 22(104):745–762, 1968.
5. ↑ M. Kraushaar. Application of the compressible and low-Mach number approaches to Large-Eddy Simulation of turbulent flows in aero-engines. PhD thesis, PhD, Institut National Polytechnique de Toulouse-INPT, 2011.
6. ↑ M. Malandain, N. Maheu, and V. Moureau. Optimization of the deflated conjugate gradient algorithm for the solving of elliptic equations on massively parallel machines. J. Comput. Phys., 238:32–47, 2013.
7. ↑ C.D. Pierce and P. Moin. Progress-variable approach for large-eddy simulation of non-premixed turbulent combustion. J. Fluid Mech., 504:73–97, 2004.
8. ↑ .G. Goodwin. An open-source, extensible software suite for CVD process simulation. Chemical Vapor Deposition XVI and EUROCVD, 14:2003–2008, 2003.
9. ↑ J. Hirschfelder, R.B. Bird, and C.F. Curtiss. Molecular Theory of Gases and Liquids. Wiley, 1964.
10. ↑ W. Sutherland. LII. the viscosity of gases and molecular force. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 36(223):507–531, 1893.
11. ↑ ^{11.0 11.1} S.D. Cohen, A.C. Hindmarsh, and P.F. Dubois. CVODE, a stiff/nonstiff ODE solver in C. Computers Phys., 10(2):138–143, 1996
12. ↑ ^{12.0 12.1} A.C. Hindmarsh and R. Serban. User documentation for CVODE. [https://computing.llnl.gov/sites/ default/files/public/cv_guide.pdf].
13. ↑ ^{13.0 13.1} A.C. Hindmarsh, P.N. Brown, K.E. Grant, S.L. Lee, R. Serban, D.E. Shumaker, and C.S. Woodward. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans. Math. Soft. (TOMS), 31(3):363–396, 2005.
14. ↑ N. Li and S. Laizet. 2DECOMP&FFT - a highly scalable 2D decomposition library and FFT interface. In Cray User’s Group 2010 conference, Edinburgh, 2010.
15. ↑ Create analytical jacobian matrix source code for chemical kinetics.
16. ↑ Cantera.
17. ↑ A. Ern and V. Giovangigli. Fast and accurate multicomponent transport property evaluation. J. Comput. Phys., 120:105–116, 1995.
18. ↑ A. Abdelsamie and D. Th´evenin. On the behavior of spray combustion in a turbulent spatially-evolving jet investigated by direct numerical simulation. Proc. Combust. Inst., 37(3):2493–2502, 2019.
19. ↑ A. Abdelsamie and D. Th´evenin. Nanoparticle behavior and formation in turbulent spray flames investigated by DNS. In M. Garc´ıa-Villalba, H. Kuerten, and M. Salvetti, editors, Direct and Large Eddy Simulation XII, volume 27 of ERCOFTAC Series. Springer, 2020.
20. ↑ C. Chi, A. Abdelsamie, and D. Th´evenin. A directional ghost-cell immersed boundary method for incompressible flows. J. Comput. Phys., 404:109122–109142, 2020.
21. ↑ A. Abdelsamie, G. Fru, F. Dietzsch, G. Janiga, and D. Th´evenin. Towards direct numerical simulations of lowMach number turbulent reacting and two-phase flows using immersed boundaries. Comput. Fluids, 131(5):123– 141, 2016.
22. ↑ C. Chi, A. Abdelsamie, and D. Th´evenin. Direct numerical simulations of hotspot-induced ignition in homogeneous hydrogen-air pre-mixtures and ignition spot tracking. Flow Turbul. Combust., 101(1):103–121, 2018.
23. ↑ T. Oster, A. Abdelsamie, M. Motejat, T. Gerrits, C. R¨ossl, D. Th´evenin, and H. Theisel. On-the-fly tracking of flame surfaces for the visual analysis of combustion processes. Comput. Graph. Forum, 37(6):358–369, 2018.
24. ↑ A. Abdelsamie and D. Th´evenin. Impact of scalar dissipation rate on turbulent spray combustion investigated by DNS. In M. Salvetti, V. Armenio, J. Fr¨ohlich, B. Geurts, and H. Kuerten, editors, Direct and Large-Eddy Simulation XI, volume 25 of ERCOFTAC Series. Springer, 2019.
25. ↑ Nek5000 version v17.0, Argonne National Laboratory, IL, U.S.A.
26. ↑ ^{26.0 26.1} A.G. Tomboulides, J.C.Y. Lee, and S.A. Orszag. Numerical simulation of low Mach number reactive flows. J. Sci. Comput., 12:139–167, 1997.
27. ↑ ^{27.0 27.1} M.O. Deville, P.F. Fischer, and E.H. Mund. High-order Methods for Incompressible Fluid Flows. Cambridge University Press, 2002.
28. ↑ A. Tanarro, F. Mallor, N. Offermans, A. Peplinski, R. Vinuesa, and P. Schlatter. Enabling adaptive mesh refinement for spectral-element simulations of turbulence around wing sections. Flow Turbul. Combust., 105(2):415– 436, 2020.
29. ↑ R.J. Kee, F.M. Rupley, and J.A. Miller. Chemkin-II: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics. SAND-89-8009, 1989.
30. ↑ M. Schmitt, C.E. Frouzakis, A.G. Tomboulides, Y.M. Wright, and K. Boulouchos. Direct numerical simulation of the effect of compression on the flow, temperature and composition under engine-like conditions. Proc. Combust. Inst., 35:3069–3077, 2015.
31. ↑ G.K. Giannakopoulos, C.E. Frouzakis, P.F. Fischer, A.G. Tomboulides, and K. Boulouchos. LES of the gasexchange process inside an internal combustion engine using a high-order method. Flow Turbul. Combust., 104:673–692, 2020.
32. ↑ B.O. Arani, C.E. Frouzakis, J. Mantzaras, and K. Boulouchos. Three-dimensional direct numerical simulations of turbulent fuel-lean H2/air hetero-/homogeneous combustion over Pt with detailed chemistry. Proc. Combust. Inst., 36(3):4355–4363, 2017.