titlepage-color

GPU Accelerated Implicit Kinetic Meshfree Method based on Modified LU-SGS

Scientific Computing Report - 01:2024
Mayuri Verma, Nischay Ram Mamidi and Anil Nemili

###### Abstract

This report presents the GPU acceleration of implicit kinetic meshfree methods
using modified LU-SGS algorithms. The meshfree scheme is based on the least
squares kinetic upwind method (LSKUM). In the existing matrix-free LU-SGS
approaches for kinetic meshfree methods, the products of split flux Jacobians
and increments in conserved vectors are approximated by increments in the
split fluxes. In our modified LU-SGS approach, the Jacobian vector products
are computed exactly using algorithmic differentiation (AD). The implicit GPU
solvers with exact and approximate computation of the Jacobian vector products
are applied to the standard test cases for two-dimensional inviscid flows.
Numerical results have shown that the GPU solvers with the exact computation
of the Jacobian vector products are computationally more efficient and yield
better convergence rates than the solvers with approximations to the Jacobian
vector products. Benchmarks are presented to assess the performance of
implicit GPU solvers compared to the explicit GPU solver and the implicit
serial LSKUM solver.

##  1 Introduction

Numerical simulation of flow fields around complex geometries of specific
interest in aerodynamics is computationally intensive. Traditionally,
computational fluid dynamics (CFD) researchers used CPU parallel codes for
such simulations [1, 2]. However, with significant technological advancements,
GPUs have evolved from computer graphics application devices to truly general-
purpose parallel processors [3]. With superior floating-point operation rates
than multicore CPUs, GPUs are more suited for highly parallel numerical
computations. In recent years, many research groups have developed finite
volume or finite element based GPU parallel solvers, and some of them applied
their codes to large-scale applications in CFD [4, 5].

In contrast to the mesh-based methods, the meshfree schemes [6, 7, 8, 9, 10,
11, 12] provide an alternative approach for the numerical solution of Euler or
Navier Stokes equations that govern the fluid flow. These schemes operate on a
distribution of points, known as a point cloud. The point cloud can be
obtained from structured, unstructured, chimera grids or quadtree algorithms.
A meshfree scheme of particular interest to us is the Least Squares Kinetic
Upwind Method (LSKUM) [6, 13, 14]. Over the past two decades, this method has
been successfully used for computing flows around realistic configurations,
moving boundaries, flutter prediction, and other applications [15, 16, 17, 18,
19, 20].

It is well-known that the explicit solvers for large-scale problems require
more iterations to converge the numerical solutions. A convergence
acceleration technique that is widely used to reduce the computational time is
the lower upper symmetric Gauss-Seidel (LU-SGS) algorithm [21]. This algorithm
requires the evaluation of split flux Jacobians in the implicit operator,
which is computationally expensive. To mitigate these costs, a matrix-free
version of the LU-SGS scheme has been developed [22]. In this approach, the
products of split flux Jacobians with incremental conserved vectors are
approximated by incremental changes in the split flux vectors. Note that these
Jacobian vector products are part of the off-diagonal terms of the implicit
operator. Furthermore, the split flux Jacobians that appear in the diagonal
matrix of the implicit operator are approximated by a function of unsplit flux
Jacobians and their spectral radii [21]. Typically, in finite difference and
finite volume methods, this approximation transforms the block diagonal matrix
into a scalar diagonal matrix. However, when implemented in the meshfree
LSKUM, it may not lead to a scalar diagonal matrix due to the non-symmetric
nature of the split stencils [23]. In order to extend the LU-SGS algorithm to
kinetic meshfree methods, Anandhanarayanan et al. [23, 19] developed a
modified approximation for the split flux Jacobians. In another work, Manish
et al. [24] used the original approximation in [21] but moved the unsplit
Jacobians that result in the block diagonal matrix to the explicit operator.

This research presents the modified versions of meshfree LU-SGS schemes
proposed by Anandhanarayanan et al. and Manish et al. In the modified LU-SGS,
instead of approximating the Jacobian vector products, they are computed
exactly using algorithmic differentiation (AD) techniques [25]. In the
numerical results, we show that the LSKUM solvers with modified LU-SGS
approaches yield a faster residual convergence rate than the solvers with
approximate Jacobian vector products.

In the works mentioned earlier on kinetic meshfree schemes, numerical
simulations were performed using explicit or implicit CPU parallel solvers.
Although an effort has been made to develop explicit kinetic meshfree GPU
solvers [26], to our knowledge, research on implicit kinetic meshfree solvers
on GPUs is yet to be pursued. This research presents the development of
implicit LSKUM solvers on GPUs using original and modified LU-SGS approaches.
The programming model CUDA is used to develop the GPU solvers. Benchmark
simulations are performed to show that the implicit GPU solvers with AD are
computationally more efficient than the solvers without AD.

In the numerical simulation of fluid flows, the changes in the flow field
occur at the wall boundary points first and then propagate inside the whole
computational domain. Therefore, the wall points are updated first in the
forward and backward sweeps of LU-SGS, followed by interior and outer boundary
points. We analyse the effect of such an implementation on the convergence and
run-time of an implicit LSKUM GPU solver and compare it with looping over all
points in the cloud.

This report is organised as follows. Section 2 presents the basic theory of
the meshfree LSKUM. Section 3 presents implicit schemes for LSKUM based on
modified LU-SGS algorithms. Details on developing a GPU accelerated implicit
LSKUM solver are given in Section 4. In Section 5, residual convergence
histories are shown for the standard test cases for $2D$ inviscid flow.
Section 6 presents benchmarks to assess the computational efficiency of the
implicit GPU solvers. Finally, conclusions are drawn in Section 7.

##  2 Basic Theory of LSKUM

In this section, we present the basic theory of the Least Squares Kinetic
Upwind Method (LSKUM) with reference to $2D$ Euler equations,

$\frac{\partial\boldsymbol{U}}{\partial
t}+\frac{\partial\boldsymbol{Gx}}{\partial
x}+\frac{\partial\boldsymbol{Gy}}{\partial y}=0$ (1)

where $\boldsymbol{U}$ is the conserved vector, $\boldsymbol{Gx}$ and
$\boldsymbol{Gy}$ are the flux vectors along the coordinates $x$ and $y$,
respectively. These vectors are defined by

$\boldsymbol{U}=\begin{bmatrix}\rho\\\ \rho u_{1}\\\ \rho u_{2}\\\ \rho
e\end{bmatrix},\;\;\boldsymbol{Gx}=\begin{bmatrix}\rho u_{1}\\\ p+\rho
u_{1}^{2}\\\ \rho u_{1}u_{2}\\\ \left(p+\rho
e\right)u_{1}\end{bmatrix},\;\;\;\;\boldsymbol{Gy}=\begin{bmatrix}\rho
u_{2}\\\ \rho u_{1}u_{2}\\\ p+\rho u_{2}^{2}\\\ \left(p+\rho
e\right)u_{2}\end{bmatrix}$ (2)

Here, $\rho$ is the fluid density, $p$ is the pressure, $e$ is the internal
energy per unit mass, $u_{1}$ and $u_{2}$ are the fluid velocities along $x$
and $y$ coordinate directions, respectively. The above conservation laws can
be obtained as the moments of the $2D$ Boltzmann equation with the velocity
distribution function being Maxwellian,

$\left<\boldsymbol{\Psi},\frac{\partial F}{\partial t}+v_{1}\frac{\partial
F}{\partial x}+v_{2}\frac{\partial F}{\partial
y}\right>=0,\;\;\left<\boldsymbol{\Psi},v_{i}F\right>=\int\limits_{\mathcal{R^{+}}\times\mathcal{R}^{2}}v_{i}\boldsymbol{\Psi}Fdv_{1}dv_{2}dI,\;i=1,2$
(3)

where, $F$ is the Maxwellian velocity distribution function and
$\boldsymbol{\Psi}$ is the moment function vector, defined by

$F=\frac{\rho}{I_{0}}\frac{\beta}{\pi}e^{-\beta\left\\{\left(v_{1}-u_{1}\right)^{2}+\left(v_{2}-u_{2}\right)^{2}\right\\}-\frac{I}{I_{0}}},\;\;\;\boldsymbol{\Psi}=\left[1,v_{1},v_{2},I+\frac{v_{1}^{2}+v_{2}^{2}}{2}\right]$
(4)

Here, $v_{1}$ and $v_{2}$ are the molecular velocities along the coordinates
$x$ and $y$, respectively. $\beta=\rho/{2p}$, $I$ is the internal energy
variable, and $I_{0}$ is the internal energy due to non-translational degrees
of freedom.
Using the Courant-Issacson-Rees (CIR) splitting [27] of the molecular
velocities, an upwind scheme for the $2D$ Boltzmann equation in eq. (3) is
given by

$\frac{\partial F}{\partial
t}+\frac{v_{1}+\left|v_{1}\right|}{2}\frac{\partial F}{\partial
x}+\frac{v_{1}-\left|v_{1}\right|}{2}\frac{\partial F}{\partial
x}+\frac{v_{2}+\left|v_{2}\right|}{2}\frac{\partial F}{\partial
x}\frac{v_{2}-\left|v_{2}\right|}{2}\frac{\partial F}{\partial x}=0$ (5)

Taking $\boldsymbol{\Psi}$-moments of the above equation results in the upwind
kinetic flux vector split (KFVS) Euler equations [28],

$\frac{\partial\boldsymbol{U}}{\partial
t}+\frac{\partial\boldsymbol{Gx}^{+}}{\partial
x}+\frac{\partial\boldsymbol{Gx}^{-}}{\partial
x}+\frac{\partial\boldsymbol{Gy}^{+}}{\partial
y}+\frac{\partial\boldsymbol{Gy}^{-}}{\partial y}=0$ (6)

In order to find a meshfree upwind scheme for the above equation, the spatial
derivatives in eq. (5) are first discretised using the least-squares (LS)
principle. Suitable moments of the LS formulae with appropriate stencils lead
to the least squares kinetic upwind method (LSKUM). We first derive a second-
order accurate approximation for the spatial derivatives $F_{x}$ and $F_{y}$
at a point $P_{0}$ using the defect correction procedure [14]. Consider the
Taylor series expansion of $F$ up to quadratic terms at a point $P_{i}\in
N\left(P_{0}\right)$,

$\begin{split}\Delta F_{i}&=\Delta x_{i}F_{x_{0}}+\Delta
y_{i}F_{y_{0}}+\frac{\Delta x_{i}}{2}\left(\Delta x_{i}F_{xx_{0}}+\Delta
y_{i}F_{xy_{0}}\right)\\\ &+\frac{\Delta y_{i}}{2}\left(\Delta
x_{i}F_{xy_{0}}+\Delta y_{i}F_{yy_{0}}\right)+O\left(\Delta x_{i},\Delta
y_{i}\right)^{3},\;i=1,\ldots,n\end{split}$ (7)

where $\Delta x_{i}=x_{i}-x_{0}$, $\Delta y_{i}=y_{i}-y_{0}$, $\Delta
F_{i}=F_{i}-F_{0}$. $N\left(P_{0}\right)$ is the set of neighbours or the
stencil of $P_{0}$. Here, $n$ denotes the number of neighbours of the point
$P_{0}$. Replacing the second-order terms in the above equation using the
Taylor series expansions of $F_{x}$ and $F_{y}$ to linear terms at a point
$P_{i}$, we get

