Automatic Creation of High-bandwidth Memory Architectures from Domain-specific Languages: The Case of Computational Fluid Dynamics

In numerical mathematics, spectral element methods (SEM) are common in solving partial differential equations, like the Navier–Stokes equations [38], which cannot be solved analytically. spectral element methods (SEM) approximates the solution by transforming the unknown physical quantities of the problem into spectral coefficients. To reduce the numerical complexity, the simulated volume is divided into \(N_{eq}\) smaller elements. To further reduce the error, SEM uses an approximation based on polynomials of a higher degree (\(p\gt 1\)). The solution is expressed as a linear system of equations that can be solved locally for each element. An element solution \(e\) can be represented in three dimensions as a tensor \(v_{ijk, e}\) with \(i,j,k \in \lbrace 0, \dots , p\rbrace\). Often, the polynomial degree \(p\) is the same for all spatial dimensions. All elements operate on independent tensors and can be elaborated in parallel.

In our example, we focus on solving the Helmholtz equation \(\lambda u - \nabla ^2 u = f\) for quadrilateral elements. The Inverse Helmholtz operator subsumes simpler operators (e.g., interpolation) that are similarly relevant in CFD simulations [4]. The operator can be formulated as

The CFDlang DSL for tensor operations [1] enables us to encode these kinds of operators concisely. We have chosen this particular DSL because of its excellent domain capture and limited optimization scope. CFDlang is target agnostic, and its user interface is close to the mathematical problem specification, reducing the programming burden on application developers. Figure 2 shows a CFDlang program implementing the Inverse Helmholtz operator, where lines 7–9 are the direct transcriptions of the expressions (1a)–(1c). The CFDlang description does not define exact implementation details, allowing the compiler to optimize the operations for a given target. Another implicit part of the given DSL program is that it is assumed to be applied to all the independent elements in an implicit outer “element loop.”

From the application developer viewpoint, once the kernel is specified in CFDlang, it can be integrated into larger Fortran or C++ code applications via their respective foreign function interface mechanisms where we can also embed the FPGA runtime library calls.

Previous work found that typical software implementations achieve performances between 1 and 16 GFLOPS based on the polynomial degree \(p\), with an average power consumption of at least 100 W [1]. GPUs, in turn, do not feature as good a scaling behavior for these kernels [26].

HBM is a novel memory architecture that enables high-performance and adaptability for memory-bound applications [20]. HBM is a three-dimensional (3D)-stacked DRAM that offers high-bandwidth and energy-efficient data movements. The logic die is connected to the HBM die(s) with through-silicon vias. The rest of this work focuses on the Xilinx Alveo U280 data center accelerator card. However, other HBM-based FPGA devices (like the Intel Stratix 10 MX) are conceptually similar [31]. The Alveo U280 is built on the Xilinx 16nm UltraScale+ architecture and offers a rich set of memory solutions, as shown in Figure 3. The Alveo U280 card features the XCU280 FPGA, which combines three super logic regions (SLRs). An super logic region (SLR) is a physical section of the FPGA with a specific amount of resources and connections, as shown in Table 1. SLR0 integrates an HBM controller to interface with the HBM2 subsystem through 32 pseudo-channels (PCs), each with direct access to 256 MB of storage (8 GB in total). Each 256-bit PC operates at 450 MHz, yielding a maximum bandwidth of 14.4 GB/s. The full system can thus achieve a theoretical bandwidth of 460.8 GB/s. SLR0 also connects to the host via 16 lanes of the PCI Express (PCI-e) interface. SLR0 and SLR1 each connect to 16 GB of DDR4 each. Finally, each region has up to 8 MB of Programmable Logic RAM (PLRAM) for fast access to small datasets. In the following, we will use the term global memory for the set of memories available on the board. The host must transfer data into the device global memory before they can be accessed by the FPGA logic.

The target system for the Alveo U280 is composed of multiple compute units (CUs). Each CU is a user-defined hardware module that can be attached to any of the PCs through independent AXI interfaces, while the built-in HBM controller and switch have access to all physical channels. The CU can be described in C++ and synthesized with HLS or specified directly in RTL. Multiple CUs allow parallel execution but must be connected to different HBM channels. The system configuration file describes the connections between the CU ports and the HBM channels. The required logic is automatically generated during system synthesis.

Xilinx offers a unified software platform, called Vitis, to develop FPGA applications. Vitis includes a rich set of hardware-accelerated open source libraries optimized for Xilinx FPGA and the Xilinx Runtime library (XRT) to facilitate communication between the host application (running on the host CPU) and the accelerator deployed on the reconfigurable portion of the card, which is connected via PCI-e. It also includes user-space libraries and APIs, kernel drivers, and board utilities that can be used to measure performance and monitor power consumption. In this work, we aim to automate the generation of CU descriptions and the associated configuration file directly on top of the existing Xilinx libraries.

In modern FPGA data center cards like the Xilinx Alveo U280, DDR4 memory is excellent for accessing large datasets with modest latency, but the transfer bandwidth is limited to 36 GB/s and no more than two parallel accesses [25]. Conversely, HBM has slightly higher latency, but it can reach a much higher theoretical bandwidth (460.8 GB/s) when fully using all its 32 PCs [31, 52]. HBM-based architectures are becoming popular for memory-bound applications [11, 28]. They have been used for convolutional neural networks, graph analytics, or weather prediction [36, 50, 53, 60]. These platforms also feature PLRAM, which can be used for small amounts of frequently accessed data. In the case of the Alveo U280, the designer can configure these connections in the system configuration file. However, several factors can affect the performance of the resulting system.

