On the implementation of flux limiters in algebraic frameworks

Computer Physics Communications 271 (2022) 108230 Contents lists available at ScienceDirect Computer Physics Communications www.elsevier.com/locate/cpc On the implementation of flux limiters in algebraic frameworks✩ Nicolás Vallea, Xavier Álvarez-Farréa, Andrey Gorobetsb, Jesús Castroa,∗, Assensi Olivaa, F. Xavier Triasa aHeat and Mass Transfer Technological Center, Technical University of Catalonia, Carrer Colom 11, 08222 Terrassa (Barcelona), Spain bKeldysh Institute of Applied Mathematics, Russian Academy of Sciences, Miusskaya Sq. 4, 125047 Moscow, Russia a r t i c l e i n f o a b s t r a c t Article history: Received 1 June 2021 Received in revised form 5 October 2021 Accepted 9 November 2021 Available online 16 November 2021 Keywords: Flux limiter Parallel CFD Heterogeneous computing Portability Mimetic The use of flux limiters is widespread within the scientific computing community to capture shock discontinuities and are of paramount importance for the temporal integration of high-speed aerodynamics, multiphase flows and hyperbolic equations in general. Meanwhile, the breakthrough of new computing architectures and the hybridization of supercomputer systems pose a huge portability challenge, particularly for legacy codes, since the computing subroutines that form the algorithms, the so-called kernels, must be adapted to various complex parallel programming paradigms. From this perspective, the development of innovative implementations relying on a minimalist set of kernels simplifies the deployment of scientific computing software on state-of-the-art supercomputers, while it requires the reformulation of algorithms, such as the aforementioned flux limiters. Equipped with basic algebraic topology and graph theory underlying the classical mesh concept, a new flux limiter formulation is presented based on the adoption of algebraic data structures and kernels. As a result, traditional flux limiters are cast into a stream of only two types of computing kernels: sparse matrix-vector multiplication and generalized pointwise binary operators. The newly proposed formulation eases the deployment of such a numerical technique in massively parallel, potentially hybrid, computing systems and is demonstrated for a canonical advection problem. ©2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction The evolution in hardware technologies enables scientific computing to advance incessantly and reach further aims. Nowadays, the use of high-performance computing (HPC) systems is rather common in the solution of both industrial and academic scale problems. However, many algorithms employed in scientific computing have a very low arithmetic intensity (AI), which is the ratio of computing work in floating-point operations (FLOP) to memory traffic in bytes, hence numerical simulation codes are usually memory-bounded, making processors suffer from serious data starvation [1,2]. To top it off, the calculations often result in irregular, non-coalescing memory access patterns, reducing the memory access efficiency. Ironically, the memory bandwidth of computing hardware grows much slower than its peak performance, aggravating the problem. All of this motivates the introduction of new parallel architectures with faster and more complex memory configurations into HPC systems. ✩The review of this paper was arranged by Prof. N.S. Scott. *Corresponding author. E-mail address: jesus.castr[email protected] (J. Castro). To take advantage of the increasing variety of computing architectures and the hybridization of HPC systems, the computing subroutines that form the algorithms, the so-called kernels, must be adapted to complex paradigms such as distributed-memory (DM) and shared-memory (SM) multiple instruction, multiple data (MIMD) parallelism, and stream processing (SP). It also encourages the demand for portable and sustainable implementations of scientific simulation codes [3]. While portability is an intangible characteristic of software, it may be easy for a developer to have an idea of how difficult it is to rewrite, debug and verify a specific code on its adaptation to a new architecture. On the other hand, sustainability refers to developing reusable and resilient codes. The way a code is conceived at its inception enormously determines the degree to which both properties can be attained. Traditionally, the development of scientific computing software is based on calculations in iterative stencil loops (ISL) over a discretized geometry—the mesh. This implementation approach is referred to as stencil-based. Despite being intuitive and versatile, the interdependency between algorithms and their computational implementations in stencil applications usually results in a large number of subroutines and introduces an inevitable complexity when it comes to portability and sustainability [4]. https://doi.org/10.1016/j.cpc.2021.108230 0010-4655/©2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Regarding portability, the complexity of stencil applications motivates the adoption of conservative strategies, which consist of porting (rewriting) the most time-consuming part of an existing code, or even the entire code, to a new architecture but minimizing the structural modifications. In other words, it leads to a partial or complete reimplementation of an existing code. These strategies were common during the rise of general-purpose computing on graphics processing units (GPUs) because they allow for direct comparison studies of both numerical and performance results versus the legacy versions. Well-known commercial computational fluid dynamics (CFD) codes and open-source platforms offer GPU extensions for solvers of systems of linear algebra equations (SLAE), which represent a significant part of the overall computing time. This can provide substantial acceleration with compactly localized changes in the code. Such an example can be found in [5], where the authors coupled a GPU-accelerated library for solving large sparse SLAE with the OpenFOAM platform and demonstrated performance on up to 128 nodes of a GPU-based cluster. The use of GPUs in scientific computing is nowadays rather mature, and there are many successful examples in the literature [6–9]. For instance, the early GPU implementations in [10], extended in [11], proved to be two orders of magnitude faster than its central processing unit (CPU) counterpart. Moreover, the solution of two-phase flows on multi-GPU systems [12]was not only faster but also more energyefficient. An example of a GPU porting of an open-source Navier– Stokes solver, the AFiD code, is found in [13]. Further examples of multi-GPU simulations of supersonic and hypersonic flows can be found in [14]. One of the most impressive GPU-based simulations is found in [15], after [16], on the solution of turbulent flows, reporting a sustained performance of 13.7PFLOPS. Regarding sustainability, the implementation of new physics or numerical methods in a stencil-based framework, or its specialization for different mesh types, usually requires the design of new computing subroutines and data structures. This represents the main drawback of such an approach because the effort is not necessarily accumulative and thus reduces the software’s sustainability. To address this, some authors propose domain-specific tools to generalize the stencil computations for specific fields. For instance, a framework that automatically translates stencil functions written in C++ to both CPU and GPU codes is proposed in [17]. However, these generalizations are still heavily restricted by the shape of the stencil they target. An alternative to stencil implementations is to break the aforementioned interdependency between algorithm and implementation so that the calculations are cast into a minimalist set of kernels. In other words, the idea is to use the classical ISL just for building data and leave the calculations to a reduced set of basic operations; thus, legacy codes may be maintained indefinitely as preprocessing tools, and the calculation engines become easy to port and optimize. By casting discrete operators and mesh functions into sparse matrices and vectors, it has been shown that nearly 90% of the calculations in a typical CFD algorithm for the direct numerical simulation (DNS) and large eddy simulation (LES) of incompressible turbulent flows boil down to the following basic linear algebra subroutines: sparse matrix-vector product (SpMV), linear combination of vectors (axpy) and dot product (dot) [18]. Moreover, after the generalizations detailed in Section 3.2 this value will be raised to 100%. Hereinafter, we refer to this implementation based on algebraic subroutines as algebraic or algebra-based. In this algebraic approach, the kernel code shrinks to dozens of lines; the portability becomes natural, and maintaining OpenMP, OpenCL, or CUDA implementations takes minor effort. Besides, standard libraries optimized for particular architectures (e.g.,cuSPARSE[19], clSPARSE [20]) can be easily linked in addition to specialized inhouse implementations. A similar approach is found in PyFR [16], where the majority of operations are cast in terms of matrixmatrix multiplications linking with appropriate BLAS libraries. In the context of the DNS, the preconditioned conjugate gradient (CG) method following such an algebraic approach was implemented in [21], and its potential was exploited in [22]to perform petascale CFD simulations. In our previous work, we proposed a heterogeneous implementation of an algebraic framework [23], the hierarchical parallel code for high-performance computing (HPC2), as a portable solution with many potential applications in the fields of computational physics and mathematics. Later, the minimalist design allowed us to easily optimize our framework for HPC systems with non-uniform memory access (NUMA) configurations [24], just by changing the definition of the three kernels and the initialization of vectors and matrices. While linear schemes fit naturally in the proposed algebraic approach, the implementation of non-linear schemes is not evident. However, as a consequence of Godunov’s theorem [25,26], those are required for the treatment of shock discontinuities to avoid unstable discretizations and the onset of wiggles. Such discontinuities are present in many industrial applications and appear in both compressible and multiphase flows, representing a classical problem in numerical analysis since, at least, the 1960’s. The construction of stable, second-order (and higher) discretizations, then, requires the adoption of non-linear schemes to exhibit a total variation diminishing (TVD) [27]behavior. Among them, flux limiters are a mature and robust method, which has been adopted in a diversity of applications. Sweby [28] generalized several limiters and stated the conditions for second-order TVD schemes in a onedimensional homogeneous mesh in its well-known Sweby diagram. Despite the known inconsistencies that arise when departing from the one-dimensional homogeneous case [29,30]these techniques have been ported to non-homogeneous Cartesian [30] and unstructured [31]meshes as well. Advances in this field have also been exploited by the multiphase flow community, particularly for the advection of the marker function [32,33]. Both the analysis and the implementation of flux limiters are typically performed from the aforementioned stencil-based perspective. However, the growing interest of the community in mimetic methods [34]unveils an alternative to the implementation of flux limiters. Mimetic methods construct discrete operators directly from the inherent incidence matrices that define the mesh. Adopting such an approach presents an important advantage from both theoretical and practical points of view. On the one hand, this allows for a flawless discrete mimicking of the continuum operators, facilitating the exact conservation of important secondary properties, such as energy [35,36], among others. On the other hand, this reduces the implementation to the right combination of a reduced set of algebraic subroutines. Therefore, the present work is devoted to the formulation of flux limiter schemes and their implementation into algebraic frameworks. The proposed implementation will be analyzed from a computational point of view and compared with the classical stencil counterpart, and the benefits of each option will be discussed in depth. The rest of the paper is organized as follows. In section 2a review of basic concepts of graph theory is briefly summarized in order to provide some context. Section 3develops a generalization of flux limiters from an algebraic perspective and introduces the matrix-based calculation of the gradient ratio. Section 4 highlights the capabilities of the method on a well-known threedimensional deformation benchmark. Finally, the discussion and conclusions are stated in Sections 5and 6. 2. Algebraic topology By using concepts from algebraic topology, mimetic methods preserve the inherent structure of the space, leading to stable and 2 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Fig. 1. 2D mesh composed of cells ci, which are bounded by faces fjoriented in the direction ˆ nj. Faces, which collapse to the same entities as edges, are bounded by the set of vertices vk. robust discretizations [37,34]. However, the development of such techniques is out of the scope of this paper, where we rather focus on exploiting the relationships between the different entities of the mesh for the construction of flux limiters. The interested reader is referred to [38] and references therein. Given whatever space of interest , we can equip it with a partition of unity, namely a mesh M, by bounding the group of cells, C, with faces, F; those with the set of edges, E, and finally those with the set of vertices, V. In this sequence, groups are related to the next element of the sequence by means of the boundary operator ∂. This is known as a chain complex [37,39]. A two-dimensional example can be seen in Fig. 1, where faces and edges collapse into the same entity. In every space Q, discrete variables are arranged in arrays such as θq∈R|Q|. The relationship between the bounding elements of a geometric entity can be cast in oriented incidence matrices, EE→V, EF→Eand EC→F, corresponding to edge-to-vertex, face-to-edge and cell-to-face, respectively. The corresponding incidence matrices for the mesh depicted in Fig. 1read: EF→V= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0+1+1+10 0 0−10 000 −1+10000 −10−10−1−10 0 0 00000 +1−10 0 +1−10 0 0 0−10+1 0000000 +1−1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,(1) EC→F= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0−1+10 0−10+1 −1+10 0 +10 0 0 +10 0 0 00 −10 00 +10 000 +1 000 +1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(2) Conversely, we can define EF→C, EE→Fand EV→E, for the face- to-cell, edge-to-face and vertex-to-edge incidence matrix. Such converse incidence matrices are obtained by transposition: EF→C=EC→FT, EE→F=EF→ET, EV→E=EE→VT. Incidence matrices represent the boundary operator between one element of the chain and the next one. Following the example of Fig. 1, EF→Cand EV→Fprovide with the orientation of the boundary faces fjfor every cell ciand the boundary vertices vk for every face fj. Incidence matrices also play an essential role in preserving properties of the discrete space. In particular, they form an exact sequence. Exact sequences are those such that the application of the boundary operator twice results in 0. This can be verified by checking EF→VEC→F=0. This property is shared by its continuum counterpart, the de Rahm cohomology [34], which is the ultimate responsible of the following vector calculus identities [37]: ∇×∇≡0,(3) ∇·∇×≡0.(4) These are powerful identities that mimetic methods preserve by construction. For an extended review of the relationship between the continuum and the discrete counterparts, the reader is referred to [37,34] and references therein. In addition to provide a suitable platform for the construction of appropriate mimetic methods, the relations contained in incidence matrices can be studied from a graph theory perspective. A straightforward use of incidence matrices allows to compute differences across faces. The fact that differences lie in a different space (faces) than variables (cells) is an inherent property of such an approach: θc=EC→Fθc.(5) Particularly useful is the construction of undirected incidence matrices (BQ→S), which are built by taking the absolute value of the elements of the directed ones (EQ→S). Considering the index notation between a generic space Q(e.g., cells, faces) and its boundary S(e.g.,faces, edges), we could proceed as follows: [BQ→S]sq =bsq =|esq|,(6) where esq =[EQ→S]sq. Similarly, one can proceed to compute the degree matrix of the graph, which accounts for the number of connections that an entity has (e.g.,the number of cells in contact with a face). Degree matrices are always diagonal and the value of the diagonal elements is obtained as follows: WQQ =diag(BS→Q1S). (7) In particular, undirected incidence matrices can be used to construct suitable shift operators [40]: C→F=WFF−1BC→F.(8) This provides with a simple face-centered interpolation, weighted with the number of adjacent faces. Note that by taking this approach, boundaries are inherently included from the graph information. The use of such C→Fis restricted to scalar fields. However, following [40], this can be readily extended to vector fields as follows. First, the discretization of a continuum vector field,  u= (u1, ..., ud)T, is arranged in a single column vector, uc∈Rd|C|×1= (u1, ..., ud)T, where ui=((ui)1, (ui)2, ..., (ui)|C|)Tare the vectors containing the components corresponding to the ith spatial direction. Note that dis the number of spatial dimensions of the problem. Next, the interpolator can be extended component-wise by applying the Kronecker product with the identity matrix of size Rd×d. The final ensemble is as follows: C→F=Id⊗C→F.(9) 3 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Fig. 2. Classical stencil for the computation of the gradient ratio at face fj. U, Cand Dcorrespond to the upstream, centered and downstream nodes. Similarly, normal vectors can be cast into a R|F|×d|F|matrix by arranging ddiagonal matrices, corresponding to every component of the face vector, next to each other as NF=(N1| ...|Nd)[40]. The two-dimensional NFmatrix corresponding to the mesh depicted in Fig. 1reads: NF= ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ n1x0... 0n1y0... 0 0n2x... 00n2y... . . . . . .0.... . .. . .0.... . . 00... n9x00... n9y ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ .(10) In such a way, it is straightforward to either project a discrete vector as NFxf, or to vectorize a scalar quantity as NFTsf, provided that both are stored at the faces. An accurate discussion about the construction of this matrix can be found in [40]. Other basic matrices derived from the graph are the graph Laplacian (LCC) and the adjacency matrix (ACC ): LCC =BF→CBC→F,(11) ACC =WCC −LCC.(12) Both are constructed based on the incidence matrices and provide information about the propagation of information along the graph. They are constructed by connecting cells to its neighbors through its bounding faces. In summary, the constructor of such operators provides with tools able to relate different elements of the graph between each others. Equipped with such basic concepts, the development of higher level operators can proceed as in the following section. 3. Flux limiters The solution of hyperbolic problems in finite volume methods when sharp discontinuities are present requires the use of high resolution schemes in order to attain second order approximations. In turn, the construction of such schemes is reduced to the appropriate reconstruction of the flux at the faces. Prone to introduce numerical instabilities, such a reconstruction requires an appropriate flux reconstruction strategy in order to guarantee TVD behavior (i.e.,such that no new minima or maxima are introduced). This is attained by limiting the flux at cell’s boundaries by means of a flux limiter function. Typically, flux limited schemes are stated in the following form [28]: θf=θC+(r)θD−θC 2,(13) where θCand θDstand for the centered and downwind values of θ according to the velocity field uand (r)is the flux limiter function. Fig. 2depicts this situation. From a physical point of view, this is equivalent to the introduction of some sort of artificial diffusion which stabilizes the method at the expense of smearing out its profile. This can be easily seen by rewriting the classical stencil-based formulation stated in equation (13)into: θf=θC+θD 2+(r)−1 2(θD−θC),(14) where ((r) −1)/2 stands for the artificial diffusion added to a classical symmetry-preserving scheme. The limiting approach has been used by several authors [41,26] who over the years developed several discontinuity sensors in order to limit the dissipation to the region near the shock. Among all discontinuity sensors, the most popular is the use of the gradient ratio. Following the nomenclature in Fig. 2, this is defined as follows [28]: rf=Uθ uθ=θC−θU θD−θC,(15) where Uθis the gradient of θat the upwind face while uθcorrespond with the gradient at the face of interest. Both differences are taken as positive in the flow direction, defined by the sign of the velocity field uf. This provides with an intuitive description of where discontinuities are and its order of magnitude, at the time that keeps a compact stencil. In addition, it allows to, after proper manipulation by the flux limiter itself, limit the flux in a way which can be interpreted as a diffusion-like term. This is known as “upwinding” as it has the same effect as recovering the 1st order upwind discretization near shocks. TVD conditions in terms of the gradient ratio were stated by Harten in 1983 [27]for one-dimensional homogeneous grids. These conditions were used by Sweby in 1984 [28]to state 2nd order and TVD conditions for different forms of flux limiters. This idea has been extended, with several degrees of accuracy, to multidimensional and irregular grids [29,30], among others. 3.1. Algebraic formulation As stated previously, the discretization of equation (14)may benefit from the adoption of an algebraic approach. In this regard, it can be easily extended to the whole computational domain as: θf=(C→F+F(rf)Q(uf)EC→F)θc,(16) where θf∈R|F|and θc∈R|C|are the vectors holding all the values of θfand θc, C→F∈R|C|×|F|is the cell-to-face interpolation defined in equation (8), F(rf)∈R|C|×|F|is the diagonal matrix absorbing the artificial diffusion introduced in equation (14), Q(uf)∈R|F|×|F|is the diagonal matrix taking the proper sign of the velocity at the faces and EC→Ftakes the difference across them as in equation (5). At this point, we may be tempted to analyze the construction of new flux limiters by means of basic algebra concepts. In particular, to bound its spectrum by means of Gershgorin’s theorem or to check its entropy conditions [42], among others. The interested reader is referred to Báez et al. [43], where a similar approach is taken for spatial filters. While both C→Fand EC→Fare readily available from the background stated in section 2, the construction of F(rf)by means of basic algebraic operations solely is addressed next. Because flux limiter functions F(rf)depend only on the local value of rfand we defer the details on the implementation of the pointwise operations to section 3.2, the problem is turned into the accurate computation of rfat faces. There has been several approaches [31,29]to the construction of rfin terms of a leastsquares reconstructed gradient. However, the implementation of such schemes can be cumbersome and may not, eventually, recover the one-dimensional homogeneous solution when a homogeneous structured mesh is used. 4 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 The construction of the gradient ratio will proceed first by the separate calculation of both the numerator (Uθ) and the denominator (uθ) of equation (15), then computed as: [rf]i=[dUθ]i [duθ]i,(17) where duθ∈R|F|is the face-centered vector holding the difference across the face taken in the direction stated by uf, while dUθ∈ R|F|holds the upstream differences according, again, to uf. In this approach we propose to employ symmetry-preserving gradients (see [35]) into the calculation of both face-centered and upstream gradients in order to preserve, as much as possible, the mimetic properties of the approach. In addition, we aim at recovering the Cartesian formulation as in [33]. Before any calculation, the sign matrix (Q(uf)) is constructed by assigning to a R|F|×|F|diagonal matrix +1 for a positive velocity and -1 for a negative one: [diag (Q(uf))]i=sign([uf]i). (18) This allows a straightforward calculation of the velocity-oriented gradient at the face as follows: duθ=Q(uf)EC→Fθc,(19) where Q(uf)is used to provide the right direction in which the difference is taken according to the velocity field. The construction of the upwind difference dUθis more involved. The idea is to construct a partial adjacency matrix which only considers upstream faces, namely the upstream adjacency matrix, AU FF(u), which is responsible to garner upstream information and will be defined further in this paper. We proceed as follows: EC→Fis used to compute the difference across every face according to equation (5). In order to assess the contribution of every neighboring face to the face of interest, differences are vectorized with its corresponding normal using NFT and added all together with AU FF(u). Finally, the resulting value is projected over the normal of the face of interest by means of NF. The overall construction of the operator is then: dUθ=NFId⊗AU FF(u)NFTEC→Fθc,(20) where, similarly to equation (9), we reuse the AU FF(u)operator for all spatial dimensions. Face normal matrices NF, defined as in equation (10), are used for both vectorization and projection of the neighboring differences. In this way, orthogonal meshes recover the original one-dimensional formulation, whereas unstructured ones are handled inherently by the incidence matrix. We are now left with construction of the upstream adjacency matrix, AU FF(u), which may look cumbersome at a first glance. It can be assembled from other simpler matrices: AU FF(u)=1 2AFF −Q(uf)AD FF,(21) where AFF is the face adjacency matrix, AD FF is a “directed adjacency matrix”, which will be introduced below, and Q(uf)is the already defined velocity sign matrix. The strategy for the construction of AU FF(u)is to add the contribution of all neighboring faces irrespective of the flow direction and then use Q(uf)AD FF, to remove the downwind ones depending on the values of Q(uf). Note that both AFF and AD FF are constant matrices and that the only matrix that needs to be updated according to ufis Q(uf). The construction of AFF proceeds similarly to equation (12)as follows: AFF =WFF −BC→FBF→C.(22) Fig. 3. Upstream (red) and downstream (blue) adjacent faces for face f1with respect to the positive component of velocity at face f1. The selection of the right ones in AU FF(u)will ultimately depend on Q(uf). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.) As seen in section 2, the adjacency matrix is symmetric and contains non-negative entries only. Following on the example depicted in Fig. 3, the corresponding AFF reads: AFF = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0+1+10 0+1+10 0 +10+10000+1+1 +1+10+1+10 0 0 0 00 +10+10 0 0 0 00 +1+100000 +100000+10 0 +10 0 0 0+10 0 0 0+1000000+1 0+100000+10 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(23) On the other hand, the construction of AD FF allows us to distinguish neighboring faces which lie, according to the face normal, behind or ahead of the face in question. This requires the inclusion of the directed incidence matrix EC→Finto the calculation of the adjacency matrix as: AD FF =EC→FBF→C,(24) which provides with the following matrix: AD FF = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0+1+10 0−1−10 0 +10+10000−1−1 −1−10+1+10 0 0 0 00 −10−10 0 0 0 00 −1−100000 +100000+10 0 −10 0 0 0−10 0 0 0−1000000−1 0−100000−10 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .(25) Note that AD FF has the same pattern as AFF but entries corresponding to faces located upstream (with respect to the face normal direction) contain −1 whereas those located downstream contain +1, as shown in Fig. 3. However, the choice of upstream/downstream faces should depend on the faces local velocity and not on its arbitrary choice of face normal. The product Q(uf)AD FF, corrects this by inverting the sign of the rows corresponding to the faces whose velocity component is not aligned with the face normal. The result is a correct choice of upstream and downstream faces according to the local face velocity. Finally, the combination of AD FF and AFF in equation (20)results in: 5 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Table 1 Particularizations of the operator ◦and pointwise function f(xi)to represent various kernels using the generalized binary operator described in Equation (28). The superbeexty corresponds to the SUPERBEE flux limiter [28]. ◦f(xi)AI Resulting Kernel +axi1/12 axpy ×axi1/12 axty /axi1/12 axdy ×xi>0? +1:−11/12 signxty ×0.5(max(0,max(min(1,2xi),min(xi,2))) −1)7/24 superbeexty dUθ=NFId⊗1 2Q(uf)AD FF −AFFNFTEC→Fθc.(26) While equation (26) succeeds at selecting the proper upstream faces, its direct implementation involves many redundant operations that may result in an unnecessary overhead. For this reason, the computation of dUθis rearranged as follows: dUθ=(Q(uf)SC→F+TC→F)θc,(27) where we introduce the new matrices SC→F=1 2NF(Id⊗ AD FFNFTEC→Fand TC→F=1 2NF(Id⊗AFF)NFTEC→F, which can be precomputed at the beginning of the simulation. 3.2. Algebraic implementation In previous works of Oyarzun et al. [18] and Álvarez et al. [23], an algebra-based implementation model was proposed for the DNS and LES of incompressible turbulent flows such that the algorithm of the time-integration phase reduces to a set of only three algebraic kernels: SpMV,axpy and dot. However, a close look at Equations (17) and (18), for instance, reveals that this set is insufficient to fulfill the implementation of the flux limiter because it comprises non-linear operations. Nevertheless, instead of being an inconvenience, this encourages us to demonstrate the high potential of our algebraic strategy again. We propose the generalization of the axpy via the introduction of a generalized binary operator (kbin) that performs any given pointwise arithmetic calculation such that: yi←yi◦f(xi). (28) This binary operator can easily map to the required kernels by defining ◦and f(xi)as outlined in Table 1. Similarly, the dot kernel can be turned into a generalized reduction operator (kred), r←r◦f(xi), (29) which can easily represent any required reduction operation such as the calculation of the norm of a vector or the Courant- Friedrichs-Lewy (CFL) condition. However, the details of the kred are out of the scope of this work because it is not required for the implementation of the flux limiter. From a computational point of view, this kernel generalization does not alter the implementation of the original axpy: it still performs simple, pointwise arithmetic operations over the vector elements and provides uniform, aligned and coalescing memory accesses which suits the simple instruction, multiple data (SIMD) and SP paradigms perfectly. Therefore, having already efficient implementations of axpy for different architectures, the implementation of kbin is straightforward (e.g., consider the use of function pointers, templates, macros, among others). On the other hand, the arithmetic intensity of this new kernel is not a fixed value anymore, as shown in Table 1. While the AI of the axpy was 1/12 FLOP per byte (one product and one addition per three 8-byte values), that of the kbin will depend on the specific arithmetic calculations involved in the function f(xi). This allows us to significantly increase the AI in our calls by means of kernel fusion, reduce the number of intermediate results, and thus reduce the time-to-solution. The final algorithm for the deployment of a flux limiter in the reconstruction of the variable at faces, θf, within our algebra-based framework is described in Algorithm 1. Note that because we are Algorithm 1 Algorithm for reconstruction of a scalar field at faces, θf, using the algebraic implementation of a flux limiter. Require: θc, uf, EC→F, SC→F, TC→F, C→F Ensure: θf 1: duθ←EC→FθcSpMV 2: duθ←Q(uf)duθsignxty 3: rf←SC→FθcSpMV 4: rf←Q(uf)rfsignxty 5: rf←rf+TC→FθcSpMV 6: rf←rf/duθaxdy 7: θf←F(rf)duθsuperbeexty 8: θf←θf+C→FθcSpMV actually interested in the evaluation of dUθrather than in the construction of the operator itself, matrix-matrix products are avoided and successive matrix-vector products are performed. Similarly, we are not interested in the construction of the diagonal matrices Q(uf)and F(rf). We will rely on the generalized binary operator instead. Therefore, the evaluation of Q(uf)ywill be done using the signxty in Table 1as follows: [y]i←[y]i×([uf]i>0? +1:−1). (30) Finally, a particular flux limiter function must be chosen on the evaluation of F(rf)y. Our framework allows to easily switch between different functions. Hereinafter, we will proceed with the superbee [28], represented by superbeexty: [y]i←[y]i×0.5(max(0,max(min(1,2[rf]i),min([rf]i,2))) −1). (31) In conclusion, the evaluation of a scalar field at faces, θf, as in Algorithm 1can be fitted in our algebra-based framework by combining two types of computing kernels: SpMV and kbin. 3.3. Comparison with stencil-based implementations Our algebraic implementation for the appropriate reconstruction of the variables at faces, θf, is now compared with classical, stencil-based approaches. This comparison is conducted firstly from a theoretical point of view, assessing the minimum number of FLOP and memory traffic (in bytes) required in different scenarios. Note that the actual number of memory accesses during kernel execution depends not only on the algorithm but also on hardware and software features. Therefore, regarding the memory traffic, two different values are estimated considering the full-hit and full-miss caseloads. The former refers to the best scenario with an ideal temporal locality: multiple accesses to a particular data element are so close in time that its value is always reused from cache. Conversely, the latter considers the worst scenario with a null temporal locality so that every repeated access results in cache-miss and requires a memory load from memory. Thus, these two values result in the interval of effective AI of each kernel. For the sake of clarity, in this comparison we only consider k-regular, periodic meshes composed of convex polygons which accomplish the following equality: |F|≈k 2|C|,(32) 6 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Fig. 4. Example of stencils in a one-dimensional Cartesian grid according to the classical flux limiter approach. The stencil topology varies according to the sign of the velocity field,  v. where |C|and |F|represent the number of cells and faces, respectively, and kis the degree of the mesh elements (i.e.,the number of neighboring cells). The resulting requirements are listed in Table 2, which is described throughout the section. Let us start from the analysis of the simplest case: the stencilbased calculation of θfin a one-dimensional Cartesian grid, depicted in Fig. 4and described in Algorithm 2. Every ith face is surrounded by two cells, ciand ci+1. For each ith face in |F|, the sign of the velocity determines the upstream, centered and downstream values of θ, as well as the centered and upstream distances. Following the algorithm, 5|F|FLOP are required for computing the gradient ratio in line 9, 1|F|for computing the limiter function in line 10 (this value may vary depending on the limiter function; in this example we have considered the superbee limiter [28]), and 7|F|for computing the value at the face in line 11. In the algorithm, three discrete fields are required for the computations: the initial scalar field, θc∈R|C|, the velocity field, uf∈R|F|and the distances between cells, df∈R|F|. Besides, the algorithm ensures the calculation of the discrete scalar field at faces, θf∈R|F|. Thus, considering double-precision values, and given that the number of faces is equal to the number of cells in the one-dimensional case (k =2 →|F| =|C|), the minimum FLOP and bytes required are 5|F| +1|F| +7|F| =13|C|and 8(|C| +|F| +|F| +|F|) =32|C|, respectively. The total memory traffic in the full-miss caseload would rise to 48|C|because two different values of dfand θcare accessed for every face. Note that these values are slightly reduced in the particular case of uniform meshes: neither the distances array nor its quotient is required, thus omitting 2|F|FLOP in line 9 and the access to df. On the other hand, the generalization of Algorithm 2for three-dimensional Cartesian grids (k =6 →|F| =3|C|) is straightforward and rises the computational requirements to 39|C|FLOP and 80–336|C|bytes for non-uniform meshes, or 33|C| FLOP and 56–216|C|bytes in the uniform case. Algorithm 2 Stencil-based calculation of θfin a onedimensional Cartesian grid. Require: θc, uf, df8C + 8F + 8F bytes Ensure: θf8F bytes 1: for i ←1to |F|do 2: if uf[i] >0then 3: θC=θc[i],θ U=θc[i −1], θD=θc[i +1] 4: di=df[i],dU=df[i −1] 5: else 6: θC=θc[i +1], θU=θc[i],θ D=θc[i +2] 7: di=df[i],dU=df[i +1] 8: end if 9: ri←(di/dU)(θC−θU)/(θD−θC)5F flops 10: i←limiter(rf)1F flops 11: θf[i] ←(θC+θD)/2 +(i−1)(θD−θC)/27F flops 12: end for A generalization of the stencil calculation for unstructured meshes is outlined in Algorithm 3. In contrast with the structured algorithm, the indices of neighboring nodes are not predictable, so the incidence graphs are required. For each ith face in |F|, the sign of the velocity determines the indices of the centered and downstream cells, cCand cD, according to the cell-to-face incidence graph. Then, for each jth face incident to cC(except fi), its contribution to the upstream gradient, projected over the normal of fi, is accumulated, accounting for 5(k −1)|F|and 4(k −1)|F| FLOP in lines 20 and 21, respectively. The calculation of θfin lines 24 to 26 follows similarly as in Algorithm 2, and adds 13|F|operations. In this case, six discrete fields, one integer list and two incidence graphs are required for the computations. The specific requirements for two k-regular unstructured meshes, k =2(onedimensional mesh) and k =6 (three-dimensional hexaedral mesh), are listed in Table 2. Algorithm 3 Stencil-based calculation of θfin a generic unstructured grid based on incidence graphs. Require: θc, uf, df, nx, ny, nz, 8C + 40F bytes kc, EC→F, EF→C4C + 8F + 8F bytes Ensure: θf8F bytes 1: for i ←1to |F|do 2: if uf[i] >0then 3: cC=EC→F[i][0] 4: cD=EC→F[i][1] 5: else 6: cC=EC→F[i][1] 7: cD=EC→F[i][0] 8: end if 9: U=0 10: for k ←1to kc[cC]do 11: j =EF→C[cC][k] 12: if j = ithen 13: if EC→F[j][0] = cCthen 14: cU=EC→F[j][0] 15: else 16: cU=EC→F[j][1] 17: end if 18: dj=df[j] 19: θU=θc[cU] 20: nk←nx[i]nx[j] +ny[i]ny[j] +nz[i]nz[j]5(k-1)F flops 21: U←U+nk(θC−θU)/dj4(k-1)F flops 22: end if 23: end for 24: ri←U/((θD−θC)/di)3F flops 25: i←superbee(ri)1F flops 26: θf[i] ←(θC+θD)/2 +(i−1)(θD−θC)/27F flops 27: end for Finally, in the algebraic implementation outlined in Algorithm 1, it can be observed how it is completely independent of the mesh type and the numerical method: these characteristics only affect the matrices. The number of calls to SpMV and kbin kernels is readily deduced from the algorithm: 4 times each. Four matrices are required for the computations: the differences at faces, EC→F, with 2|F|non-zero elements, the oriented and unoriented differences, SC→Fand TC→F, with 2k|F|each, and the cell-to-face interpolation, C→F, with 2|F|. In this example, we consider the use of the ELLPACK format [44]in which each nonzero element accounts for 12 bytes (i.e.,8 bytes for the coefficient and 4 bytes for the column index). The specific requirements for two k-regular meshes are listed in Table 2. 4. Numerical study 4.1. Three-dimensional deformation problem Next, the application of this technique is applied to a canonical benchmark. In particular, the deformation (advection) of a sharp 7 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Fig. 5. Contours for θ=0.5inx=0.35, y=0.35 and z=0.35 planes after 3.0 time-units for meshes of 723, 1443, 2883, 4323and 5763cells. Table 2 Minimum number of FLOP and memory traffic (in bytes) required per mesh cell for computing the variable at the faces in different scenarios: stencil- and algebra-based implementations on uniform, non-uniform and unstructured one-dimensional (k = 2) and three-dimensional (k =6) meshes. kApproach FLOP Bytes AI full-hit full-miss full-hit full-miss 2 Uniform 11 32 48 0.344 0.229 2 Non-uniform 13 40 56 0.325 0.232 2 Unstructured 20 84 132 0.238 0.152 2 Algebraic 28 288 352 0.097 0.080 6 Uniform 33 80 240 0.413 0.138 6 Non-uniform 39 104 360 0.375 0.108 6 Unstructured 168 228 1020 0.737 0.165 6 Algebraic 180 1408 1984 0.128 0.091 profile, which has been tested on three-dimensional hexahedral meshes of 723, 1443, 2883, 4323and 5763cells following Algorithm 4, where M ∈R|C|×|F|is the divergence operator [40] and U(uf) ∈R|F|×|F|is a diagonal matrix containing the velocities at faces. Recall we evaluate the products by diagonal matrices by means of kbin calls. Algorithm 4 Algorithm for the advection of a scalar field with a 1st order Euler method, using the algebraic implementation of a flux limiter. Require: θn c, uf, dt, EC→F, SC→F, TC→F, C→F, M Ensure: θn+1 c 1: duθ←Q(uf)EC→Fθn c 2: rf←(Q(uf)SC→F+TC→F)θn c/duθ 3: θf←C→Fθn c+F(rf)duθ 4: θn+1 c←θn c+dtMU(uf)θf The sharp profile is initialized in a physical domain of [0, 1] × [0, 1] ×[0, 1]as a sphere of radius r=0.15, located at (0.35, 0.35, 0.35)and subject to a divergence-free velocity field: u=2sin2(πx)sin(2πy)sin(2πz)cos(πt/T), (33) v=−sin(2πx)sin2(πy)sin(2πz)cos(πt/T), (34) w=−sin(2πx)sin(2πy)sin2(πz)cos(πt/T), (35) during 3.0 time-units, T[45]. The results of the profile on meshes of 723, 1443, 2883, 4323 and 5763cells are shown in Fig. 5for the slices in x =0.35, y = 0.35 and z=0.35 planes. The three-dimensional temporal evolution of the sphere on a mesh of 5763cells is shown in Fig. 6. As in [45], the resulting shapes after the deformation are satisfactory, and mass is exactly conserved. 4.2. Performance analysis The benchmark described in Section 4.1 has been deployed on the HPC2[23]framework, designed for the efficient implementation of algebraic algorithms on hybrid supercomputers. To demonstrate its portability, the simulations have been run on different supercomputers. Before going into details, we define the theoretical achievable performance of a particular kernel, πk, as: πk=min(π,AIkβ), where πis the peak performance of the computing device in double-precision, βis the peak memory bandwidth and AIkis the maximum AI of the kernel, taking the full-hit scenario as described in Section 3.3. Then, we define the performance and memory efficiency as the ratio of measured performance to πkand measured memory traffic to full-hit, respectively. The simulations on meshes of 723–4323cells have been executed on up to 64 nodes (3,072 cores) of the CPU-based MareNostrum 4 supercomputer at the Barcelona Supercomputing Center. Its nodes are equipped with two Intel Xeon 8160 CPUs (24 cores, 2.1 GHz, 6 DDR4-2666 memory channels, 128 GB/s memory bandwidth, 33 MB L3 cache), interconnected through the Intel Omni- Path network (12.5 GB/s). The application achieved a sustained performance of up to 1.6 TFLOPS, corresponding to nearly 80% of performance efficiency. The simulation on a mesh of 5763cells has been executed on 27 nodes of the Lomonosov-2 hybrid supercomputer at Lomonosov Moscow State University. Its hybrid nodes are equipped with one Intel Xeon E5-2697 v3 CPU (14 cores, 2.6 GHz, 4 DDR4-2133 memory channels, 68 GB/s memory bandwidth, 35 MB L3 cache) and one NVIDIA Tesla K40M GPU (12 GB of GDDR5 memory, 288 GB/s, PCIe 3.0 x16 – 16 GB/s), interconnected via InfiniBand FDR network (7 GB/s). The application achieved a sustained performance of 0.9 TFLOPS, corresponding to nearly 75% of performance efficiency and 98.5% of the heterogeneous efficiency (i.e.,the ratio of the heterogeneous performance to the sum of CPU-only and GPU- only performances). Finally, we complete the performance analysis, and the theoretical comparison in Section 3.3, by comparing the actual performance of the algebraic and stencil approaches. Following Algorithm 3, a stencil kernel has also been implemented in HPC2. The tests have been carried out on the CPU-based JFF cluster at the Heat and Mass Transfer Technological Center. Its nodes are equipped with two Intel Xeon Gold 6230 CPUs (20 cores, 2.1 GHz, 6 DDR4-2933 memory channels, 140 GB/s memory bandwidth, 27.5 MB L3 cache). Considering that the DM parallelization is equivalent in both approaches since the data exchanges are the same, only single-node comparisons have been conducted. The results are shown in a roofline plot in Fig. 7for both single- and dual-socket executions using a NUMA-aware, SM parallelization on meshes of 1443(A) and 2883(B) cells. We have discarded 8 N. Valle, X. Álvarez-Farré, A. Gorobets et al. Computer Physics Communications 271 (2022) 108230 Fig. 6. Time evolution of the θ=0.5 contour for t=0, 0.25, 0.75, 1.5, 2.25, 2.75 and 3.0 time-units for mesh of 5763cells. Fig. 7. Roofline models representing the ranges of achievable performance for the stencil- and algebra-based implementations of the flux limiter on three-dimensional unstructured meshes. Results on Intel Xeon Gold 6230 are shown for meshes of 1443(A) and 2883(B). the smallest mesh of 723to ensure a memory-bounded behavior and the biggest meshes of 4323and 5763because they do not fit in a single node. In the plot, two vertical lines represent the minimum and maximum values of the AI for each kernel as estimated in Section 3.3 and outlined in Table 2. Then, the roofline curve is calculated as follows: (π,β)=min(π,AIβ). In this particular test, βand πare 140 GB/s and 1344 GFLOPS per socket, respectively. The real number of memory accesses to main memory have been measured using a profiling tool to calculate the effective AI of each execution. In single-socket execution, the stencil kernel performs nearly twice faster than the algebraic one (stencil: 33.02 (A) and 24.52 (B) GFLOPS, algebraic: 13.89 (A) and 13.54 (B) GFLOPS), even though its lower performance efficiency (stencil: 32% and 24%, algebraic: 78% and 76%). However, this gap is reduced to approximately ×1.35 in the dual-socket case (stencil: 37.01 and 33.82 GFLOPS, algebraic: 27.16 and 25.17 GFLOPS). The algebraic kernels feature a regular unit-strided memory access everywhere except the input vector in SpMV. In contrast, the stencil kernel leads to irregular accesses to EF→C, EC→F, df, θcand  n, resulting in higher cache miss rates and reducing the memory efficiency, especially in dualsocket configurations (stencil: 36% and 33%, algebraic: 96% and 94%). Thus, the actual performance gap is far from the ×5.75 of the (worst) theoretical scenario. For further details of the implementation of our framework and a detailed performance and scalability analysis on different types of supercomputing facilities, the reader is referred to Álvarez-Farré et al. [24]. 5. Discussion The advantages and disadvantages of the algebra- and stencilbased implementations of the flux limiter are discussed below. The algebraic formulation of flux limiters proposed in this work allows for fitting the calculation of such high-resolution schemes into algebra-based frameworks. Following this approach, the DNS and LES of turbulent multiphase flows, for instance, is reduced to a minimal set of algebraic subroutines, leading to fully portable and sustainable implementations. The challenges associated with the introduction of new architectures, and the ongoing hybridization of HPC systems make these two properties of the uttermost importance in the development of modern scientific computing codes. Moreover, algebraic kernels are so widespread that optimized libraries are available for virtually all the existing architectures. From the results listed in Table 2, a significant (theoretical) overhead is revealed in the algebraic implementation. This overhead is mainly induced from Equation (27), where we are enforced to calculate the contribution of all neighboring faces twice as described in Fig. 8, according to the oriented and unoriented matrices, to cancel the downwind values. Note that TC→Fand SC→Fare the matrices with the largest number of non-zero elements. On the other hand, the algebraic implementation requires more kernel calls which reduces the maximum AI of the algorithm. However, the memory access patterns in algebraic kernels are regular and unit-strided everywhere except the input vector in SpMV. 9