$\Delta F_{i}-\frac{1}{2}\left\\{\Delta x_{i}\Delta F_{x_{i}}+\Delta
y_{i}\Delta F_{y_{i}}\right\\}=\Delta x_{i}F_{x_{0}}+\Delta
y_{i}F_{y_{0}}++O\left(\Delta x_{i},\Delta y_{i}\right)^{3},\;\;i=1,\ldots,n$
(8)

Denoting the left hand side as
$\Delta\widetilde{F}_{i}=\widetilde{F}_{i}-\widetilde{F}_{0}$, the above
equation reduces to

$\Delta\widetilde{F}_{i}=\Delta x_{i}F_{x_{0}}+\Delta
y_{i}F_{y_{0}}+O\left(\Delta x_{i},\Delta y_{i}\right)^{3},\;i=1,\ldots,n$ (9)

For $n\geq 3$, eq. (9) leads to an over-determined linear system of equations.
Using the least-squares principle, the second-order approximations to $F_{x}$
and $F_{y}$ at the point $P_{0}$ are given by

$\begin{bmatrix}F_{x}\\\ F_{y}\end{bmatrix}_{P_{0}}={\begin{bmatrix}\sum\Delta
x_{i}^{2}&\sum\Delta x_{i}\Delta y_{i}\\\ \sum\Delta x_{i}\Delta
y_{i}&\sum\Delta y_{i}^{2}\end{bmatrix}}^{-1}\begin{bmatrix}\sum\Delta
x_{i}\Delta\widetilde{F}_{i}\\\ \sum\Delta y_{i}\Delta\widetilde{F}_{i}\\\
\end{bmatrix}_{P_{i}\in N\left(P_{0}\right)}$ (10)

Taking $\boldsymbol{\Psi}$ \- moments of eq. (5) along with the above least-
squares formulae, we get the semi-discrete second-order upwind meshfree scheme
based on LSKUM for $2D$ Euler equations. The least-squares formulae for the
spatial derivatives of $\boldsymbol{Gx}^{\pm}$ at a point $P_{0}$ are given by

$\frac{\partial\boldsymbol{Gx}^{\pm}}{\partial x}=\left[\frac{\sum\Delta
y_{i}^{2}\sum\Delta
x_{i}\Delta\widetilde{\boldsymbol{Gx}}_{i}^{\pm}-\sum\Delta x_{i}\Delta
y_{i}\sum\Delta y_{i}\Delta\widetilde{\boldsymbol{Gx}}_{i}^{\pm}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right]$ (11)

Here, the perturbations $\Delta\widetilde{\boldsymbol{Gx}}_{i}^{\pm}$ are
defined by

$\Delta\widetilde{\boldsymbol{Gx}}_{i}^{\pm}=\Delta{\boldsymbol{Gx}}^{\pm}_{i}-\frac{1}{2}\left\\{\Delta
x_{i}\frac{\partial}{\partial x}\Delta\boldsymbol{Gx}_{i}^{\pm}+\Delta
y_{i}\frac{\partial}{\partial y}\Delta\boldsymbol{Gx}_{i}^{\pm}\right\\}$ (12)

The split flux derivatives in eq. (11) are approximated using the stencils

$N_{x}^{\pm}\left(P_{0}\right)=\left\\{P_{i}\in N\left(P_{0}\right)\mid\Delta
x_{i}\lessgtr 0\right\\}$ (13)

Similarly, using appropriate split stencils, we can write the formulae for the
spatial derivatives of $\boldsymbol{Gy}^{\pm}$. A drawback of the defect
correction approach using the Maxwellian distributions is that the numerical
solution may not be positive as $\Delta\widetilde{F}_{i}$ is not the
difference between two Maxwellians. Instead, it is the difference between two
perturbed Maxwellian distributions $\widetilde{F}_{i}$ and $\widetilde{F}_{0}$
[14, 29]. For preserving the positivity of the solution, instead of the
Maxwellians, $\boldsymbol{q}$-variables [30] can be used in the defect
correction procedure [29]. Note that $\boldsymbol{q}$-variables can represent
the fluid flow at the macroscopic level as the transformations
$\boldsymbol{U}\longleftrightarrow\boldsymbol{q}$ and
$F\longleftrightarrow\boldsymbol{q}$ are unique. The
$\boldsymbol{q}$-variables in $2D$ are defined by

$\boldsymbol{q}=\begin{bmatrix}\ln\rho+\frac{\ln\beta}{\gamma-1}-\beta\left(u_{1}^{2}+u_{2}^{2}\right),\;2\beta
u_{1},\;2\beta u_{2},\;-2\beta\end{bmatrix},\;\beta=\frac{\rho}{2p}$ (14)

Where $\gamma$ is the ratio of the specific heats. Using
$\boldsymbol{q}$-variables, second-order accuracy can be achieved by replacing
$\Delta\widetilde{\boldsymbol{Gx}}_{i}^{\pm}$ in eq. (12) with
$\Delta\boldsymbol{Gx}_{i}^{\pm}\left(\widetilde{\boldsymbol{q}}\right)=\boldsymbol{Gx}^{\pm}\left(\widetilde{\boldsymbol{q}}_{i}\right)-\boldsymbol{Gx}^{\pm}\left(\widetilde{\boldsymbol{q}}_{0}\right)$.
Here, $\widetilde{\boldsymbol{q}}_{i}$ and $\widetilde{\boldsymbol{q}}_{0}$
are the modified $\boldsymbol{q}$-variables, defined by

$\widetilde{\boldsymbol{q}}_{i}=\boldsymbol{q}_{i}-\frac{1}{2}\left(\Delta
x_{i}{\boldsymbol{q}_{x}}_{i}+\Delta
y_{i}{\boldsymbol{q}_{y}}_{i}\right),\;\;\widetilde{\boldsymbol{q}}_{0}=\boldsymbol{q}_{0}-\frac{1}{2}\left(\Delta
x_{i}{\boldsymbol{q}_{x}}_{0}+\Delta y_{i}{\boldsymbol{q}_{y}}_{0}\right)$
(15)

Here $\boldsymbol{q}_{x}$ and $\boldsymbol{q}_{y}$ are evaluated to second-
order using the least-squares formuale

$\begin{bmatrix}\boldsymbol{q}_{x}\\\
\boldsymbol{q}_{y}\end{bmatrix}={\begin{bmatrix}\sum\Delta
x_{i}^{2}&\sum\Delta x_{i}\Delta y_{i}\\\ \sum\Delta x_{i}\Delta
y_{i}&\sum\Delta y_{i}^{2}\end{bmatrix}}^{-1}\begin{bmatrix}\sum\Delta
x_{i}\Delta\widetilde{\boldsymbol{q}}_{i}\\\ \sum\Delta
y_{i}\Delta\widetilde{\boldsymbol{q}}_{i}\\\ \end{bmatrix}_{P_{i}\in
N\left(P_{0}\right)}$ (16)

The above formulae for $\boldsymbol{q}_{x}$ and $\boldsymbol{q}_{y}$ are
implicit and must be solved iteratively. These iterations are known as inner
iterations.

##  3 Implicit LSKUM based on Modified LU-SGS

Consider the kinetic flux vector split $2D$ Euler equations in the implicit
form.

$\left(\frac{\partial\boldsymbol{U}}{\partial
t}\right)^{n}+\frac{\partial}{\partial
x}\left(\boldsymbol{Gx}^{+}\right)^{n+1}+\frac{\partial}{\partial
x}\left(\boldsymbol{Gx}^{-}\right)^{n+1}+\frac{\partial}{\partial
y}\left(\boldsymbol{Gy}^{+}\right)^{n+1}+\frac{\partial}{\partial
y}\left(\boldsymbol{Gy}^{-}\right)^{n+1}=0$ (17)

Linearising the split fluxes $\left(\boldsymbol{Gx}^{\pm}\right)^{n+1}$ and
$\left(\boldsymbol{Gy}^{\pm}\right)^{n+1}$ in time, we get

$\left(\boldsymbol{Gx}^{\pm}\right)^{n+1}=\left(\boldsymbol{Gx}^{\pm}\right)^{n}+\left(\boldsymbol{Ax}^{\pm}\right)^{n}\delta\boldsymbol{U}^{n},\;\;\left(\boldsymbol{Gy}^{\pm}\right)^{n+1}=\left(\boldsymbol{Gy}^{\pm}\right)^{n}+\left(\boldsymbol{Ay}^{\pm}\right)^{n}\delta\boldsymbol{U}^{n}$
(18)

where, $\delta\boldsymbol{U}^{n}=\boldsymbol{U}^{n+1}-\boldsymbol{U}^{n}$ is
the solution increment vector. $\boldsymbol{Ax}^{\pm}$ and
$\boldsymbol{Ay}^{\pm}$ are the split flux Jacobians along the coordinate
directions $x$ and $y$ respectively. Using the above expressions in eq. (17),
we obtain

$\begin{split}\left(\frac{\partial\boldsymbol{U}}{\partial
t}\right)^{n}+&\frac{\partial}{\partial
x}\left(\boldsymbol{Ax}^{+}\delta\boldsymbol{U}\right)^{n}+\frac{\partial}{\partial
x}\left(\boldsymbol{Ax}^{-}\delta\boldsymbol{U}\right)^{n}+\frac{\partial}{\partial
y}\left(\boldsymbol{Ay}^{+}\delta\boldsymbol{U}\right)^{n}+\frac{\partial}{\partial
y}\left(\boldsymbol{Ay}^{-}\delta\boldsymbol{U}\right)^{n}\\\
=&-\left(\frac{\partial\boldsymbol{Gx}^{+}}{\partial
x}+\frac{\partial\boldsymbol{Gx}^{-}}{\partial
x}+\frac{\partial\boldsymbol{Gy}^{+}}{\partial
y}+\frac{\partial\boldsymbol{Gy}^{-}}{\partial y}\right)^{n}\\\
=&-\boldsymbol{R}^{n}\end{split}$ (19)

Here, $\boldsymbol{R}^{n}$ is the second-order accurate flux residual based on
meshfree LSKUM. In the above equation, the time derivative is approximated
using a first-order forward difference formula. The spatial derivatives on the
left-hand side are discretised using the first-order least-squares
approximations with appropriate split stencils.

To simplify the left-hand side of the above equation, consider the least-
squares formulae for the spatial derivatives of
$\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}$ at the point
$P_{0}$

$\frac{\partial}{\partial
x}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}=\left\\{\frac{\sum\Delta
y_{i}^{2}\sum\Delta
x_{i}\Delta\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\sum\Delta
x_{i}\Delta y_{i}\sum\Delta
y_{i}\Delta\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}$ (20)

Using the expression
$\Delta\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}=\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}$,
the above equation reduces to