Challenge 1: CPU-Host Communication Cost. Moving data across the PCI-e interface to the global memory has a high cost compared to the data access from the logic regions to the global memory. So, such movements must be optimized to reduce their latency. First, the designers must move enough data to perform significant computation, for example by operating on multiple elements in sequence before sending back the results. In many cases, the kernel must read the same data multiple times. For example, the Inverse Helmholtz operator must read matrix S several times during the two tensor contractions. To avoid multiple global memory accesses to the same data, designers can use internal buffers that gather enough data.

Challenge 2: Read/Write Burst Transactions. Vitis can automatically infer the size of data transfers to create more efficient burst transactions. However, frequently switching between read and write transactions is inefficient due to memory controller timing parameters. So, transactions should be ordered to maximize the data movement in one direction before switching to moving data in the other direction. The alternative is to partition reads and writes into separate HBM channels, at the cost of increasing the number of channels required for each CU.

Challenge 3: Full Bandwidth Utilization. The AXI interfaces to the HBM channels can be configured with wider data busses. For example, the designer can use a bus width of 256 bits running at 450 MHz or even 512 bits running at 225 MHz. This allows the CU to exchange more data in each clock cycle, reducing the read/write cost from/to the HBM channels. For example, when configuring the 32 parallel HBM channels with a bus width of 256 bits, the potential bandwidth between the CPU and the FPGA is about 460 GB/s at 450 MHz. However, such data must be processed efficiently and in parallel to avoid performance bottlenecks [28].

Challenge 4: HBM-Data Allocation. To avoid congestion in the switch, it is important to partition the data into the different HBM memory regions such that each CU uses as few channels as possible and shares these channels with as few other CUs as possible. Otherwise, the designer must introduce an efficient crossbar to mitigate congestion [15].

Challenge 5: Synthesis-Related Issues. FPGA devices include hard macros like DSP and BRAM for efficient computation and storage, respectively. However, the physical location of these FPGA resources can affect routing. Also, the XCU280 FPGA, which is the device at the basis of the Alveo U280, is partitioned into three SLRs, but only SLR0 is connected to the HBM channels (cf. Section 2.2). When the CUs cannot fit into a single SLR, they are automatically split over multiple SLRs using special resources called super long line routes. These resources allow routing between CUs in any SLR, enabling access to all memory resources, but they introduce performance overhead. So, if a CU requires access to the HBM, then it should be allocated in SLR0. Also, hardware modules that access multiple HBM channels should be close together to reduce wire length.

MLIR [35] is a recent compiler infrastructure that allows defining custom abstractions and high-level transformations. It has received considerable traction as a framework for developing domain-specific compilers for heterogeneous systems, especially for machine learning frameworks. MLIR is in itself not a fixed intermediate representation (IR), but an infrastructure with a central plug-in mechanism to extend the vocabulary of the intermediate representation (IR) using dialects. A dialect implements an intermediate abstraction as a set of operations, types, and attributes. Examples of prominent dialects include linalg for linear algebra operations, affine for nested loop programs, and the llvm dialect that enables a seamless transition into the LLVM compiler framework at the lowest level of abstraction. Custom dialects allow compiler designers to develop analyses and transformations at the right level of abstraction. Dialects can be integrated into larger language stacks, fostering the reuse of abstractions and code transformations. A typical way of integrating a dialect is using lowering, where a more abstract dialect is converted into a more concrete one through specific transformations. Arguably the most complete flow built on top of MLIR to date is presented in Reference [58], where authors describe how MLIR allows for a decentralized and composable compiler.

In this article we leverage MLIR to create a composable and reusable domain-specific compiler. We implement the abstractions for tensors described in Section 2.1 and propose a lowering pipeline that integrates with existing dialects in MLIR, while exposing profitable transformations for HBM-optimized designs. We also show how hardware-specific abstractions, e.g., for number representations, are integrated into the compilation flow.

As FPGA devices are entering data centers, HBM architectures are becoming extremely popular to accelerate a variety of data-intensive applications [27]. Benchmarking these architectures defined that HBM provides higher bandwidth and access latency than traditional DDR [25, 37, 39, 60]. Instead, comparing FPGA performance with both CPUs and GPUs in HPC workloads concluded that there is potential in this technology. They also point out that although FPGAs struggle to compete in absolute terms with GPUs, they achieve much greater energy efficiency [11, 30, 42]. References [15, 20] argue that one of the main problems in FPGA is the absence of a comprehensive memory hierarchy, like the one in CPUs. For this reason, they propose two different custom crossbar solutions that efficiently access the HBM and create an additional memory layer using the BRAMs available inside the FPGA. However, these solutions restrict the amount of BRAMs available to implement the user’s application and they do not fully exploit the FPGA pipeline parallelism. Other works like References [24, 31, 53] define a series of guidelines to achieve high performance. Specifically, the designer must always instantiate the maximum number of computing units to exploit HBM parallelism, use 256-bit packets to maximize the bandwidth of the HBM AXI port, exploit the algorithm parallelism through hardware pipelining, and finally execute the task in a dataflow manner [28]. However, none of these works evaluated these optimizations together, along with the effects when scaling the number of computing units. None of the papers presented so far defines methodologies at the HLS level that can be applied to design to improve its performance, especially for tensor applications. Reference [17] proposes three buffer restructuring approaches, coarse-grained, fine-grained, and hybrid, to adjust the BRAM and lookup table (LUT) utilization as needed. Similarly, tensor optimizations were proposed in Reference [54]. Reference [46] studies the effect of HLS pragma directives on the design. In particular, they indicate that loop unrolling, and array partitioning directives can cause overhead in off-chip memory performance due to non-burst access. Finally, they provide a code transformation technique to exploit the data width of the memory controller. The Merlin compiler¹ provides an infrastructure for source-to-source transformation to accelerate applications on FPGA that are already described in C++ [16]. Instead, operating at the DSL level allows us to implement high-level transformations to the source code, including exploring different memory layouts.

