1 of 29

Development of CUDA-enabled SU2 Code for GPU Accelerated CFD

An ongoing implementation and study of GPU Support in SU2 Code

SEPTEMBER 2024

SU2 CONFERENCE

2 of 29

    • This presentation is based on the work done during the Google Summer of Code ‘24 for the SU2 Code by Areen Raj, Leonardo Cavanha and Ole Burghardt.

    • The presentation will broadly be divided into two categories

Introduction

SEPTEMBER 2024

SU2 CONFERENCE

What we’ve done

What we’ve learnt

    • Changes made to the current code
    • Description of Algorithm
    • Its capabilities
    • Performance
    • Advantage-Disadvantage Analysis
    • Reaction of the code to GPU Support
    • Areas of focus for improvements
    • Skeleton Algorithms that can be used
    • Approach for future implementations

3 of 29

Ryzen 7 5800H (Mobile)

8C/16t, 45W, 16 GB DDR4-3200 RAM

RTX 3060 (Mobile)

80 W (Cap), 6GB GDDR6 VRAM, x8 PCIe 3.0

Before we begin, a few things

    • All changes were made to the SU2 ver. 8.0.1 "Harrier”. The entire pull request can be found at - PULL REQUEST.

    • All the performance tests were done with the following hardware specifications :-

Introduction

SEPTEMBER 2024

SU2 CONFERENCE

    • Members are welcome to run the code on their own machines to correlate results and benchmarks that can offer us better insights as to the performance of the code.

4 of 29

Analysis of the Solution Path of SU2 Solvers

SU2 CONFERENCE

SEPTEMBER 2024

5 of 29

    • Advent of General Purpose Computing for GPU (GPGPU) methodology allows for much faster working of certain codes by switching the hardware on which it executes.

    • GPUs excel at parallel computation, compared to the raw monolithic power of a CPU.

What is the use of such a technique in CFD?

    • PDE solvers spend most of their execution time directly in the linear algebra computations that actually solve the equations. Most of the matrices are sparse in nature.

    • These computations include both level 2 and 3 calls such as Matrix Vector Product, Matrix Matrix Product etc.These calculations can be carried out at much higher speeds on the GPU as it is specialized to handle them.

    • Many CFD programs have employed the use of GPU support - Ansys, Siemens, OpenLB, and specifically FUN3D (7x speedup!)

GPU Programming

SEPTEMBER 2024

SU2 CONFERENCE

6 of 29

    • The following was a run of the FVM NSS Solver with Flat Plate Case which was laminar

SU2 CPU Code

SEPTEMBER 2024

SU2 CONFERENCE

Profiling done through Tracy Profiler

7 of 29

FGMRES - Solution Methodology

SEPTEMBER 2024

SU2 CONFERENCE

    • Time and Space Integration subroutines take up the most of the execution time in a single iteration.

Space Integration

46%

Time Integration

38%

    • The majority of the time spent in the Time Integration subroutine is that in the Linear Solver itself - 96 %

    • And the majority of the time spent in FGMRES is in the Preconditioning and Matrix Vector Product - 98%

8 of 29

FGMRES - Solution Methodology

SEPTEMBER 2024

SU2 CONFERENCE

This is the flame graph for the execution of the solver. Preconditioning takes up a major chunk of the time followed by the numerous Matrix Vector Products.

Our target for the change should satisfy the following criteria

    • Extensive Use in Solution Process
    • Straightforward Implementation
    • Parallel Nature

Our current candidate for treatment will be Matrix Vector Products.

9 of 29

Implementing CUDA Support in SU2 - Matrix Vector Products

SU2 CONFERENCE

SEPTEMBER 2024

10 of 29

Memory and Thread Management

SEPTEMBER 2024

SU2 CONFERENCE

    • The Jacobian Matrix in SU2 is a sparse block matrix.

    • The number of rows of the matrix is equal to the number of points in the domain of the mesh.

    • The number of rows of each block is equal to the variables being solved for and the number of columns is equal to the number of equations. (Block Size = nVar x nEqn)

    • The CUDA API deploys a grid on the GPU which comprises of a number of thread blocks. Each thread block has a fixed number of threads in all three dimensions.

Images taken from CUDA Documentation

11 of 29

12 of 29

Memory and Thread Management

SEPTEMBER 2024

SU2 CONFERENCE

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block 1

...............................

Total Grid and Problem Size

Each thread in a block represents a single point in the domain.

These threads are all put into blocks of 1024 threads each (Max Block Size)

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block 0

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block 1

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block 2

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block N-3

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block N-2

..........

THREAD

0

THREAD

1

THREAD

2

THREAD

3

THREAD

XDim-3

THREAD

XDim-2

THREAD

XDim -1

Block N-1

XDim = 1024/(nVar*nEqn)

N = nPointDomain/xDim

Grid and Block Initialization

13 of 29

Algorithm

T0, I0

T1, I0

T2, I0

T3, I0

T4, I0

T5, I0

T6, I0

T7, I0

T8, I0

T0, I1

T1, I1

T2, I1

T3, I1

T4, I1

T5, I1

T6, I1

T7, I1

T8, I1

T9, I0

T10, I0

T11, I0

T12, I0

T13, I0

T14, I0

T15, I0

T16, I0

T17, I0