$\begin{split}\frac{\partial}{\partial
x}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}&=\left\\{\frac{\sum\Delta
y_{i}^{2}\sum\Delta
x_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\sum\Delta
x_{i}\Delta y_{i}\sum\Delta
y_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\ &-\left\\{\frac{\sum\Delta y_{i}^{2}\sum\Delta
x_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}-\sum\Delta
x_{i}\Delta y_{i}\sum\Delta
y_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\end{split}$ (21)

Regrouping the first term on the right-hand side of the above equation based
on $P_{i}<P_{0}$ and $P_{i}>P_{0}$, we arrive at

$\begin{split}\frac{\partial}{\partial
x}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}&=\left\\{\frac{\sum\Delta
y_{i}^{2}\sum\limits_{i<0}\Delta
x_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\sum\Delta
x_{i}\Delta y_{i}\sum\limits_{i<0}\Delta
y_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\ &+\left\\{\frac{\sum\Delta
y_{i}^{2}\sum\limits_{i>0}\Delta
x_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\sum\Delta
x_{i}\Delta y_{i}\sum\limits_{i>0}\Delta
y_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\ &-\left\\{\frac{\sum\Delta y_{i}^{2}\sum\Delta
x_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}-\sum\Delta
x_{i}\Delta y_{i}\sum\Delta
y_{i}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\end{split}$ (22)

Denoting the first, second and third terms on the right-hand side of the above
equation as
$LSx\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i,i<0}$,
$LSx\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i,i>0}$ and
$LSx\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}$
respectively, we get

$\frac{\partial}{\partial
x}\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}=LSx\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i,i<0}+LSx\left(\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i,i>0}-LSx\left(\boldsymbol{Ax}^{\pm}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}$
(23)

Similarly, the spatial approximations of
$\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}$, evaluated
using the split stencils $Ny^{\pm}_{0}=\left\\{P_{i}\in
N\left(P_{0}\right)\mid\Delta y_{i}\lessgtr 0\right\\}$ are given by

$\begin{split}\frac{\partial}{\partial
y}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}&=\left\\{\frac{\sum\Delta
x_{i}^{2}\sum\limits_{i<0}\Delta
y_{i}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\sum\Delta
x_{i}\Delta y_{i}\sum\limits_{i<0}\Delta
x_{i}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\ &+\left\\{\frac{\sum\Delta
x_{i}^{2}\sum\limits_{i>0}\Delta
y_{i}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}-\sum\Delta
x_{i}\Delta y_{i}\sum\limits_{i>0}\Delta
x_{i}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\ &-\left\\{\frac{\sum\Delta x_{i}^{2}\sum\Delta
y_{i}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}-\sum\Delta
x_{i}\Delta y_{i}\sum\Delta
x_{i}\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{0}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\
&=LSy\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i,i<0}+LSy\left(\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}\right)^{n}_{i,i>0}-LSy\left(\boldsymbol{Ay}^{\pm}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}\end{split}$
(24)

Using the expressions obtained in eqs. (23) and (24), eq. (19) reduces to

$\begin{split}&\left\\{\frac{\mathcal{I}}{{\Delta
t}_{0}}-LSx\left(\boldsymbol{Ax}^{+}\right)^{n}_{0}-LSx\left(\boldsymbol{Ax}^{-}\right)^{n}_{0}-LSy\left(\boldsymbol{Ay}^{+}\right)^{n}_{0}-LSy\left(\boldsymbol{Ay}^{-}\right)^{n}_{0}\right\\}\delta\boldsymbol{U}^{n}_{0}\\\
+&\left\\{LSx\left(\boldsymbol{Ax}^{+}\delta\boldsymbol{U}\right)^{n}_{i,i<0}+LSx\left(\boldsymbol{Ax}^{-}\delta\boldsymbol{U}\right)^{n}_{i,i<0}+LSy\left(\boldsymbol{Ay}^{+}\delta\boldsymbol{U}\right)^{n}_{i,i<0}+LSy\left(\boldsymbol{Ay}^{-}\delta\boldsymbol{U}\right)^{n}_{i,i<0}\right\\}\\\
+&\left\\{LSx\left(\boldsymbol{Ax}^{+}\delta\boldsymbol{U}\right)^{n}_{i,i>0}+LSx\left(\boldsymbol{Ax}^{-}\delta\boldsymbol{U}\right)^{n}_{i,i>0}+LSy\left(\boldsymbol{Ay}^{+}\delta\boldsymbol{U}\right)^{n}_{i,i>0}+LSy\left(\boldsymbol{Ay}^{-}\delta\boldsymbol{U}\right)^{n}_{i,i>0}\right\\}\\\
=&-\boldsymbol{R}^{n}_{0}\end{split}$ (25)

Denote the first, second, and third terms on the left-hand side of the above
equation as $\widetilde{\mathcal{D}}_{0}$, $\mathcal{L}_{0}$ and
$\mathcal{U}_{0}$ respectively. Here, $\widetilde{\mathcal{D}}$ is the block
diagonal matrix, $\mathcal{L}$ and $\mathcal{U}$ are strictly the lower and
upper triangular matrices. For all points in the given computational domain,
eq. (25) leads to a linear system of equations of the form

$\left(\mathcal{L}^{n}_{i}+\widetilde{\mathcal{D}}^{n}_{i}+\mathcal{U}^{n}_{i}\right)\delta\boldsymbol{U}^{n}_{i}=-\boldsymbol{R}^{n}_{i},\;\;i=1\ldots
N$ (26)

Here $N$ is the number of points in the domain. Using the implicit lower upper
symmetric Gauss-Seidel (LU-SGS) algorithm [21], the linear system can be
solved using a two-step procedure as

$\begin{split}\delta\boldsymbol{U}^{*}_{i}&=-{\widetilde{\mathcal{D}}}^{-1}_{i}\left(\boldsymbol{R}^{n}_{i}+\mathcal{L}_{i}\delta\boldsymbol{U}^{*}_{i}\right),\;\;i=1\ldots
N\;\;\;\text{(Forward sweep)}\\\
\delta\boldsymbol{U}^{n}_{i}&=\delta\boldsymbol{U}^{*}_{i}-\widetilde{\mathcal{D}}^{-1}_{i}\left(\mathcal{U}\delta\boldsymbol{U}^{n}\right),\;\;i=N\ldots
1\;\;\;\text{(Reverse sweep)}\end{split}$ (27)

We observe that the computation of $\mathcal{L}\delta\boldsymbol{U}^{*}$ and
$\mathcal{U}\delta\boldsymbol{U}^{n}$ in the forward and reverse sweeps
involves the evaluation of split flux Jacobian matrices and their products
with the solution increment vectors,
$\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}$ and
$\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}$. The computation of expensive
Jacobians can be avoided by approximating the Jacobian vector products with
the increments in the numerical flux functions,

$\begin{split}\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}&=\boldsymbol{Gx}^{\pm}\left(\boldsymbol{U}+\delta\boldsymbol{U}\right)-\boldsymbol{Gx}^{\pm}\left(\boldsymbol{U}\right)=\delta\boldsymbol{Gx}^{\pm}\\\
\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}&=\boldsymbol{Gy}^{\pm}\left(\boldsymbol{U}+\delta\boldsymbol{U}\right)-\boldsymbol{Gy}^{\pm}\left(\boldsymbol{U}\right)=\delta\boldsymbol{Gy}^{\pm}\end{split}$
(28)

This matrix-free approach proposed for finite volume methods [22] is followed
in the earlier works on implicit meshfree solvers based on LSKUM [23, 24, 19].
The LU-SGS algorithm also requires the inverse of split flux Jacobians to
compute $\widetilde{\mathcal{D}}^{-1}$. In order to ease the computation, the
split flux Jacobians are approximated using spectral radii [21] as

$\boldsymbol{Ax}^{\pm}=\frac{1}{2}\left(\boldsymbol{Ax}\pm\rho\left(\boldsymbol{Ax}\right)\boldsymbol{I}\right),\;\;\boldsymbol{Ay}^{\pm}=\frac{1}{2}\left(\boldsymbol{Ay}\pm\rho\left(\boldsymbol{Ay}\right)\boldsymbol{I}\right)$
(29)

Here, $\boldsymbol{Ax}$ and $\boldsymbol{Ay}$ are the unsplit flux Jacobians.
$\rho\left(\boldsymbol{Ax}\right)$ and $\rho\left(\boldsymbol{Ay}\right)$ are
the spectral radii of the flux Jacobians, given by

$\rho\left(\boldsymbol{Ax}\right)=\left|u_{1}\right|+a,\;\;\rho\left(\boldsymbol{Ay}\right)=\left|u_{2}\right|+a$
(30)