The core of our DSL abstraction consists in modeling tensor expressions. There are several such DSLs, like TensorComprehensions (TC) [59] and TensorFlow Eager [6], that are capable of representing tensor operations in a more general domain. However, many are based on backing software libraries, such as TensorFlow [5] and Theano [9], which are less flexible in terms of their back end. TVM [13] is a compiler focusing on machine learning (ML) that is not as limited in that regard, even providing an FPGA back-end. Its focus on machine learning (ML) makes it hard to apply TVM to other application domains. A similar application-specific flow is provided by ESP4ML [22], which explores a bigger design space of systems-on-chip using HLS but does not provide automatic generation from tensor expressions. More generally applicable, the Halide [45] ecosystem provides extensive support for heterogeneous systems. Spatial [29] provides automatic compilation and system generation from high-level descriptions. The framework offers high-level communication libraries, but the application development burden is still on the designer.

Thanks to MLIR [55], as a framework for designing intermediate languages, several projects have recently proposed programming flows for heterogeneous systems. In the MLIR infrastructure, compiler internals are composed using a plug-in system, which allows connecting various front-, middle- and back-ends, which may be entirely orthogonal. For example, with Teckyl,² the aforementioned TensorComprehensions (TC) has received a front-end for the MLIR compiler infrastructure, enabling a variety of back-end consumers. As an alternative to just lowering to LLVM [34], IREE³ offers perhaps the first functional end-to-end flow entirely designed within MLIR. Closer to our domain is the recent Open Earth Compiler [23] for efficient GPU code generation for stencil operations. We also rely on the plug-and-compile MLIR philosophy, which will allow others to reuse our abstractions and allows us to profit from existing language stacks and compilation flows.

In the context of MLIR, projects such as Polygeist [40] and Phism [64] use polyhedral modelling to ingest, optimize, and even emit C code for HLS. Our compiler still employs polyhedral analysis and rescheduling as described in Reference [19]. Streaming implementations are a challenge in FPGA implementations, which can tie into the polyhedral model via consecutivity constraints [2].

In conclusion, to address productivity issues in scientific computing, designers need a comprehensive framework that (1) allows them to use high-level tensor expressions to specify the behavior, (2) automates the hardware generation process on top of commercial HLS solutions, and (3) efficiently targets modern HBM architectures.

To simulate the complete volume, the CFD application applies the Inverse Helmholtz operator to \(N_{eq}\) independent elements. To exploit such intrinsic parallelism and address the challenges discussed in Section 2.3, we aim to build a target system like the one in Figure 4.

From the system-level perspective, the CFD simulation runs on the host CPU, which sends the data to the FPGA HBM via PCI-e. Once the data are in the HBM, the CU can fetch them in parallel through AXI channels for several simulation elements. Since each HBM interface is 256 bits wide, the channel can be conceptually divided into multiple lanes based on the data bitwidth. Each lane can serve an independent accelerator, called Kernel, derived from the high-level DSL description, i.e., the Inverse Helmholtz operator of Figure 2. The CU is the fundamental unit that can be directly attached to the HBM channels and can be internally composed of several kernels with a wide range of communication patterns depending on the high-level description of the functionality to be implemented. To manage data exchanges between the HBM interfaces and the parallel kernels, the CU uses two additional hardware modules, i.e., Read and Write in Figure 4, that execute in parallel to the kernels and communicate with a dataflow model. The Read module fetches the input data of one element (matrices S, D, and u) from the HBM memory into internal buffers, while the Write module transfers the corresponding results (matrix v) into the HBM memory. To avoid conflicts and reduce the complexity of the AXI interfaces, we assign these modules to independent HBM channels. The data read and write modules within the CU split the 256-bit HBM data into 32-bit or 64-bit data for the computation modules based on the data format.

To hide the communication cost, the host computes the maximum number of sets of input data that can fit in one HBM channel. We thus define a batch as the number of elements \(E\) that the CU kernels will elaborate before the host retrieves the output. Executing a batch of several elements maximizes the size of the host data transfers, minimizing their associated overhead. In this way, the data exchanges between the CPU and the global memory can be overlapped with the CU execution using a double-buffer method. While the host is exchanging data with the Ping (Pong) channel, CUs can operate on the Pong (Ping) channel. Given the CU structure, we can compute the number of batches \(N_b=N_{eq}/E\) to be executed based on the total number of elements \(N_{eq}\) to be simulated.

To further exploit parallelism, we can then replicate the CU structure multiple times, each of them operating on independent HBM channels. Let \(N_{cu}\) be the number of CUs that can be instantiated based on the resource constraints, the host application executes \(I=N_b / N_{cu}\) iterations, evenly distributing the number of batches over the available CUs.

In the following, we describe our DSL-to-bitstream flow to automatically create this architecture on top of the existing HLS platform (e.g., Vitis for the Alveo U280). We also describe how to automate the integration of the optimizations to increase performance and resource efficiency and how to modify the host code to match the memory layout required by the hardware modules.