T9, I1

T10, I1

T11, I1

T12, I1

T13, I1

T14, I1

T15, I1

T16, I1

T17, I1

SEPTEMBER 2024

SU2 CONFERENCE

temp_var += matrix[matrix_index + (j * nEqn + k)] * vec[vec_index + k]

atomicAdd(&prod[prod_index + j], temp_var)

14 of 29

Yea sure, I totally know whats going on here

SEPTEMBER 2024

SU2 CONFERENCE

15 of 29

    • We introduce a member function for the CSysMatrix Class that carries out the GPU Matrix Vector Product.

    • The input vector and changed matrix is copied to the GPU and then the resultant product vector is copied back onto the host. The matrix is held in pinned memory that has almost twice the bandwidth compared to normal paged transfers.

    • Doxygen documentation is available in the PR.

Implementation

SEPTEMBER 2024

SU2 CONFERENCE

    • To compile with CUDA please use the option -Denable-cuda=true. You will also need to specify the environment variable CUDA_PATH with the location of the installed CUDA Folder - usually found in the /usr/local directory. The decision between GPU and CPU execution is decided at runtime with the help of a new config flag ENABLE_CUDA=YES.

16 of 29

Results

SU2 CONFERENCE

SEPTEMBER 2024

17 of 29

Results

SEPTEMBER 2024

SU2 CONFERENCE

    • All of the computations were carried out for Flat Plate Laminar FVM Case on the NSS Solver with triangular mesh elements.

GPU = 246.477s

CPU = 201.952

Single Core CPU is 18% faster

GPU = 20.72 ms

CPU = 5.478 ms

Single Core CPU is four times faster

18 of 29

Results

SEPTEMBER 2024

SU2 CONFERENCE

19 of 29

So does this mean that there is no feasibility of GPU Acceleration in SU2?

The answer is no, there is a lot of potential.

SEPTEMBER 2024

SU2 CONFERENCE

20 of 29

Performance Analysis and Solution Proposal

SU2 CONFERENCE

SEPTEMBER 2024

21 of 29

MemCpy 1

Jacobian Matrix

1.2 GB/s

MemCpy 2

Input Vector

1.2 GB/S

MemCpy3

Resultant Vector

2 GB/S

Profiling

SEPTEMBER 2024

SU2 CONFERENCE

    • Nsight Systems reveals that memory copy to and from the CPU kills any speedup we get with the call of the CUDA Kernel. (583.26 GFLOPS)

    • In this case the kernel is around 10 times faster than the CPU function. The continuous exchange and transfer of memory slows down the entire process.

22 of 29

Solution

SEPTEMBER 2024

SU2 CONFERENCE

    • Carry out a single copy of matrix and initial vector data at the start.
    • Both the Preconditioning and Matrix Vector Products will take place on the GPU with the device variables.

    • The Gram Schmidt Orthogonalization will also be done on the GPU.

    • What gets copied back is the Hessenberg Matrix per iteration.

    • Finally, at the end, a Memcpy is done to transfer the final resultant vectors back.

23 of 29

SEPTEMBER 2024

SU2 CONFERENCE

New Algorithm

24 of 29

SEPTEMBER 2024

SU2 CONFERENCE

New Algorithm

There are two approaches we can take to port this entire code

Single Kernel

Multiple Kernel

    • Single Kernel handling both Preconditioning and Matrix Operations
    • Uses __syncthreads() barrier synchronization along with extensive shared memory
    • Harder to Debug and Extremely Low Readability
    • Speculation - May have an edge over multiple Kernel
    • Multiple Kernel for each operation
    • CPU Synchronization occurs after each launch as Kernels end
    • Easier to Debug and Extremely High Readability
    • Speculation - May have a slight performance hit due to multiple kernel launch overhead

25 of 29

SEPTEMBER 2024

SU2 CONFERENCE

New Algorithm

We propose a combined approach of the previous two methods

    • First and Second Symmetric Iteration of the Gauss Siedel Preconditioner have been condensed into their own kernels. However, Gauss Elimination has been made its own kernel
    • Mat Vec continues to be its own Kernel
    • Gram Schmit is composed of multiple small kernel launches of Norm, Add and Dot

Ease of porting

    • Preconditioner - Medium
    • Mat Vec - Easy
    • Gram Schmit - HARD

All of these changes and code have been written in a build version dated September 24, 2024 - which has a faster algorithm, more memory transfer control and better readability

26 of 29

SEPTEMBER 2024

SU2 CONFERENCE

So where is this version?

27 of 29

SEPTEMBER 2024

SU2 CONFERENCE

It doesn’t compile...

28 of 29

SEPTEMBER 2024

SU2 CONFERENCE

Summary

    • GPU Acceleration seems very plausible in the SU2 code, the issue being frequent Memory Transfers required that will kill the performance of the code.

    • Both Space and Time integration are the major execution bottlenecks of the code and can be under the inspection for being accelerated. These subroutines also occur one after the another, so there is a possibility of minimizing data transfer.

    • Preconditioning and General Linear Algebra computations remain the most expensive protocols occurring under the hood, accelerating both holds good promise and currently seems to be possible.

29 of 29

SEPTEMBER 2024

SU2 CONFERENCE

See you (hopefully) for the next SU2 Conference