where $a$ is the speed of sound. Using these approximations the expression for
$LSx\left(\boldsymbol{Ax^{+}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}$ in
eq. (25) reduces to

$\begin{split}LSx\left(\boldsymbol{Ax^{+}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}&=\frac{1}{2}\left\\{\frac{\sum\Delta
y_{i}^{2}\sum\Delta
x_{i}\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}-\sum\Delta
x_{i}\Delta y_{i}\sum\Delta
y_{i}\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}}{\sum\Delta
x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\\\
&+\frac{1}{2}\rho\left(\boldsymbol{Ax}\right)\left\\{\frac{\sum\Delta
y_{i}^{2}\sum\Delta x_{i}-\sum\Delta x_{i}\Delta y_{i}\sum\Delta
y_{i}}{\sum\Delta x_{i}^{2}\sum\Delta y_{i}^{2}-\left(\sum\Delta x_{i}\Delta
y_{i}\right)^{2}}\right\\}\delta\boldsymbol{U}^{n}_{0}\\\
&=\frac{1}{2}LSxp\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}+\frac{1}{2}\rho\left(\boldsymbol{Ax}\right)LSxp(1)\delta\boldsymbol{U}^{n}_{0}\end{split}$
(31)

Similarly, the expressions for
$LSx\left(\boldsymbol{Ax^{-}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}$, and
$LSy\left(\boldsymbol{Ay^{\pm}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}$
can be written as

$\begin{split}LSx\left(\boldsymbol{Ax^{-}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}&=\frac{1}{2}LSxn\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}-\frac{1}{2}\rho\left(\boldsymbol{Ax}\right)LSxn(1)\delta\boldsymbol{U}^{n}_{0}\\\
LSy\left(\boldsymbol{Ay^{+}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}&=\frac{1}{2}LSyp\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}+\frac{1}{2}\rho\left(\boldsymbol{Ay}\right)LSyp(1)\delta\boldsymbol{U}^{n}_{0}\\\
LSy\left(\boldsymbol{Ay^{-}}\right)^{n}_{0}\delta\boldsymbol{U}^{n}_{0}&=\frac{1}{2}LSyn\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}-\frac{1}{2}\rho\left(\boldsymbol{Ay}\right)LSyn(1)\delta\boldsymbol{U}^{n}_{0}\\\
\end{split}$ (32)

From eqs. (31) and (32), $\widetilde{\mathcal{D}}\delta\boldsymbol{U}^{n}_{0}$
in eq. (25) reduces to

$\begin{split}\widetilde{\mathcal{D}}\delta\boldsymbol{U}^{n}_{0}&=-\frac{1}{2}\left\\{LSxp\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}+LSxn\left(\boldsymbol{Ax}\delta\boldsymbol{U}\right)^{n}_{0}+LSyp\left(\boldsymbol{Ay}\delta\boldsymbol{U}\right)^{n}_{0}+LSyn\left(\boldsymbol{Ay}\delta\boldsymbol{U}\right)^{n}_{0}\right\\}\\\
&+\left\\{\frac{\mathcal{I}}{\Delta
t}-\frac{1}{2}\rho\left(\boldsymbol{Ax}\right)\left\\{LSxp(1)-LSxn(1)\right\\}-\frac{1}{2}\rho\left(\boldsymbol{Ay}\right)\left\\{LSyp(1)-LSyn(1)\right\\}\right\\}\delta\boldsymbol{U}^{n}_{0}\\\
&=\mathcal{S}_{0}^{n}+\mathcal{D}\delta\boldsymbol{U}^{n}_{0}\end{split}$ (33)

In the above equation, $\mathcal{D}$ is a scalar diagonal matrix. The term
$\mathcal{S}$ on the right-hand side vanishes when the point distribution
around the point $P_{0}$ is symmetric, thus transforming the block diagonal
matrix $\widetilde{\mathcal{D}}$ to a scalar matrix. However, for an arbitrary
distribution of points, ensuring a symmetric stencil for each point may not
always be possible. To avoid the term $\mathcal{S}$, instead of eq. (29),
Anandhanarayanan et al. [23] approximated the split flux Jacobians using the
spectral radii,

$\boldsymbol{Ax}^{\pm}=\pm\rho\left(\boldsymbol{Ax}^{\pm}\right)\boldsymbol{I},\;\;\boldsymbol{Ay}^{\pm}=\pm\rho\left(\boldsymbol{Ay}^{\pm}\right)\boldsymbol{I}$
(34)

where $\rho\left(\boldsymbol{Ax}^{\pm}\right)$ and
$\rho\left(\boldsymbol{Ay}^{\pm}\right)$ are the spectral radii of
$\boldsymbol{Ax}^{\pm}$ and $\boldsymbol{Ay}^{\pm}$ respectively. They are
defined by

$\rho\left(\boldsymbol{Ax}^{\pm}\right)=\frac{1}{2}\left|\left(u_{1}\pm
a\right)\pm\left|u_{1}\pm
a\right|\right|,\;\;\rho\left(\boldsymbol{Ay}^{\pm}\right)=\frac{1}{2}\left|\left(u_{2}\pm
a\right)\pm\left|u_{2}\pm a\right|\right|$ (35)

In a later research work, Manish et al. [24] retained the term $\mathcal{S}$
and moved it to the explicit part of eq. (25). It was done to increase the
diagonal dominance of the linear system of equations and thus increase the
convergence rate. The linear system in eq. (26) then reduces to

$\left(\mathcal{L}^{n}_{i}+\mathcal{D}^{n}_{i}+\mathcal{U}^{n}_{i}\right)\delta\boldsymbol{U}^{n}_{i}=-\boldsymbol{R}^{n}_{i}+\mathcal{S}^{n}_{i},\;\;i=1,\ldots,N$
(36)

The LU-SGS algorithm is then given by

$\begin{split}\delta\boldsymbol{U}^{*}_{i}&=-{\mathcal{D}}^{-1}_{i}\left(\boldsymbol{R}^{n}_{i}-\mathcal{S}^{n}_{i}+\mathcal{L}_{i}\delta\boldsymbol{U}^{*}_{i}\right),\;\;i=1,\ldots,N\;\;\;\text{(Forward
sweep)}\\\
\delta\boldsymbol{U}^{n}_{i}&=\delta\boldsymbol{U}^{*}_{i}-\mathcal{D}^{-1}_{i}\left(\mathcal{U}\delta\boldsymbol{U}^{n}\right),\;\;i=N,\ldots,1\;\;\;\text{(Reverse
sweep)}\end{split}$ (37)

Since $\delta\boldsymbol{U}^{n}$ is the unknown vector, the term $\mathcal{S}$
can be computed using the available flow solution at the $(n-1)^{th}$
iteration. The Jacobian vector products $\boldsymbol{Ax}\delta\boldsymbol{U}$
and $\boldsymbol{Ay}\delta\boldsymbol{U}$ present in $\mathcal{S}$ are
evaluated using the matrix-free approximations,

$\begin{split}\boldsymbol{Ax}\delta\boldsymbol{U}&=\boldsymbol{Gx}\left(\boldsymbol{U}+\delta\boldsymbol{U}\right)-\boldsymbol{Gx}\left(\boldsymbol{U}\right)=\delta\boldsymbol{Gx}\\\
\boldsymbol{Ay}\delta\boldsymbol{U}&=\boldsymbol{Gy}\left(\boldsymbol{U}+\delta\boldsymbol{U}\right)-\boldsymbol{Gy}\left(\boldsymbol{U}\right)=\delta\boldsymbol{Gy}\end{split}$
(38)

In the numerical results, we show that the implicit LSKUM solver with the LU-
SGS algorithm in eq. (37) does not improve the residual convergence rate over
the LU-SGS scheme in eq. (27).

In this research, the Jacobian vector products present in the term
$\mathcal{S}$ of eq. (37), and also in the computation of
$\mathcal{L}\delta\boldsymbol{U}^{*}$ and
$\mathcal{U}\delta\boldsymbol{U}^{n}$ in eqs. (27) and (37) are evaluated
exactly using the algorithmic differentiation (AD) techniques [25]. In the
numerical results, we show that the exact computation of the Jacobian vector
products notably enhances the convergence rate.

We illustrate the exact computation of the flux Jacobian vector product
$\boldsymbol{Ax}\delta\boldsymbol{U}$ using AD. Listing 1 shows a Fortran
subroutine, $\tt flux\\_Gx$ for computing the full flux vector
$\boldsymbol{Gx}$ along the coordinate direction $x$. Given the conserved
vector $\boldsymbol{U}$ as the input, this routine yields the flux vector
$\boldsymbol{Gx}$ as the output. Listing 2 shows the tangent linear code
generated by algorithmically differentiating the routine $\tt flux\\_Gx$ in a
black-box fashion using the AD tool Tapenade [31]. Here, the derivative
routine is named by appending $\tt\\_d$ to the primal routine $\tt flux\\_Gx$.
By specifying the input vectors $\tt U$ and $\tt Ud$, the tangent code
computes the flux vector $\tt Gx$ and its derivative $\tt Gxd$ in the
direction of $\tt Ud$. Here, $\tt Gxd$ is the Jacobian vector product
$\boldsymbol{Ax}\delta\boldsymbol{U}$. Since we are interested in finding $\tt
Gxd$ alone, the black-box tangent code can be modified to omit the computation
of the flux vector $\tt Gx$. Listing 3 shows the optimised version of the
tangent code. Given the input vectors $\tt U$ and $\tt Ud$, this code gives
the desired Jacobian vector product $\tt AxdU$ as the output. Following the
same approach, the split flux Jacobian vector products,
$\boldsymbol{Ax}^{\pm}\delta\boldsymbol{U}$ and
$\boldsymbol{Ay}^{\pm}\delta\boldsymbol{U}$ too are computed exactly using AD.
Note that the computation of a split flux Jacobian conserved vector product
using incremental fluxes requires two function calls, one with
$\boldsymbol{U}+\delta\boldsymbol{U}$ and the other with $\boldsymbol{U}$ as
input vectors. Since the expressions for kinetic split fluxes involve an error
function, invoking two split flux functions doubled the calls to the error
function. On the other hand, computing a split flux Jacobian vector product
using AD involves only one call to the error function. In Section 6, we show
that the implicit meshfree GPU solvers using AD are computationally more
efficient than those employing matrix-free approximations for the flux
Jacobian conserved vector products.

##  4 GPU Accelerated Implicit LSKUM Solver

This section presents the development of a GPU accelerated implicit meshfree
LSKUM solver based on LU-SGS. Algorithm 1 presents the general structure of a
serial implicit solver for steady-state flows. The solver consists of a fixed
point iterative scheme, represented by the $\tt for$ loop. The subroutine $\tt
q\\_variables()$ computes the $\boldsymbol{q}$-variables defined in eq. (14).
$\tt q\\_derivatives()$ performs the inner iterations in the defect correction
step to evaluate $\boldsymbol{q}_{x}$ and $\boldsymbol{q}_{y}$ using the
formulae in eq. (16). $\tt flux\\_residual()$ computes the second-order
accurate approximation of the kinetic split flux derivatives in eq. (19) using
least-squares formulae. The routine $\tt LUSGS()$ performs the forward and
backward sweeps of the implicit scheme. Since the solution develops at the
wall boundary, the wall points are first updated in the forward and backward
sweeps, followed by interior and outer boundary points. Algorithm 2 presents
the corresponding forward and backward sweep subroutines. Finally, $\tt
state\\_update()$ updates the flow solution. All input and output operations
are performed in the subroutines $\tt preprocessor()$ and $\tt
postprocessor()$, respectively. The parameter $N$ represents the fixed point
iterations required to achieve the desired convergence in the flow solution.

subroutine  _LSKUM()_ call  preprocessor()  for  _$n\leftarrow 1$ to $n\leq
N$_ do call q_variables()  call q_derivatives()  call flux_residual()  call
LUSGS()  call state_update()  end for call  postprocessor()  end subroutine
Algorithm 1 Implicit LSKUM.

subroutine  _LUSGS_forward_sweep()_

for  _$i\leftarrow 1$ to _wall_points_ _ do

$p={\tt wall\\_point(i)}$

call forward_sweep(p)

end for

for  _$i\leftarrow 1$ to _interior_points_ _ do

$p={\tt interior\\_point(i)}$

call forward_sweep(p)

end for

for  _$i\leftarrow 1$ to _outer_points_ _ do

$p={\tt outer\\_point(i)}$

call forward_sweep(p)

end for

end subroutine

subroutine  _LUSGS_backward_sweep()_

for  _$i\leftarrow{\tt wall\\_points}$ to $i\geq 1$ _ do

$p={\tt wall\\_point(i)}$

call backward_sweep(p)

end for

for  _$i\leftarrow{\tt interior\\_points}$ to $i\geq 1$ _ do

$p={\tt interior\\_point(i)}$

call backward_sweep(p)

end for

for  _$i\leftarrow{\tt outer\\_points}$ to $i\geq 1$ _ do

$p={\tt outer\\_point(i)}$

call backward_sweep(p)

end for

end subroutine

Algorithm 2 Forward and Backward sweeps of the implicit LU-SGS algorithm.

Algorithm 3 presents the structure of a GPU accelerated implicit LSKUM solver.
The GPU solver mainly consists of the following sequence of operations:
transfer the input meshfree data structure from host to device, perform the
fixed-point iterations on the GPU, and finally transfer the converged flow
solution from device to host. In the baseline GPU code, for each subroutine in
Algorithm 1, equivalent CUDA kernels are constructed. The GPU code is then
profiled using NVIDIA Nsight Compute [32]. The profiler metrics have shown
that the register usage of the flux_residual and LUSGS kernels is very high as
these kernels are too large. When the kernel requires more registers than what
is available on the GPU, it causes register spilling. This results in
performance degradation as the excess data is stored in the slower global
memory, leading to an increase in memory latency [26]. In order to reduce the
register pressure, these kernels are split into smaller kernels. The
flux_residual kernel is split into four kernels that compute the spatial
derivatives of the split fluxes $\boldsymbol{Gx}^{\pm}$ and
$\boldsymbol{Gy}^{\pm}$. From eq. (19), the sum of these spatial derivatives
gives the flux residual $\boldsymbol{R}^{n}$. Similarly, the LUSGS kernel has
been split into two kernels that perform the forward and backward sweeps of
the implicit scheme. These changes have reduced the register pressure and
improved the performance of the GPU solver.

From Algorithm 2, we observe that the i-loops in the forward and backward
sweeps must be executed sequentially. Numerical experiments have shown that
the implicit GPU solver based on Algorithm 2 yields a marginal speedup of
around three compared to the serial solver. This is due to the inherent
structure of the LU-SGS algorithm, which prevents the full utilisation of GPU
resources. On the other hand, launching a CUDA kernel by ignoring the loop
structure of the forward sweep would update $\delta\boldsymbol{U}^{*}$ at a
point $P_{i}$. However, it could lead to an undesirable scenario where threads
operating simultaneously on other points might access
$\delta\boldsymbol{U}^{*}_{i}$, which has not been updated. This race
condition is caused by the presence of the point $P_{i}$ in the connectivity
set of other points.

subroutine  _LSKUM_cuda()_ call  preprocessor()  cudaHostToDevice(CPU_data,
GPU_data)  for  _$n\leftarrow 1$ to $n\leq N$_ do kernel $<<<$ grid, block
$>>>$ q_variables_cuda()  kernel $<<<$ grid, block $>>>$ q_derivatives_cuda()
kernel $<<<$ grid, block $>>>$ spatial_derivative_Gx_positive_cuda()  kernel
$<<<$ grid, block $>>>$ spatial_derivative_Gx_negative_cuda()  kernel $<<<$
grid, block $>>>$ spatial_derivative_Gy_positive_cuda()  kernel $<<<$ grid,
block $>>>$ spatial_derivative_Gy_negative_cuda()  kernel $<<<$ grid, block
$>>>$ LUSGS_forward_sweep_cuda()  kernel $<<<$ grid, block $>>>$
LUSGS_backward_sweep_cuda()  kernel $<<<$ grid, block $>>>$
state_update_cuda()  end for cudaDeviceToHost(GPU_data, CPU_data)  call
postprocessor()  end subroutine Algorithm 3 GPU accelerated implicit LSKUM
solver.

To address this challenge and to extract parallelism, modifications must be
made to the forward and reverse sweep algorithms without affecting the
connectivity related data dependency. It can be achieved by dividing the
points in the computational domain into several disjoint sets. Here, each set
consists of points that do not lie in the connectivity of other points. To
distinguish these sets, a unique color is assigned to the points in each set.
Typically, the point colors are represented by positive integers. An advantage
of this approach is that the forward and reverse sweeps of the points with the
same color can be executed in parallel using CUDA threads.

The coloring scheme employed in the present work is shown in Algorithm 4. The
basic idea of this scheme is that a point and its neighbours should have
different colors. A necessary condition for the coloring scheme is that if a
point $P_{i}$ is in the neighbourhood of $P_{0}$, then
$C\left(P_{i}\right)\cap
C\left(P_{0}\right)\supseteq\left\\{P_{i},P_{0}\right\\}$. Here,
$C\left(P_{i}\right)$ and $C\left(P_{0}\right)$ are the connectivity sets of
the points $P_{i}$ and $P_{0}$, respectively. Violating this criterion may
result in a point and some of its neighbours having the same color. However,
this condition on the connectivity sets is not required for the least-squares
approximation of the spatial derivatives in the meshfree numerical algorithm.

subroutine  _point_coloring_scheme()_

Comment: Initialise the colors of all points.

for  _$i\leftarrow 1$ to ${\tt max\\_points}$_ do

$color(i)=0$

end for

$color(1)=1$

Comment: Find the color of subsequent points.

for  _$i\leftarrow 1$ to $\tt max\\_points$_ do

for  _$j\leftarrow 1$ to $\tt neighbours(i)$_ do

$p$ = connectivity(i,j)

if  _$\left(color(p)=0\right)$_ then

$S=\text{set of colors of the neighbours of }p$

$color(p)=\text{smallest integer }k>0{\text{ and }}k\notin S$

end if

end for

end for

end subroutine

Algorithm 4 Coloring scheme for the point distribution.

Algorithm 5 shows the modified forward and backward sweeps using the point
colors. Here, the i-loops, which run over the point colors, must be executed
sequentially. On the other hand, the j-loops correspond to the points with the
same color. For these points, the forward sweeps can be performed
simultaneously on a GPU using CUDA threads. The more the number of points in a
j-loop, the higher the number of CUDA threads executed parallelly, leading to
an increase in the speedup. The performance of the GPU solver can be enhanced
further if only one i-loop exists in the forward sweep. Note that the three
i-loops in this algorithm are due to the grouping of the point distribution
into wall, interior, and outer boundary points.

Algorithm 6 shows a computationally more efficient implementation of forward
and backward sweeps. Here, the given point distribution is treated as one
group, which results in only one i-loop per sweep. Since the number of points
with a given color is more than the number of wall, interior, or outer
boundary points with the same color, it allows for more thread parallelism on
a GPU. Although this approach may result in a slower convergence rate, it is
well compensated by a significant speedup over Algorithm 5. The corresponding
CUDA kernels are shown in Algorithm 7. The computational efficiency of this
approach over Algorithm 5 is demonstrated in the numerical results.

subroutine  _LUSGS_forward_sweep()_

for  _$i\leftarrow 1$ to _wall_point_colors_ _ do

$n_{w}={\tt wall\\_points\\_with\\_color(i)}$

if  _$\left(n_{w}\neq 0\right)$_ then

for  _$j\leftarrow 1$ to $n_{w}$ _ do

$p$ = wall_point(i,j)

call forward_sweep(p)

end for

end if

end for

for  _$i\leftarrow 1$ to _interior_point_colors_ _ do

$n_{i}={\tt interior\\_points\\_with\\_color(i)}$

if  _$\left(n_{i}\neq 0\right)$_ then

for  _$j\leftarrow 1$ to $n_{i}$ _ do

$p$ = interior_point(i,j)

call forward_sweep(p)

end for

end if

end for

for  _$i\leftarrow 1$ to _outer_point_colors_ _ do

$n_{o}={\tt outer\\_points\\_with\\_color(i)}$

if  _$\left(n_{o}\neq 0\right)$_ then

for  _$j\leftarrow 1$ to $n_{o}$ _ do

$p$ = outer_point(i,j)

call forward_sweep(p)

end for

end if

end for

end subroutine

subroutine  _LUSGS_backward_sweep()_

for  _$i\leftarrow$ wall_point_colors to $i\geq 1$ _ do

$n_{w}={\tt wall\\_points\\_with\\_color(i)}$

if  _$\left(n_{w}\neq 0\right)$_ then

for  _$j\leftarrow 1$ to $n_{w}$ _ do

$p$ = wall_point(i,j)

call backward_sweep(p)

end for

end if

end for

for  _$i\leftarrow$ interior_point_colors to $i\geq 1$_ do

$n_{i}={\tt interior\\_points\\_with\\_color(i)}$

if  _$\left(n_{i}\neq 0\right)$_ then

for  _$j\leftarrow 1$ to $n_{i}$ _ do

$p$ = interior_point(i,j)

call backward_sweep(p)

end for

end if

end for

for  _$i\leftarrow$ outer_point_colors to $i\geq 1$_ do

$n_{o}={\tt outer\\_points\\_with\\_color(i)}$

if  _$\left(n_{o}\neq 0\right)$_ then

for  _$j\leftarrow 1$ to $n_{o}$ _ do

$p$ = outer_point(i,j)

call backward_sweep(p)

end for

end if

end for

end subroutine

Algorithm 5 Modified forward and backward sweeps.

subroutine  _LUSGS_forward_sweep()_

for  _$i\leftarrow 1$ to _point_colors_ _ do

$n_{p}={\tt points\\_with\\_color(i)}$

for  _$j\leftarrow 1$ to $n_{p}$ _ do

$p$ = point(i,j)

call forward_sweep(p)

end for

end for

end subroutine

subroutine  _LUSGS_backward_sweep()_

for  _$i\leftarrow{\tt point\\_colors}$ to $i\geq 1$ _ do

$n_{p}={\tt points\\_with\\_color(i)}$

for  _$j\leftarrow 1$ to $n_{p}$ _ do

$p$ = point(i,j)

call backward_sweep(p)

end for

end for

end subroutine

Algorithm 6 Efficient forward and backward sweeps.

subroutine  _LUSGS_forward_sweep_cuda()_

for  _$i\leftarrow 1$ to _point_colors__ do

$n_{p}={\tt points\\_with\\_color(i)}$

kernel $<<<$ grid, block $>>>$ forward_sweep_cuda($n_{p}$)

end for

end subroutine

subroutine  _LUSGS_backward_sweep_cuda()_

for  _$i\leftarrow{\tt point\\_colors}$ to $1$ _ do

$n_{p}={\tt points\\_with\\_color(i)}$

kernel $<<<$ grid, block $>>>$ backward_sweep_cuda($n_{p}$)

end for

end subroutine

Algorithm 7 CUDA kernels for efficient forward and backward sweeps.

##  5 Numerical Results

In this section, we illustrate the performance of the GPU accelerated implicit
meshfree LSKUM solvers with exact computation of the flux Jacobian conserved
vector products in LU-SGS algorithms presented in Section 3. The test cases
considered are the inviscid subsonic, transonic, and supersonic flows over the
NACA 0012 airfoil and the subsonic flow over a three-element airfoil. For
residual comparisons, the following GPU solvers are considered.

  1. (a)

Implicit GPU solver based on the LU-SGS scheme in eq. (37) with incremental
flux approximations for the flux Jacobian conserved vector products. We refer
to it as the implicit solver based on Manish et al.’s approach.

  2. (b)

Implicit GPU solver based on the LU-SGS scheme in eq. (37) with the exact
computation of the flux Jacobian conserved vector products. We denote this
approach as the implicit solver based on Manish et al.+AD.

  3. (c)

Implicit GPU solver based on the LU-SGS scheme in eq. (27) with incremental
flux approximations for the flux Jacobian conserved vector products. We refer
to it as the implicit solver based on Anandhanarayanan et al.’s approach.

  4. (d)

Implicit GPU solver based on the LU-SGS scheme in eq. (27) with the exact
computation of the flux Jacobian conserved vector products. We refer to this
approach as Anandhanarayanan et al.+AD.

  5. (e)

Explicit GPU LSKUM solver.

All the GPU solvers are written in CUDA Fortran. These codes are compiled with
nvcc 23.3 using the flags: -O3 and -mcmodel=large [33]. Table 1 shows the
hardware configuration of the GPU node used for the simulations.

| CPU | GPU
---|---|---
Model | AMD EPYCTM $7532$ | Nvidia Ampere $A100$ SXM4
Cores |  $64$ $\left(2\times 32\right)$ | $5120$
Core Frequency |  $2.40$ GHz |  $1.23$ GHz
Global Memory |  $1$ TiB |  $80$ GiB
$L2$ Cache |  $16$ MiB |  $40$ MiB
Table 1: Hardware configuration used to perform numerical simulations.

###  5.1 Subsonic flow over the NACA 0012 airfoil

The test case under investigation is the subsonic flow over the NACA 0012
airfoil with a freestream Mach number, $M_{\infty}=0.63$ and angle of attack,
$AoA=2^{o}$. Numerical simulations are performed on coarse and very fine point
clouds, with $38,400$ and $1,228,800$ points.

We first investigate the performance of the forward and reverse sweeps
strategies presented in Algorithms 5 and 6. For this purpose, these algorithms
are integrated into the implicit GPU solver based on Anandhanarayanan et al.’s
approach. Note that in the forward and backward sweeps of Algorithm 5, the
wall points with a given color are first updated, followed by interior and
outer boundary points. On the other hand, in Algorithm 6, all points with the
same color are updated together. Figure 1 shows the density residual history
of the implicit GPU solvers on the coarse and fine point distributions. We
observe that the implicit solver with wall points first update in the LU-SGS
sweeps yields a slightly faster rate of convergence compared to the solver
with all points updated together. However, numerical experiments have shown
that the run-time per iteration of the implicit GPU solver with Algorithm 5 is
found to be higher compared to the solver using Algorithm 6. On the coarse
distribution, the run-time is increased by a factor of $1.61$ times, while on
the fine distribution, it is increased by $1.12$ times.

(a) Coarse distribution.  (b) Fine distribution.

Figure 1: Subsonic flow over the NACA 0012 airfoil at $M_{\infty}=0.63$ and
$AoA=2^{o}$. Comparison of residues based on the implicit GPU solvers with the
forward and reverse sweep strategies in Algorithm 5 and 6.

This increase in runtime can be attributed to two main factors, one of which
is the overhead costs. In executing the implicit GPU solver, Algorithm 5 has
six i-loops in forward and backward sweeps together compared to two i-loops in
Algorithm 6. Correspondingly, Algorithm 5 requires six small CUDA kernels
while Algorithm 6 uses only two kernels. Launching a CUDA kernel incurs
overhead, which is notably higher than calling the same function on a CPU. It
is important to highlight that the overhead costs of loop execution and kernel
launches are fixed and do not depend on the point cloud size.

Another source for an increase in the run-time is the number of points with a
given color. Table 2 shows the color groups and the number of points in each
group for both the coarse and fine distributions. Note that the number of
wall, interior, or outer boundary points with a given color is less than the
total number of points with the same color. On the coarse distribution, the
number of points in each color group of Algorithm 5 is small. Due to this, the
implicit solver could not fully utilise the available GPU resources.
Furthermore, the overhead costs are also significant compared to the cost of
floating point operations within a kernel. Both these factors resulted in a
higher run-time factor on the coarse distribution. However, on the fine
distribution, with more points in each color group, the overheads are
relatively smaller than the floating point operations cost. The efficient
utilisation of GPU resources on the fine distribution reduced the run-time
factor significantly. Given these findings, the computational efficiency of
Algorithm 6 makes it the preferred choice for implementation in all versions
of the implicit GPU meshfree LSKUM solver.

| Color group |
---|---|---
| $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$
Level | Number of points in each color group |
Coarse | $9441$ | $9383$ | $9441$ | $9382$ | $433$ | $316$ | $3$ | $1$
Fine | $303491$ | $302393$ | $300101$ | $296968$ | $22526$ | $2912$ | $396$ | $13$

Table 2: Number of points in each color group for the coarse and fine point
distributions.

Figures 2 and 3 compare the residual histories for all the four versions of
the implicit GPU solver and the explicit GPU solver. Figures 2a and 3a show
that the rate of residue fall is higher for the implicit GPU solver based on
Manish et al.+AD approach. From Figures 2b and 3a, we observe that the LU-SGS
algorithm of Manish et al. do not affect the rate of convergence as the
residue plots are almost the same as those from Anandhanarayanan et al.’s
approach. Figure 2c shows that the Anandhanarayanan et al.+AD approach has a
very minimal effect on the convergence rate on the coarse distribution. From
Figure 2a, convergence acceleration on the coarse distribution can be
attributed to the diagonal dominance induced by the exact computation of the
flux Jacobian conserved vector products in the LU-SGS algorithm of eq. (37).
Figure 3a shows the exact computation of Jacobian vector products in eqs. (27)
and (37) has significantly improved the convergence rate on the fine
distribution. Figures 2d and 3b show that a very good speedup in convergence
is achieved using the implicit GPU solver over the explicit solver.

Figure 4 shows the surface pressure distribution $\left(C_{p}\right)$ and the
Mach contours on the finest distribution. The suction peak in the $C_{p}$
distribution is accurately captured, and the Mach contours are smooth.
Furthermore, the $C_{p}$ plots for the implicit and explicit GPU solvers are
identical.

(a)

(b)

(c)

(d)

Figure 2: Subsonic flow over the NACA 0012 airfoil at $M_{\infty}=0.63$ and
$AoA=2^{o}$. Comparison of residues on the coarse distribution with $38,400$
points.

(a)

(b)

Figure 3: Subsonic flow over the NACA 0012 airfoil at $M_{\infty}=0.63$ and
$AoA=2^{o}$. Comparison of residues on the fine distribution with $1,228,800$
points.

(a) Surface pressure distribution

(b) Mach contours

Figure 4: Subsonic flow over the NACA 0012 airfoil at $M_{\infty}=0.63$ and
$AoA=2^{o}$. Surface pressure distribution and Mach contours on the fine
distribution with $1,228,800$ points.

###  5.2 Transonic flow over the NACA 0012 airfoil

This test case represents a transonic flow over the NACA 0012 airfoil with a
freestream Mach number, $M_{\infty}=0.85$, and angle of attack, $AoA=1^{o}$.
The simulations are performed on coarse and fine point distributions with
$38,400$ and $1,228,800$ points.

Figures 5 and 6 show a comparison of the residual histories using the implicit
and explicit GPU solvers on the coarse and fine distributions, respectively.
These plots show that the implicit GPU solver based on the Manish et al.+AD
approach yields a faster convergence rate than the other variants of the
implicit GPU solver. Furthermore, the residual histories based on the
approaches of Manish et al. and Anandhanarayanan et al. are almost identical.
Again, we observe that the exact computation of the flux Jacobian conserved
vector products improved the diagonal dominance of the linear system in eq.
(26), which enhanced the convergence rate. Figures 5d and 6b show a good
speedup is obtained using the implicit solver compared to the explicit solver.

(a)

(b)

(c)

(d)

Figure 5: Transonic flow over the NACA 0012 airfoil at $M_{\infty}=0.85$ and
$AoA=1^{o}$. Comparison of residues on the coarse distribution with $38,400$
points.

(a)

(b)

Figure 6: Transonic flow over the NACA 0012 airfoil at $M_{\infty}=0.85$ and
$AoA=1^{o}$. Comparison of residues on the fine distribution with $1,228,800$
points.

In this test case, the flow results in two shocks of different strengths, one
on the suction side and the other on the pressure side of the airfoil. The
surface pressure distribution in Figure 7a indicates that the discontinuities
are resolved accurately. The pressure contours in Figure 7b show that the
shocks are captured crisply. The computed lift and drag force coefficients are
compared in Table 3. We observe that the force coefficients agree with the
AGARD workshop results [34].

(a) Surface pressure distribution

(b) Pressure contours

Figure 7: Transonic flow over the NACA 0012 airfoil at $M_{\infty}=0.85$ and $AoA=1^{o}$. Surface pressure distribution and pressure contours on the fine distribution with $1,228,800$ points.  Force coefficient | Implicit GPU solver | AGARD
---|---|---
Lift $\left(C_{L}\right)$ | $0.3310$ | $0.3300-0.3890$
Drag $\left(C_{D}\right)$ | $0.0545$ | $0.0464-0.0590$
Table 3: Transonic flow over the NACA 0012 airfoil at $M_{\infty}=0.85$ and
$AoA=1^{o}$. Comparison of the lift and drag coefficients on the fine
distribution with $1,228,800$ points.

###  5.3 Supersonic flow over the NACA 0012 airfoil

The supersonic flow over the NACA 0012 airfoil is computed with a freestream
Mach number, $M_{\infty}=1.2$ and zero angle of attack. Numerical simulations
are performed on coarse and fine point distributions, as used in Section 5.2.

Figures 8 and 9 show the residue plots based on the implicit and explicit
meshfree GPU solvers. Similar to the subsonic and transonic flows, in the
supersonic case, the implicit GPU solver based on the Manish et al.+AD
approach exhibited a superior convergence rate on the coarse and fine
distributions. From Figure 8a, it is evident that this approach resulted in a
better rate of convergence on the coarse distribution compared to subsonic and
transonic flows. Again, the residual histories based on the approaches of
Manish et al. and Anandhanarayanan et al. are almost the same. Figures 8d and
9b show the speedup of the implicit GPU solver over the explicit solver.

(a)

(b)

(c)

(d)

Figure 8: Supersonic flow over the NACA 0012 airfoil at $M_{\infty}=1.2$ and
$AoA=0^{o}$. Comparison of residues on the coarse distribution with $38,400$
points.

(a)

(b)

Figure 9: Supersonic flow over the NACA 0012 airfoil at $M_{\infty}=1.2$ and
$AoA=0^{o}$. Comparison of residues on the fine distribution with $1,228,800$
points.

For this test case, the flow results in a bow shock upstream of the leading
edge and a fishtail shock at the trailing edge of the airfoil. Since the flow
is symmetric, the lift coefficient is zero. Figure 10 shows the surface
pressure distribution $\left(C_{p}\right)$ and the pressure contours on the
finest point distribution. The contour plot shows that the shocks are resolved
very crisply. The $C_{p}$ plot shows the symmetry of the solution. Table 4
compares the lift and drag coefficients on the finest point distribution. We
observe that the computed force coefficients agree with AGARD values [34].

(a) Surface pressure distribution (b) Pressure contours

Figure 10: Supersonic flow over the NACA 0012 airfoil at $M_{\infty}=1.2$ and $AoA=0^{o}$. Surface pressure distribution and pressure contours on the fine distribution with $1,228,800$ points.  Force coefficient | Implicit GPU solver | AGARD
---|---|---
Lift $\left(C_{L}\right)$ | $1.7414\times 10^{-5}$ | $0$
Drag $\left(C_{D}\right)$ | $0.0975$ | $0.0946-0.0960$
Table 4: Supersonic flow over the NACA 0012 airfoil at $M_{\infty}=1.25$ and
$AoA=0^{o}$. Comparison of the lift and drag coefficients on the fine
distribution with $1,228,800$ points.

###  5.4 Subsonic flow over the three-element airfoil

The next test case under investigation is the subsonic flow over the McDonnell
Douglas Aerospace (MDA) $30P$-$30N$ three-element high-lift configuration
[35]. Numerical simulations are performed with a freestream Mach number,
$M_{\infty}=0.2$, and angle of attack, $\alpha=16^{o}$. The computational
domain consists of $52,997$ points. The main airfoil has $696$ points, while
the slat and flap contain $298$ points each. This test case is chosen to see
the effect of AD based modified LU-SGS approaches on the residue fall at low
Mach numbers.

Figure 11a shows the residue plots based on the implicit GPU solvers. We
observe that the implicit solver based on Manish et al.+AD approach exhibited
a slightly superior convergence rate up to $6$ decade residue fall. Figures
11b shows a good speedup in convergence is achieved using the implicit GPU
solver over the explicit solver. Figure 12 shows the surface pressure
distribution $\left(C_{p}\right)$ and the Mach contours. The computed $C_{p}$
is in good agreement with the experimental values [35].

(a)

(b)

Figure 11: Subsonic flow over the MDA three lement airfoil at $M_{\infty}=0.2$
and $AoA=16^{o}$. Comparison of residues based on implicit and explicit GPU
solvers.

(a) Surface pressure distribution (b) Mach contours

Figure 12: Subsonic flow over the MDA three lement airfoil at $M_{\infty}=0.2$
and $AoA=16^{o}$. Surface pressure distribution and pressure contours.

##  6 Performance Analysis of GPU Solvers

In this section, we assess the performance of the implicit and explicit
meshfree GPU solvers. For the benchmarks, the test case of the subsonic flow
over the NACA 0012 airfoil is considered. Numerical simulations are performed
on seven levels of point distributions ranging from $38,400$ to $2,457,600$
points.

In order to measure the performance of the GPU solvers, we define a cost
metric known as the Rate of Data Processing (RDP). The RDP of a GPU solver is
defined as the total wall clock time in seconds per iteration per point in the
computational domain. This research evaluates the RDP values by specifying the
number of fixed point iterations $N$ in Algorithm 3 to $20,000$. Table 5
compares the RDP values on all levels of point distributions.

The table shows that the RDP of the GPU codes decreases with continuous
refinement in the point distribution. We also observe that the implicit
solvers with the exact computation of the products of flux Jacobians and
increments in conserved vectors are computationally more efficient as their
RDP values are smaller than those with approximations to the Jacobian vector
products. Moreover, the efficiency increases with the size of the point
distribution. In particular, the efficiency of the implicit solver based on
the Anandhanarayanan et al.+AD approach increased from $1.05$ times on the
coarse distribution to $1.15$ times on the finest distribution over
Anandhanarayanan et al.’s approach. It is worth noting that the exact
computation of the flux Jacobian conserved vector product requires only one
CUDA kernel. On the other hand, approximating it with an increment in the flux
function requires two kernel calls. From eqs. (28) and (38), we observe that
one of the kernels computes the flux function with the conserved vector
$\boldsymbol{U}$ as the input. The other kernel computes the flux function
with $\boldsymbol{U}+\delta\boldsymbol{U}$ as input. Repeatedly calling these
two kernels over a large number of times increased the overhead costs and thus
resulted in higher RDP values for the implicit GPU solvers based on
Anandhanarayanan et al. and Manish et al. approaches.

Table 5 indicates that the RDP of the implicit GPU solvers is higher than the
RDP of the explicit GPU solver. This can be attributed to the additional CUDA
kernels in implicit solvers that perform the forward and backward sweeps of
the LU-SGS algorithm. From Algorithm 7, it is clear that the more the number
of color groups, the more the calls to the CUDA kernels and thus an increase
in the overheads of the implicit GPU solvers.

| Implicit GPU solver with LU-SGS scheme based on |
---|---|---
Level | Points | Colors | Manish+AD | Manish | Anandhanarayanan+AD | Anandhanarayanan | Explicit |
RDP $\times$ $10^{-8}$ (Lower is better) |  |
$1$ | $38,400$ | $8$ | $20.8373$ | $21.4989$ | $19.5199$ | $20.4602$ | $6.0998$ |
$2$ | $76,800$ | $7$ | $11.0308$ | $11.5563$ | $10.5074$ | $10.9451$ | $3.4449$ |
$3$ | $1,53,600$ | $7$ | $6.5088$ | $6.8088$ | $6.2158$ | $6.4191$ | $2.2051$ |
$4$ | $3,07,200$ | $7$ | $4.0029$ | $4.2636$ | $3.8017$ | $4.0663$ | $1.5822$ |
$5$ | $6,14,400$ | $7$ | $2.7715$ | $3.0830$ | $2.6278$ | $2.8991$ | $1.1963$ |
$6$ | $12,28,800$ | $8$ | $2.2984$ | $2.6486$ | $2.1567$ | $2.4495$ | $1.0273$ |
$7$ | $24,57,600$ | $7$ | $1.8639$ | $2.1985$ | $1.7518$ | $2.0205$ | $0.9266$ |

Table 5: RDP comparison of implicit and explicit meshfree LSKUM GPU solvers.

We also observe that the ratio of the RDP of implicit and explicit solvers
decreases from course to fine point distributions. To investigate this
behaviour in the RDP ratio, we analyse the streaming multiprocessors (SM)
utilisation of the forward and backward sweeps kernels. Tables 6 and 7 compare
the SM utilisation on level $1$ and level $6$ point distributions for the
implcit GPU solver based on Anandhanarayanan et al.+AD approach. We observe
that the SM utilisation depends on the number of points in a color group. The
more the points in a color group the higher the SM utilisation. Although the
point groups $1$ and $3$ on the coarse distribution have the same number of
points, their SM utilisation is different. This is due to the point types
(wall, interior, or outer boundary point) in each group and the number of
neighbours in the connectivity set of a point. Higher values of SM indicate an
efficient usage of GPU resources [26] and thus resulted in smaller values of
RDP ratio.

We want to find an AD based implicit solver that effectively utilises the GPU
resources. Note that the implicit schemes based on Anandhanarayanan et al.+AD
and Manish et al.+AD differ only in their forward sweeps, while having the
same backward sweep algorithms. Therefore, it is enough to compare the SM
utilisation of the CUDA kernels for forward sweep in both the solvers. Table 8
shows the SM utilisation of the forward sweep kernel in the Manish et al.+AD
implicit solver on level $6$ point distribution. From Tables 7 and 8, we
observe that the SM utilisation is more for the implicit solver based on
Anandhanarayanan et al.+AD, indicating better utilisation of GPU resources.

To assess the overall performance of the implicit and explicit GPU solvers, we
define another metric called speedup. The speedup of a GPU code is defined as
the ratio of the RDP of the serial code to the RDP of the GPU code. Note that
the serial version of the implicit LSKUM solver is based on Algorithms 1 and
2. The serial code simulations are performed on a computer node consisting of
$2$ AMD EPYC $7542$ processors with $64$ cores and $256$ GB RAM.

Figure 13 shows the speedup achieved by the GPU solvers. We observe that the
implicit GPU solvers using AD achieved more speedup than the solvers with
approximations to the Jacobian vector products. Furthermore, the rate of
speedup increases as the number of points in the cloud increases. On the
finest distribution, the AD based solvers achieved around $50$ times more
speedup compared to the solvers without AD. Among the AD based GPU solvers,
the implicit approach based on Anandhanarayanan et al.+AD yielded around $25$
times more speedup than the Manish et al.+AD’s approach. We also observe that
the explicit GPU solver shows a better speedup on the coarse and medium point
distributions than the implicit GPU solvers. However, on the fine
distributions, due to the efficient utilisation of GPU resources by the CUDA
kernels for forward and backward sweeps, the AD based implicit solvers
achieved more speedup than the explicit GPU solver.

CUDA | Point color group (Number of points in each color group) |
---|---|---
Kernel | $1(9441)$ | $2(9383)$ | $3(9441)$ | $4(9382)$ | $5(433)$ | $6(316)$ | $7(3)$ | $8(1)$
| SM utilisation (shown in percentange) |
Forward sweep | $11.59$ | $7.02$ | $11.04$ | $7.26$ | $0.50$ | $0.37$ | $0.07$ | $0.16$
Backward sweep | $7.21$ | $9.90$ | $8.07$ | $9.51$ | $0.59$ | $0.31$ | $0.14$ | $0.19$

Table 6: Implicit GPU solver based on Anandhanarayanan et al.+AD: Streaming
multiprocessor (SM) utilisation of the forward and reverse sweep kernels on
the coarse distribution with $38,400$ points.

CUDA | Point color group (Number of points in each color group) |
---|---|---
Kernel | $1(303491)$ | $2(302393)$ | $3(300101)$ | $4(296968)$ | $5(22526)$ | $6(2912)$ | $7(396)$ | $8(13)$
| SM utilisation (shown in percentange) |
Forward sweep | $32.60$ | $31.09$ | $30.69$ | $27.64$ | $15.98$ | $2.76$ | $2.11$ | $2.09$
Backward sweep | $33.60$ | $33.97$ | $35.98$ | $36.54$ | $17.42$ | $2.60$ | $3.35$ | $3.45$

Table 7: Implicit GPU solver based on Anandhanarayanan et al.+AD: Streaming
multiprocessor (SM) utilisation of the forward and reverse sweep kernels on
the fine distribution with $12,28,800$ points.

CUDA | Point color group (Number of points in each color group) |
---|---|---
kernel | $1(303491)$ | $2(302393)$ | $3(300101)$ | $4(296968)$ | $5(22526)$ | $6(2912)$ | $7(396)$ | $8(13)$
| SM utilisation (shown in percentange) |
Forward sweep | $29.48$ | $28.29$ | $28.06$ | $25.19$ | $14.91$ | $2.64$ | $2.03$ | $1.92$

Table 8: Implicit GPU solver based on Manish et al.+AD: Streaming
multiprocessor (SM) utilisation of the forward sweep kernel on the fine
distribution with $12,28,800$ points.

(a) Speedup of the implicit GPU solvers

(b) Speedup of the explicit GPU solver

Figure 13: Speedup achieved by the implicit and explicit LSKUM GPU solvers.

##  7 Conclusions

This report presented the GPU acceleration of implicit kinetic meshfree
methods for the inviscid compressible flows. The meshfree scheme is based on
the Least Squares Kinetic Upwind Method (LSKUM). The implicit schemes were
based on modifications to the meshfree LU-SGS approaches proposed by
Anandhanarayanan et al. and Manish et al. In the modified LU-SGS algorithms,
instead of approximating the product of split flux Jacobians and incremental
conserved vector with the increments in the flux function, they are computed
exactly using algorithmic differentiation (AD). The programming model CUDA was
used to develop the GPU solvers.

The implicit GPU solvers with the exact and approximate computation of
Jacobian vector products in the LU-SGS algorithms for LSKUM were applied to
the standard $2D$ test cases for subsonic, transonic, and supersonic flows.
Numerical results showed that the exact computation of Jacobian vector
products yielded slightly better convergence rate than its approximation. In
particular, the implicit GPU solver based on Manish et al. with the modified
LU-SGS approach exhibited slightly superior convergence rate (in the sense of
slightly better residue fall). In all the test cases, the implicit GPU solvers
achieved a good convergence speed over the explicit GPU solver.

Two metrics, the RDP and speedup, were defined to assess the computational
performance of the implicit GPU solvers. The RDP of a solver was defined as
the total wall clock time in seconds per iteration per point in the cloud. The
speedup was defined as the ratio of the RDP of serial implicit solver to the
RDP of GPU solver. Benchmark simulations on seven levels of point
distributions showed that the implicit GPU solvers with AD based exact
computation of the Jacobian vector products are computationally more efficient
as their RDP values were smaller than the solvers with approximate Jacobian
vector products. On the finest distribution, the AD based implicit GPU solvers
achieved around $50$ times more speedup than those with approximations. Among
the AD based solvers, the implicit approach based on Anandhanarayanan et
al.+AD yielded around $25$ times more speedup than the Manish et al.+AD’s
approach.

In summary, the implicit meshfree scheme based on Anandhanarayanan et al. with
the modified LU-SGS approach is preferred due to its better utilisation of GPU
resources and superior computational efficiency (in terms of much less RDP).
We presume that this approach can yield significant speedups in three-
dimensional flows with more points in the computational domain. Research is in
progress to develop implicit GPU solvers for three-dimensional flows. We also
plan to extend these solvers to laminar and turbulent flows. In the future, we
plan to develop multi-node and multi-GPU versions of these solvers.

##  8 Appendix-A

⬇

1SUBROUTINE flux_Gx(U, Gx)

2!

3 IMPLICIT NONE

4!

5 real*8 :: U(4), Gx(4)

6 real*8 :: rho, u1, u2, pr

7!

8 rho = U(1)

9 u1 = U(2) / rho

10 u2 = U(3) / rho

11 pr = 0.4d0 * (U(4) - 0.5d0 * rho * (u1 * u1 + u2 * u2))

12!

13 Gx(1) = rho*u1

14 Gx(2) = pr+rho*u1*u1

15 Gx(3) = rho*u1*u2

16 Gx(4) = (pr+U(4))*u1

17!

18 END SUBROUTINE

Listing 1: Fortran code to compute the flux vector $\boldsymbol{Gx}$ along the
$x$-coordinate direction.

⬇

1 SUBROUTINE FLUX_GX_D(U, Ud, Gx, Gxd)

2 IMPLICIT NONE

3!

4 real*8 :: U(4), Gx(4)

5 real*8 :: Ud(4), Gxd(4)

6 real*8 :: rho, u1, u2, pr

7 real*8 :: rhod, u1d, u2d, prd

8 real*8 :: temp

9!

10 rhod = Ud(1)

11 rho = U(1)

12 temp = U(2)/rho

13 u1d = (Ud(2)-temp*rhod)/rho

14 u1 = temp

15 temp = U(3)/rho

16 u2d = (Ud(3)-temp*rhod)/rho

17 u2 = temp

18 temp = u1*u1 + u2*u2

19 prd = 0.4d0*(Ud(4)-0.5d0*(temp*rhod+rho*(2*u1*u1d+2*u2*u2d)))

20 pr = 0.4d0*(U(4)-0.5d0*(rho*temp))

21!

22 Gxd = 0.0d0

23 Gxd(1) = u1*rhod + rho*u1d

24 Gx(1) = rho*u1

25 Gxd(2) = prd + u1**2*rhod + rho*2*u1*u1d

26 Gx(2) = pr + rho*u1*u1

27 Gxd(3) = u2*(u1*rhod+rho*u1d) + rho*u1*u2d

28 Gx(3) = rho*u1*u2

29 Gxd(4) = u1*(prd+Ud(4)) + (pr+U(4))*u1d

30 Gx(4) = (pr+u(4))*u1

31!

32 END SUBROUTINE FLUX_GX_D

Listing 2: The black-box tangent linear code to compute the Jacobian vector
product $\boldsymbol{Ax}\delta\boldsymbol{U}$.

⬇

1 SUBROUTINE Jacobian_vector_product_AxdU(U, Ud, AxdU)

2 IMPLICIT NONE

3!

4 real*8 :: U(4), Ud(4), AxdU(4)

5 real*8 :: rho, u1, u2, pr

6 real*8 :: rhod, u1d, u2d, prd

7 real*8 :: temp

8!

9 rhod = Ud(1)

10 rho = U(1)

11 temp = U(2)/rho

12 u1d = (Ud(2)-temp*rhod)/rho

13 u1 = temp

14 temp = U(3)/rho

15 u2d = (Ud(3)-temp*rhod)/rho

16 u2 = temp

17 temp = u1*u1 + u2*u2

18 prd = 0.4d0*(Ud(4)-0.5d0*(temp*rhod+rho*(2*u1*u1d+2*u2*u2d)))

19 pr = 0.4d0*(U(4)-0.5d0*(rho*temp))

20!

21 AxdU(1) = u1*rhod + rho*u1d

22 AxdU(2) = prd + u1**2*rhod + rho*2*u1*u1d

23 AxdU(3) = u2*(u1*rhod+rho*u1d) + rho*u1*u2d

24 AxdU(4) = u1*(prd+Ud(4)) + (pr+U(4))*u1d

25!

26 END SUBROUTINE

Listing 3: The optimised tangent linear code to compute the Jacobian vector
product $\boldsymbol{Ax}\delta\boldsymbol{U}$.

## References

  * [1] H. G. Weller, G. Tabor, H. Jasak, and C. Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  Computers in Physics, 12(6):620–631, November 1998.
  * [2] Thomas D. Economon, Francisco Palacios, Sean R. Copeland, Trent W. Lukaczyk, and Juan J. Alonso.  Su2: An open-source suite for multiphysics simulation and design.  AIAA Journal, 54(3):828–846, March 2016.
  * [3] Mark J. Harris, William V. Baxter, Thorsten Scheuermann, and Anselmo Lastra.  Simulation of cloud dynamics on graphics hardware.  In HWWS 03: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware, pages 92–101, 2003.
  * [4] F.D. Witherden, A.M. Farrington, and P.E. Vincent.  PyFR: An open source framework for solving advection–diffusion type problems on streaming architectures using the flux reconstruction approach.  Computer Physics Communications, 185(11):3028–3040, nov 2014.
  * [5] Matteo Bernardini, Davide Modesti, Francesco Salvadore, Srikanth Sathyanarayana, Giacomo Della Posta, and Sergio Pirozzoli.  Streams-2.0: Supersonic turbulent accelerated navier-stokes solver version 2.0.  Computer Physics Communications, 285:108644, April 2023.
  * [6] S.M. Deshpande, J. C. Mandal, and A.K. Ghosh.  Least squares weak upwind method for Euler equations.  FM Report No. 1989-FM-4, Dept. of Aerospace Engg., Indian Institute of Science, Bangalore, 1989.
  * [7] J.T. Batina.  A gridless Euler/Navier-Stokes solution algorithm for complex-aircraft applications.  AIAA Paper 1993-0333, 1993.
  * [8] Rainald Löhner, Carlos Sacco, Eugenio Oñate, and Sergio Idelsohn.  A finite point method for compressible flow.  International Journal for Numerical Methods in Engineering, 53(8):1765–1779, December 2001.
  * [9] K. Morinishi.  An implicit gridless type solver for the navier-stokes equations.  Computational Fluid Dynamics Journal, pages 551–560, 2001.
  * [10] D. Sridar and N. Balakrishnan.  An upwind finite difference scheme for meshless solvers.  Journal of Computational Physics, 189(1):1–29, July 2003.
  * [11] Aaron Katz and Antony Jameson.  Multicloud: Multigrid convergence with a meshless operator.  Journal of Computational Physics, 228(14):5237–5250, August 2009.
  * [12] Edmond Kwan-yu Chiu, Qiqi Wang, Rui Hu, and Antony Jameson.  A conservative mesh-free scheme and generalized framework for conservation laws.  SIAM Journal on Scientific Computing, 34(6):A2896–A2916, January 2012.
  * [13] A. K. Ghosh and S. M. Deshpande.  Least squares kinetic upwind method for inviscid compressible flows.  AIAA paper 1995-1735, 1995.
  * [14] S. M. Deshpande, P. S. Kulkarni, and A. K. Ghosh.  New developments in kinetic schemes.  Computers Math. Applic., 35(1):75–93, 1998.
  * [15] S M. Deshpande, K. Anandhanarayanan, C. Praveen, and V. Ramesh.  Theory and application of 3-D LSKUM based on entropy variables.  Int. J. Numer. Meth. Fluids, 40:47–62, 2002.
  * [16] V. Ramesh and S.M. Deshpande.  Unsteady flow computations for flow past multiple moving boundaries using lskum.  Computers & Fluids, 36(10):1592–1608, December 2007.
  * [17] A.K. Mahendra, R.K. Singh, and G. Gouthaman.  Meshless kinetic upwind method for compressible, viscous rotating flows.  Computers & Fluids, 46(1):325–332, July 2011.
  * [18] K. Anandhanarayanan, Konark Arora, Vaibhav Shah, R. Krishnamurthy, and Debasis Chakraborty.  Separation dynamics of air-to-air missile using a grid-free euler solver.  Journal of Aircraft, 50(3):725–731, 2013.
  * [19] K. Anandhanarayanan, R. Krishnamurthy, and Debasis Chakraborty.  Development and validation of a grid-free viscous solver.  AIAA Journal, 54(10):3312–3315, 2016.
  * [20] V. Ramesh and S. M. Deshpande.  Kinetic mesh-free method for flutter prediction in turbomachines.  Sadhana, 39(1):149–164, 2014.
  * [21] Seokkwan Yoon and Antony Jameson.  Lower-upper symmetric-gauss-seidel method for the euler and navier-stokes equations.  AIAA Journal, 26(9):1025–1026, September 1988.
  * [22] Dmitri Sharov and Kazuhiro Nakahashi.  Low speed preconditioning and lu-sgs scheme for 3-d viscous flow computations on unstructured grids.  In 36th AIAA Aerospace Sciences Meeting and Exhibit. American Institute of Aeronautics and Astronautics, January 1998.
  * [23] K. Anandhanarayanan, M. Nagarathinam, and S.M. Deshpande.  An entropy variable based gridfree solver with LU-SGS accelerator.  ICAS paper 2004-2.4.3, 2004.
  * [24] Manish K. Singh, V. Ramesh, and N. Balakrishnan.  Implicit scheme for meshless compressible euler solver.  Engineering Applications of Computational Fluid Mechanics, 9(1):382–398, January 2015.
  * [25] A. Griewank and A. Walther.  Evaluating derivatives: Principles and techniques of Algorithmic Differentiation.  SIAM, 2008.
  * [26] Nischay Ram Mamidi, Dhruv Saxena, Kumar Prasun, Anil Nemili, Bharatkumar Sharma, and S. M. Deshpande.  Performance analysis of gpu accelerated meshfree q-lskum solvers in fortran, c, python, and julia.  In 2022 IEEE 29th International Conference on High Performance Computing, Data, and Analytics (HiPC), pages 156–165, 2022.
  * [27] R. Courant, E. Issacson, and M. Rees.  On the solution of nonlinear hyperbolic differential equations by finite differencesew developments in kinetic schemes.  Comm. Pure Appl. Math., 5:243–255, 1952.
  * [28] J. C. Mandal and S. M. Deshpande.  Kinetic flux vector splitting for Euler equations.  Comp. & Fluids, 23(2):447–478, 1994.
  * [29] S. M. Deshpande, K. Anandhanarayanan, C. Praveen, and V. Ramesh.  Theory and application of 3-D LSKUM based on entropy variables.  Int. J. Numer. Meth. Fluids, 40:47–62, 2002.
  * [30] S. M. Deshpande.  On the maxwellian distribution, symmetric form, and entropy conservation for the euler equations.  NASA-TP-2583, 1986.
  * [31] L. Hascoët and V. Pascual.  The Tapenade Automatic Differentiation tool: Principles, Model, and Specification.  ACM Transactions On Mathematical Software, 39(3), 2013.
  * [32] NVIDIA Corporation.  Developer Tools Documentation.  2021\.
  * [33] David Kirk.  NVIDIA CUDA software and GPU parallel computing architecture.  volume 7, pages 103–104, 01 2007.
  * [34] H Viviand.  Test cases for inviscid flow field methods.  AGARD AR-211, 1985.
  * [35] Christopher L Rumsey, Elizabeth M Lee-Rausch, and Ralph D Watson.  Three-dimensional effects in multi-element high lift computations.  Computers & Fluids, 32(5):631–657, June 2003.