We propose a modular tool flow that simplifies the creation of FPGA accelerators for tensor-based applications, like numerical simulations, that are expressed in domain-specific languages. Figure 5 shows an overview of our flow that is composed of two major steps. We start from the DSL description of the kernel and generate the corresponding optimized and HLS-friendly C code (DSL-to-C generation). The compiler performs a series of transformations at different levels of abstraction (cf. Section 3.4), inserts pragmas to guide the HLS flow, and produces metadata for the memory generation. Starting from the C code of the kernel, we create the corresponding parallel system with multiple CUs (cf. Section 3.5), each of them connected to one or more PCs (C-to-system generation). Also, each CU can instantiate one or more kernels based on the required amount of FPGA resources (cf. Section 3.6), along with the logic to move data from the global memory to the kernel on-chip buffer (Read and Write modules). We specify the CU functionality in C++, wrapping the kernel code with specific code and HLS directives to enable the generation of the memory-optimized architecture. To do so, we use a combination of distinct tools: a redesign of the CFDlang compiler (cf. Section 3.3), which integrates novel kernel-level transformations, the Mnemosyne tool [44], which improves memory sharing inside the single kernel to reduce on-chip memory requirements, and the newly developed Olympus tool, which generates the system architecture (i.e., CU description and configuration file) around the kernels based on the optimization requested by the designer. Olympus also generates the corresponding host software to control the accelerators. This code comprises multiple steps (data allocation, kernel configuration, data transfers, synchronization primitives). Each step is wrapped into a specific function so that it can be further specialized in the hardware generation flow with specific optimizations. We then use the v++ tool of the commercial Xilinx Vitis Unified Platform to create the hardware description with HLS and generate the final FPGA bitstream. The Vitis platform also builds and links the host application with the Xilinx Runtime (g++). Our flow can easily be extended to target similar HBM-based platforms (like the Intel Stratix 10 MX) with the respective toolchains.

Compiling for the template architecture in Figure 4 is fundamentally different to compiling for mainstream multicore CPUs. Instead of extending the original CFDlang compiler [1] with new target support, we devise an entirely new implementation based on the MLIR framework. With MLIR, we raise the language abstraction level while remaining target-agnostic in our front end. At the same time, our middle-end is also mostly agnostic toward the concrete DSL we chose.

This also allows us to profit from the well-engineered and widespread MLIR and LLVM ecosystems, enabling re-use of abstractions and of lowering flows to different architectures. In the following, we describe the new CFDlang compiler infrastructure, focusing on the stack of dialects we created to support CFD simulations (cf. Figure 6). Our infrastructure includes the key additions of only three new dialects: cfdlang, i.e., our language front-end, teil, i.e., our domain optimization middle-end, and base2, i.e., our example for a hardware-specific back-end task. Apart from the dialects linalg and affine mentioned in Section 2.4, the figure also includes other standard dialects, like tensor for cross-domain tensor operations, scf to model structured control flow with, e.g., explicit loops, and external tools such as the integer set library (ISL), which we use to produce the code for the downstream HLS compilation. The transformations performed within these dialects are discussed in Section 3.4.

Since we aim at leveraging a commercial HLS tool that does not currently support MLIR as input, our DSL compiler is an MLIR-based transpiler that emits C99 source code and additional metadata outputs to directly interface with Vitis HLS and Mnemosyne, respectively.

In Reference [19], we described analysis and transformations at the kernel level that allowed for memory allocation and operation scheduling for FPGA offloading. These were based on polyhedral techniques, which still apply to the affine intermediaries that our MLIR flow produces. For example, the liveness analysis required for Mnemosyne’s sharing optimization follows the same approach described in Reference [19]. To address the HBM challenges (cf. Section 2.3), we introduce new transforms at higher levels, namely operator graph partitioning and pipelining. These increase the achieved memory bandwidth required to take advantage of the high throughput of the HBM-based target architecture. We elaborate on Section 3.3 by discussing the individual steps as shown in Figure 9.

Our hardware generation flow aims at optimizing the data transfers around the kernel implementation produced by CFDlang. It receives as input the description of the kernel and the compatibility graph of the internal buffers (from CFDlang), along with information about board resources (from the user). It produces an optimized CU description (in C++) and the platform configuration file based on the number of CUs that can be instantiated. The configuration file also specifies the proper connections to the HBM channels. To implement our hardware generation flow, we developed Olympus, which creates an optimized system-level architecture (in synthesizable C++) for our accelerators. During its execution, it interfaces with Mnemosyne [44] to optimize the on-chip memory associated with each kernel and limit the number of local memory resources. Mnemosyne uses the buffer compatibility graph generated by the CFDlang compiler to determine when the physical on-chip memory banks can be reused without performance overhead [44]. After this, Olympus reads the kernel interface and determines how to connect the input/output ports to the rest of the system to efficiently exchange data with the HBM channels. Data ports are connected to HBM channels via AXI Master interfaces, while configuration ports are connected to the host via AXI-Lite, memory-mapped interfaces. Data exchanged with HBM channels (i.e., input and output matrices) are buffered on-chip to allow fast, fixed-latency access during kernel execution.

The generation of the system architecture is guided by the designer, as shown in Figure 5. Starting from the C kernel produced by the CFDlang compiler, Olympus generates a minimal wrapper to run HLS and obtain an estimate of the resources needed for the target FPGA device. Then, we apply optimizations to both the on-chip local memory of each kernel (with Mnemosyne) and the system-level memory architecture to exchange data with the HBM (with Olympus).

At the kernel level, we use Mnemosyne to generate the RTL of the on-chip memory architecture associated with each kernel. Mnemosyne is a tool that exploits sharing compatibilities, i.e., when distinct internal buffers have no overlapping lifetime, to assign them to the same physical banks based on a given cost metric. In these cases, the tool generates custom logic to manage the accesses to the same memory banks from different kernel interfaces [19, 44]. Figure 13 shows the conceptual interface, described in C, of a kernel before and after the execution of Mnemosyne. Mnemosyne requires the specification of the compatibility graph of the local arrays and, in our flow, these metadata are computed and produced by CFDlang during liveness analysis [19]. Mnemosyne wraps the RTL kernel description (produced by HLS) with the resulting RTL description of the kernel memory architecture to expose only input and output ports to the CU. This conceptual interface is then used to integrate the kernel into the CU.

