In strongly coupled uid-structure interaction simulations, the uid dynamics and solid dynamics problems are solved independently on their own meshes. Therefore, it becomes necessary to interpolate physical properties (pressure, displacement) across two meshes. In this work, we propose to accelerate the interpolation process by the method of radial basis functions using the matrix-free solution of the equation system on a GPU.
Main interpolation methods on non-matching meshes for uid-structure
interaction (FSI) simulations are overviewed in [
The choice of the function determines the conditioning and density of the matrix, and, as a result, the computational complexity of solving the system of equations. The advantages of the RBF interpolation are the following: { it does not require mesh-connectivity information; { it requires solving a small system of equations, especially with the compact basis functions; { it can be e ciently parallelized.
This paper is structured as follows. Section 2 brie y describes the RBF interpolation scheme for the FSI problem. The next section presents a new approach ? This work is supported by the Russian Foundation for Basic Research (projects: 16-37-00060-mol a, 16-01-00129-a, 17-01-00402-a). based on layer-by-layer mesh partitioning for reducing the problem size. The fourth section describes a matrix-free solution of the interpolation problem on a GPU. 2
Consider the problem of interpolation by the method of radial basis functions, in the context of pressure interpolation in the FSI problem. Let be the mesh approximation of the domain with the boundary = @ and the given pressure p . The mesh approximation of the domain with the interpolated pressure is denoted by with the boundary = @ . The domains and have a common part of boundary (interface part), i.e. \ 6= f;g. The pressure interpolation between these meshes can be expressed in the matrix form as follows: W
P T
P 0 = p 0 or
A = b ; and the target pressure vector p is obtained by the matrix-vector product (1) (2) p = [W where W , W are the n n and n n matrices, consisting of the elements with values equal ( kxi xj k2 ) and ( kxi xj k2 ) respectively; n and n are the numbers of interpolation points on the interface boundary of the domains; ; are the coe cients of the interpolant; xi | the vector of coordinates of the interpolation points.
The function (kxk2) reduce to a scalar function of the Euclidean norm kxk2 of their vector argument x, i.e.: they are radial in the sense (kxk2) = (r); x 2 IR3 for the \radius" r = kxk2 with a scalar function : IR ! IR. This makes their use for highdimensional reconstruction problems very e cient, and it induces invariance under orthogonal transformations.
There are two types of radial basis functions: basis functions with global and compact support. The basis functions with global support are Gaussian ( (r) = e c r2 ; r 0; c > 0), inverse multiquadric ( (r) = (r2 + c2) 1=2; r 0; c > 0), thin plate spline ( (r) = r2 log r; r 0), cubic ( (r) = r3; r 0). The basis functions with compact support have the form: (r) = (c r)2; r 0; c > 0, (r) = (c r)4 (4r + 1); r 0; c > 0, etc. Further, the constant c is considered from the point of view of the in uence on the error and the rate of convergence of the solution of the system (1) on the example of the global basis function (r) = e c r2 .
Solving the system of equations (1) is the most computationally expensive
part of the interpolation. In [
In this section, we demonstrate reducing the size of the system of equations for
the uid-structure interaction of a supersonic ow with a nozzle wall that has a
high geometric expansion ratio [
Let us divide the interface part of the mesh into layers. To do this, we use the neighborhood criterion where any two mesh cells are considered adjacent if they have at least one common node.
The mesh is the discrete approximation of the rotation surface with the
closed directrix. To form layers in parallel to the directrix (see Fig. 1 (b)) or
along the surface generatrix (see Fig. 1 (c)), we use the algorithm proposed in
[
The quality of interpolation is compared for the global (r) = e c r2 , using di erent partitions, numbers of layers and constant c. The quality of the pressure interpolation can be estimated as the relative error computed by the ratio of the norms of the resultant forces of the pressure on the interface boundary.
Table 1 shows the results for the pressure interpolation in parallel to the directrix (Radial partitioning) and along the surface generatrix (Longitudinal partitioning). In the last column, the evaluation of the interpolation quality is given for all possible interpolation points (28800) of . Therefore, the results shown in this column do not depend on the partitioning.
The quality of the interpolation with the data reduction depends not only on the number of interpolation points but also on the choice of the points (Fig. 1). In addition, the table 1 the convergence of the solution of the system (1) depends on the choice of interpolation points. So, when choosing points based on longitudinal partitioning, the number of iterations for solving the resulting system of equations is twice as large as for a radial partition, regardless of the number of equations. With an increase in the constant c, the number of iterations decreases, but the interpolation error increases. A similar situation is typical for local basis functions. The best interpolation is achieved for the radial partitioning of the domain with c = 0:1. It allows to reduce the number of equations in system (1) by a factor of 15 with the acceptable quality of the interpolation. 4
One of the speci c features of the system (1) is a dense matrix, which imposes
some restrictions on the GPU use due to the small capacity of the available
GPU memory. The problem can be resolved by (i) using several GPUs, thereby
increasing the total memory available for the system solution; (ii) solving the
system of equations without the formation of a matrix (Matrix-Free Algorithm).
In this case, the matrix elements are computed as they are required in the
algorithm of the system solution. The solution of the system by the RBF method
is possible without the formation of a matrix, since the matrix elements are
computed by the chosen basis function. This improves the data locality and
arithmetic intensity for matrices and vectors. The memory requirements and
CPU-GPU communications are reduced. The e ciency of the algorithm can be
improved if multi-GPUs are used in the similar way to that in [
Let us consider in more detail the MFA computing expenses. Table 2 shows the time of the sequential and parallel formation of the matrix A of the system (1). In the MFA, the formation time is excluded. For comparison, the time of the solution of the system with an assembled matrix is given. The time of copying the matrix A of the system (1) to the GPU memory is also presented. In addition, the time is given for solving the system with the use of both the algorithm with an assembled matrix and the matrix-free solution algorithm.
The CPU parallelization is carried out with OpenMP. The solution of the
system of equations on several GPUs is carried out by CUDA in conjunction with
OpenMP. The system of equations is solved by the conjugate gradient method
with the diagonal preconditioner [
When the system of equations is solved using the assembled matrix on the CPU, the step of the matrix formation is added. The use of GPU increases the cost due to the necessity of copying the data to the GPU. When the system is solved using the MFA the cost is not increased because there in no need to copy the data. In the last line of Table 2, the total time is given for each of the above approaches.
The numerical computations show that the use of eight CPU threads within one computing node reduces the solution time almost by a factor of seven. One GPU allows to speed up solving the system by a factor of 250 compared with one CPU thead and by a factor of 50 compared with 8 CPU. The GPU e ciency increases with the increase of the system size. Using two GPUs reduces the time by a factor of 1.5 compared with one GPU and by a factor of 350 compared with the CPU. With an increase in the number of GPUs, the strong scalability can be provided only when the sizes of the submatrices on each GPU are preserved.
The matrix-free solution of the system using 8 CPU reduces the solution time by a factor of 2.5. However, the solution with the assembled matrix is twice as fast as the matrix-free solution. When one GPU is used, the time for the matrix-free solution of the system of equations is 5 times larger than that for the solution with the assembled matrix, and in the case of using two GPUs, the matrix-free solution is 3 times longer. The speedup obtained at the use of one CPU thread is 55 times smaller than that when using one GPU and 110 times smaller than that when using two GPUs. It should be noted that the use of local basis functions with an introduced radius of in uence increases the MFA e ciency.
Let us estimate the maximum size of the system, which can be solved using the MFA on a one GPU. For the matrix A formation, the coordinates of the interpolation points are used. Then for interpolation in a three-dimensional space, it is necessary to allocate memory for the vector of coordinates of length equal to n 3. The required memory size for solving the system with the assembled matrix is n n . The remaining vectors participating in the conjugate gradient method coincide for both algorithms. Thus, the memory size for interpolating the mesh data in the three-dimensional space is decreased by a factor of n =3. The maximum system size solved by the MFA increases by the same factor. The algorithm of the conjugate gradient method with a diagonal preconditioner involves the use of memory to store a matrix of size n 3 (the MFA) and six vectors n 1. Thus, for solving the system using the MFA and double precision arithmetic, n (3 + 6) 8 bytes are required. Consequently, the maximum size of a system for the GPU with a 4Gb GDDR is about 6 108 equations. Using two graphics cards, the possible size of the system is increased to 1:2 109 equations. Thus, for the dense matrices obtained on the basis of global basis functions, a parallel method of conjugate gradients is constructed. The computations are distributed among several GPUs. The use of the matrixfree approach makes it possible to remove any limitations on the amount of memory. 5
For solving the problems on unstructured meshes, choosing the data points affects the quality of interpolation. The obtained results show that the interpolation on the irregularly distributed data is of high quality and applicable for meshes of super-large dimensions. The variation of the constant c allows us to choosing the optimal ratio of the interpolation error and the time of its construction.
Using a matrix-free algorithm on large meshes signi cantly reduces the memory costs associated with the formation of the interpolation matrix. At the same time, the computation locality of the matrix-vector product computations increases when solving the system of equations by iterative methods. The solution of systems with dense matrices by the MFA on the CPU does not lead to any signi cant time reductions. Since the time of the matrix formation is less than 1% of the solution time, the use of the MFA in conjunction with the CPU is ine cient. The matrix-free approach is most e ective when using a GPU, especially when it is not possible to achieve a large reduction of points without the interpolation quality loss. Using a GPU for solving larger systems of equations allows minimizing the cost of additional computations associated with the formation of the matrix elements.