At the system level, Olympus generates the C++ description of the memory architecture around the kernel description, which is integrated as custom RTL to use the results produced by Mnemosyne. Since the hardware cost of the kernel may limit the number of parallel CUs, the designer can use Olympus to understand which optimizations can be applied given the available FPGA resources. Indeed, we characterize each optimization with an estimation of the extra resources. With this information, the designer can select the most suitable optimizations, and Olympus generates the corresponding CU description around the CFDlang-generated code of the kernel and the system configuration file. In the future, this process can be further automated by combining it with state-of-the-art design space exploration frameworks. Each CU can feature multiple kernels, each of them connected to a lane to fully utilize the AXI bandwidth (cf. Section 3.1). The CU wrapper implements data-movement optimizations and is designed accordingly with changes to the host application and the configuration file. For example, the kernel may benefit from a change in the way data are written to and read from global memory and therefore the host application must be updated accordingly. The configuration file, instead, defines how each CU interfaces with the HBM. By modifying this file, Olympus defines how each CU is connected to the individual channels.

The resulting components are then passed to Vitis, i.e., the HLS platform that we use, to automatically generate the bitstream required for board configuration. When timing is not met, Vitis automatically downscales the execution frequency. While this is useful to enable proper functionality, the designer has little control over this process. Conversely, Olympus introduces some optimizations for the synthesis process (cf. Section 3.6.4).

Let \(p\) be the polynomial degree for the simulation, each Inverse Helmholtz kernel needs \(p\times p\) values for matrix S and \(p\times p \times p\) values for matrices D and u. The kernel then produces \(p\times p \times p\) values for the output matrix v. These values are useful to estimate the number of resources for input, temporary, and output buffers. The size of the total amount of data used to compute one element can be used to determine how many elements worth of data can fit into an HBM channel (max size is 256 MB). The number of elements is the batch size, i.e., the number of executions the CU can perform without interruption before needing the host to send more data. Data transfers are required between consecutive batches.

In the following, we describe the optimizations that we apply to our CFD application by means of the Olympus options (step Optimize in Figure 5). Such optimizations (and the Olympus hardware generation process) can be easily adapted to similar tensor-based applications.

To evaluate our DSL-to-bitstream flow, we extended the flow presented in [19] as described in Figure 5: CFDlang is implemented on top of the MLIR infrastructure, Mnemosyne [44] is an open source tool,⁹ and Olympus is a new in-house prototype. Olympus is built in Python on top of the Pyverilog library [57] for hardware generation (i.e., the generation of the kernel wrappers around Mnemosyne artifacts) and the Pycparser library¹⁰ for code generation. With our novel flow, we targeted a Xilinx Alveo U280 card on an AMD EPYC 7282 [7] server running Centos 7. We used Xilinx Vitis 2021.1 [63] for synthesis and bitstream creation. Unless otherwise specified, we target a synthesis frequency of 450 MHz for both the platform and the CU description. For each implementation, we evaluate performance and energy efficiency by using the GFLOPS and GFLOPS/W metrics, respectively. The GFLOPS metric is obtained by dividing the total number of floating-point operations executed by the application by the total execution time, while the GFLOPS/W metric is obtained by dividing the GFLOPS metric by the average power consumption of the system. To get accurate information about FPGA power consumption, we profile the power consumption during the system execution with Xilinx XRT, and we use the average value.

To evaluate the cumulative benefits introduced by each hardware optimization, we progressively apply them to the Inverse Helmholtz operator used as a case study and so extensively discussed in this article. We evaluated the implementations for two polynomial degrees: \(p=7\) and \(p=11\). In particular, given the polynomial degree \(p\), we assume that each contraction is composed of three loop nests that execute two floating-point operations (one addition and one multiplication) for \(p\times p \times p~ \times ~p\) times each. Similarly, the Hadamard product requires \(p\times p \times p\) multiplications. So, the entire Inverse Helmholtz operator the following number of floating-point operations:So a single element is required to execute \(N^{el}_{op}\) = 177,023 floating-point operations when \(p=11\) and \(N^{el}_{op}\) = 29,155 floating-point operations when \(p=7\). The total number of floating-point operations for a CFD simulation can be obtained asWe executed all experiments with \(N_{eq}= \hbox{2,000,000}\), i.e., we simulated 2,000,000 elements.

We executed our CFDlang on the DSL description in Figure 2 to generate the C kernel for hardware optimization. We first performed experiments to evaluate the effects of optimizations with \(p=11\). In particular, we progressively added the following optimizations:

For each of these implementations, we measured total and kernel execution times, maximum and average power consumption, and cost in terms of hardware resources. Figure 15 shows the performance (in terms of GFLOPS) achieved in each experiment when adding the specific optimization on top of the previous ones. In each experiment, the left black and white bar (CU) shows the GFLOPS of the CUs on their own, without considering host-FPGA data transfers, while the right azure bar (System) includes the entire application. Comparing the two bars allows us to evaluate the peak performance of the kernels and the effects of data transfers.

The Baseline case achieves only 2.9 GFLOPS, and the difference between the CU performance and the overall system performance is 9.2%. This is due to the serial nature of the implementation where data are transferred from the host to the HBM, then processed by the CU and sent back to the host before starting a new batch. If more data need to be transferred, then this discrepancy between CU performance and overall system performance will grow larger, as the CU needs to wait for all of the data to be sent before beginning execution.

After the Double Buffering optimization, the CU performance remains similar, with a small degradation due to overhead, while the system performance is now the same as the CU performance. This is an improvement over the Baseline implementation because now the host to HBM data transfers are happening in parallel to and are entirely hidden behind the CU execution.

We then executed two experiments for evaluating Bus Optimization. In the Serial version, we attempt to utilize the full bandwidth of the 256-bit bus by packing four doubles. The CU reads them in parallel, but then it serializes them when it needs to access its own local buffers. While this optimization is supposed to speed up data reads from the HBM, its implementation in the CU leads to a performance degradation of about 3\(\times\). This is mostly due to the complexity of aligning the data (\(p\times p\) and \(p\times p\times p\)) to multiples of fours inside the bus. To combat this, but still use the full bus bandwidth, this optimization was replaced with the Parallel implementation where four kernels are instantiated in the CU and the data for each “lane” is stored in separate buffers, one for each kernel. This led to a \(3.92\times\) speedup over the Serial implementation, as the four lanes are now all read in parallel. There is only a \(1.23\times\) speedup over the Double Buffering implementation, where a \(4\times\) speedup would be expected. This is because in both Bus Opt implementations, the HLS tool only instantiated two double-precision multipliers (rather than 11 in all other implementations). So when the innermost loops are unrolled, a resource limitation violation increases the initiation interval (II) to 4, effectively counteracting the expected \(4\times\) speedup. We use the Parallel architecture in the following experiments.

Next, we tested various forms of the Dataflow Optimization. Each implementation of this optimization separated the kernels into read, compute, and write modules and streams were used to pass data between them, allowing a pipelined structure. When using one compute subkernel (Dataflow (1 Compute)) test, the speedup was \(3.68\times\). In these cases, the HLS tool instantiated the full 11 multipliers, removing the II violation. This decision of the HLS tool along with the overlapping execution of the read, compute, and write modules allows us to effectively achieve a 4\(\times\) speed-up, i.e., to full exploit the four parallel lanes. Since the compute module dominated the execution time, it was further split into two, three, and seven modules. The Inverse Helmholtz operator comprises seven loop nests, each implementing the operations in the seven grey rows on the right side of Figure 10. To split into two modules, the kernel is divided into a first module with the first three loop nests with S and u as input and t as output, and a second module with the last four loop nests with S, D, and t as input and v as the final output. The split was made as such so that the first module does not need D as an input and the second module does not need u as an input. To split into three modules, we use the division shown on the left side of Figure 10 and in Figure 11. This is the most natural division as it matches the initial DSL representation where the first three loop nests implement the gemm operator, the fourth loop nest implements the mmult operator, and the last three loop nests implement the gemm_inv operator. Another benefit to this division is that the mmult loop nest consumes and produces data in the same order it is sent via the streams, meaning that no extra buffering is needed for this module and each data element can be immediately processed as it is received leading to a minimal latency. To split into seven modules, each loop nest is a separate module.

All three of these tests gained speedup over the 1-Compute version by breaking the total execution time of a single module down further. The 2-Compute version is \(1.7\times\) faster than the 1-Compute version. The discrepancy between this result and the ideal speedup of \(2\times\) is due to the extra data buffering that must be done in each module to allow random access. In the 1-Compute case, the input arrays are each buffered one time, which adds an overall latency equal to the total input data size. In the 2-Compute case, the S array is needed by both modules and must be buffered twice. Additionally, the output of the first module, t, is used as input to the second module and must also be buffered. This extra buffering is overlapped while the two modules execute in a pipeline fashion, but it means that the latencies of each module are not exactly half of the total latency of one unified module. However, 3-Compute modules was slower than 2-Compute modules. The overall execution time is determined by the module with the longest latency, as it is the limiting factor in the overall latency of the system. Because the loop nest implementing the gemm operator has a minimal latency, moving it to a separate module does not significantly reduce the latency of the largest module. In fact, in each case, the module with the longest latency was the same, but the extra modules and control routing caused the tools to frequency scale the 3-Compute case to execute at 266 MHz, whereas the 2-Compute case executed at 292 MHz. When this is considered, the performance of both tests is approximately the same. The 7-Compute test, however, performed the best because each of the compute modules were much smaller than the previous tests. In this case, the latencies of these modules were now slightly shorter than the latency of the read module, meaning that this is the limit of the performance increase by dividing the compute portion. The 7-Compute test gained a total speedup of \(4.03\times\) over the Bus Opt Parallel implementation.

To evaluate the efficiency of the allocated resources, we computed the “ideal” GFLOPS value for each of the double-precision floating-point implementations. We identified the total number of double-precision adders and multipliers instantiated in the CU (# Ops) by analyzing the synthesis reports generated by Vitis HLS. The ideal GFLOPS is computed by multiplying this value by the frequency of the CU and represents the performance when all operators were constantly in use concurrently. Table 2 compares these values with the measured GFLOPS of each implementation. In the last column, an “efficiency” is calculated as a ratio between the ideal and achieved GFLOPS.

This “efficiency” reflects the behavior of the allocation and scheduling of the HLS tool more than the efficiency of our system design surrounding each HLS kernel. However, we can still gain some insight into our design in the cases where the HLS decisions were the same. For instance, in the Baseline and Double Buffering cases, the same kernels are used, and therefore they each have the same # Ops. The efficiency increases because less time is “wasted” waiting on data transfers from the host. Both Bus Opt implementations reduce the # Ops because the HLS tool used a different local memory type with fewer read ports, restricting the unrolling, and therefore only used two adders and two multipliers for each kernel. The efficiency values of these implementations are also much higher because these are the only cases where the multipliers themselves are pipelined. The ideal GFLOPS metric expects each operator to produce a result every clock cycle, so in all other cases where the operators are not pipelined, there are several cycles of latency for each operation reducing the efficiency. Between each of the Dataflow implementations, the efficiency drops slightly as the computation is split into more modules because it is impossible to split the computation into equal latency modules. The module with the longest latency may constantly be computing, but the shorter-length modules must stall.

The efficiency values for all implementations (except Bus Opt) are all near 0.5 because each multiply-accumulate is implemented as eleven parallel multipliers and eleven sequential adders. Even though the additions are sequential, the tool still allocated eleven of them. Because the Bus Opt implementations are restricted to two adders, their efficiencies are higher.

At this point, we want to start replicating the CUs using the remaining area available in the FPGA fabric to maximize parallelism. For this reason, we need to evaluate the hardware cost of each implementation. The numbers of LUT, FF, BRAM, URAM, and DSP used by each case for \(p=11\) are shown in Table 3. In general, each test from Baseline to Dataflow (7 Compute) showed an increase in resource utilization. Any utilization value over 25% is shown in red. These resources are most likely to cause placement and routing issues when instantiating multiple CUs. We tested a few methods to reduce resource utilization to be able to increase the number of instantiated CUs.

The Mem Sharing optimization is applied to the Dataflow 1-Compute implementation where several arrays are used in the compute module (cf. Figure 14(d)). Mnemosyne generated an architecture to internally share arrays based on their liveness intervals. This decreased the BRAM utilization by 14.5% and the URAM utilization by 48.3% while the LUT and FF utilization only increased minimally and the DSP utilization remained the same. Also, the execution time was only slightly reduced (a slowdown of \(0.98\times\)). Conversely, this optimization cannot be applied to the Dataflow 2-Compute, 3-Compute, and 7-Compute implementations because, in these cases, each compute module only uses arrays that cannot be shared, as they are always in use during the module execution. This optimization is indeed beneficial when on-chip memory inside the CU is the limiting factor, and when replicating the CUs brings more improvements than dataflow execution.

Another method to reduce resources is to change the numerical representation. All of the previous tests used the floating-point format with double precision. In general, fixed-point representations utilize fewer resources than floating-point ones. We tested 64- and 32-bit fixed-point representations by modifying the Dataflow 7-Compute implementation. The 64-bit implementation uses 24 bits for the integer portion and 40 bits for the fractional portion. The 32-bit implementation uses 8 bits for the integer portion and 24 bits for the fractional portion. These values are provided by the user after an analysis of the algorithm. Because the 32-bit data are half the size, we instantiate eight kernels per CU and divide the 256-bit bus into eight lanes. In the Fixed Point 64 test, the LUT utilization reduced by 46.3%, the FF utilization reduced by 53.4%, the RAM utilization remained the same, and the DSP utilization increased by 44.8%. In the Fixed Point 32 test, concerning the Fixed Point 64 test, the LUT and FF utilization remained roughly the same. The DSP utilization was nearly halved. The BRAM increased by about four times while the URAM decreased to zero. This is because the data representation is half as long, so the overall size of the data structures is half as big. The arrays representing the tensors are no longer big enough for the synthesis tool to decide if it is efficient to use URAM to store them. When considering the size of the physical memories, the total memory space is approximately halved. The performance of the Fixed Point 64 test had a slight speedup of \(1.19\times\) due to the simplification of the logic allowing the frequency to be higher. The Dataflow 7-Compute test with double format was scaled to 199 MHz while the Fixed Point 64 test was scaled to 234 MHz. The performance of the Fixed Point 32 test had a speedup over the double format of \(2.37\times\) and it reaches up to 103 GFLOPS. This represents a speed up of more than 35\(\times\) over the Baseline version. The Fixed Point 64 test exhibited a mean square error of \(9.39\times 10^{-22}\) while the Fixed Point 32 test had a mean square error of \(3.58\times 10^{-12}\). It is up to the application designer to determine what an acceptable error is and decide on an appropriate number format, and our flow can help facilitate a design space exploration of these parameters.

Another method to reduce resource utilization for this kernel is to vary the input parameter \(p\). We tested the Dataflow 7-Compute implementation using 64-bit double, 64-bit fixed point, and 32-bit fixed point with \(p=7\) and \(p=11\). The results are summarized in Figure 16 and Table 4.

Compared to their \(p=11\) counterparts, the \(p=7\) implementations performed slightly slower. This is because the actual hardware implementation does not scale the same way as the conceptual floating point operations per kernel (used to compute the GFLOPS). However, the resource reduction between \(p=11\) and \(p=7\) is enough to allow for more replication of the CUs. For instance, the Fixed Point 32 implementation uses 66.4% of the available BRAM for \(p=11\) while it only uses 21.7% for \(p=7\), allowing \(4\times\) replication.

To further facilitate instantiating multiple CUs, we reduced the stream FIFOs from a naive full size to a small enough version to save space and still prevent deadlock. This led to a small performance reduction due to stalls but significantly reduced the total number of BRAMs. Also, because the DSP utilization was exceptionally high in some cases, we used pragmas to guide the HLS tool on using LUTs instead of DSPs to implement fixed-point multipliers. We used this pragma in one of the seven compute modules to shift some of the resource load off of DSPs and onto LUTs.

We were able to instantiate two parallel CUs for the cases of Double with \(p=11\), Fixed Point 64 with \(p=11\), and Fixed Point 64 with \(p=7\), three CUs for the cases of Double with \(p=7\) and Fixed Point 32 with \(p=11\), and four CUs for the case of Fixed Point 32 with \(p=7\). The performance results for these implementations are shown in Figure 17, and the area results are shown in Table 5. All of these implementations were built targeting 225 MHz, as most of their 1 CU counterparts could not even achieve this.

In most cases, replicating the CUs actually led to a slowdown. This is because the extra logic and routing caused the maximum frequency to be reduced thereby slowing down everything in the system. However, most cases did show speedup in terms of the CU execution time. In particular, the Fixed Point 32 implementations achieved up to 172 GFLOPS for the kernel but around 87 GFLOPS for the system. This huge discrepancy is because even though several CUs are now executing in parallel, all of the data must still be sent from the host to the HBM in series. Host data transfers are now the dominating factor by far, so it is not recommended to replicate CUs until the host data transfer time can be reduced. Otherwise, the overall system will have a slowdown from the extra logic.

From the resource utilization results, it can be seen that both 64-bit data types are constrained by resources used for computation, namely LUTs and DSPs. The 32-bit fixed-point implementation is also somewhat constrained by DSPs. In any case, this application is composed of almost entirely floating or fixed point multiplications, and performance-optimized designs will quickly use most of the available DSPs. The 32-bit cases are also constrained by the on-chip memories, the 3 CU implementation of Fixed Point 32 with \(p=11\) even uses 100% of the URAM, but both Fixed Point 32 implementations were able to be replicated more than their Fixed Point 64 counterparts, due to the data width reduction. The \(p=7\) tests were also, in general, able to be replicated more than their \(p=11\) counterparts due to the effect \(p\) has on the amount of computations and the array sizes. Fixed Point 64 was only able to be replicated twice in both cases of \(p\), as the reduction of DSPs between \(p=11\) and \(p=7\) was not enough to allow for a third CU.

Figure 18 shows the power consumption of the different implementations and a comparison of the energy efficiency (GFLOPS/W for floating-point operations and GOPS/W for fixed-point operations). The bars report the average power consumption measured with the XRT infrastructure. We also include the results of the multiple-CU implementations to show the effects of replication on both power consumption (W bars) and energy efficiency (GFLOPS/W and GOPS/W bars).

As expected, the fixed-point implementations are more efficient than the floating-point ones. Also, reducing the bitwidth from 64 to 32 bits allows us to achieve maximum efficiency. This is because these implementations are much faster and use fewer hardware resources. The \(p=7\) implementations have lower average power consumption than their \(p=11\) counterparts due to their smaller resource utilization. However, in most cases, the efficiency of the \(p=7\) cases is lower due to their longer overall execution time. The multiple-CU implementations are generally less efficient than their single-CU counterparts, both because of the increased work occurring in parallel, yielding a higher average power, and because of longer execution times from frequency scaling.

This section presents the results for two additional kernels. The first kernel performs an interpolation, which maps from \(\mathbf {u}\in \mathbb {R}^{N \times N \times N}\) to \(\mathbf {u}^\prime \in \mathbb {R}^{M \times M \times M}\) via an isotropic operator \(\mathbf {A}\in \mathbb {R}^{M \times N}\). We implemented the Interpolation kernel with \(M=N=11\). The second kernel computes \(\nabla u\), the gradient of \(u\) in all three dimensions. We implemented the Gradient kernel with dimensions \(8\times 7\times 6\). In all cases, we compute the total number of floating-point operations needed for simulating 2,000,000 elements, and we use these numbers to compute the GFLOPS and GFLOPS/W metrics.

Figure 19(a) shows the performance results for the Inverse Helmholtz, Interpolation, and Gradient kernels on the AMD EPYC 7282 and the FPGA. These results include the AMD execution (black bars), baseline (no optimizations) FPGA implementation (green bars), and the fully optimized FPGA implementation (azure bars). The fully optimized kernels all use double-precision floating-point data and implement the Double Buffering, Bus Opt (Parallel), and Dataflow (where each loop nest is a subkernel module) optimizations. The baseline implementations achieved \(10.7\times\)–\(38.3\times\) speedup over their software execution on the AMD EPYC 7282. The optimized FPGA implementations, however, achieved \(36.4\times\)–\(160.2\times\) speedup over the AMD execution.

To compare our results with state-of-the-art software implementations, Figure 19(a) also includes the performance of highly optimized Intel implementations for the Inverse Helmholtz and Interpolation kernels [1] (red bars). These implementations are generated by the original CFDlang compiler using the process described in Reference [1], which was found to outperform expert-crafted manually optimized kernels. The executables are compiled with the Intel Compiler and use the Math Kernel Library 2017.2.174 and were profiled on a 24-core Intel Xeon E5-2680 v3 CPU (Haswell), running at 2.50 GHz. The FPGA-optimized Inverse Helmholtz and Interpolation kernels achieved \(2.7\times\) and \(1.4\times\) speedup over the optimized Intel execution, respectively.

The average power and efficiency for the optimized kernels are shown in Figure 19(b). The estimated efficiencies (GFLOPS/W) of the Intel executions for the Inverse Helmholtz and Interpolation kernels are also shown. These estimations are calculated using the GFLOPS results of the kernel and a conservative estimate of the average power (100 W), assuming the CPU would be operating under a lower load than the thermal design power (120 W). The Interpolation and Inverse Helmholtz kernels are \(4.8\times\) and \(7.0\times\) more efficient than the Intel CPU execution, respectively. Recalling the results from Figure 18, the most power-efficient implementation of the Inverse Helmholtz (32-bit fixed-point with \(p=11\) and 1 CU) is \(24.5\times\) more efficient than the Intel execution.