Parallel Optimization with Boundary Elements and Kernel Independent Fast Multipole Method

Topological-shape optimization (layout optimization) problems are ubiquitous in modern engineering. The problem of topological-shape optimization can be stated as follows - for a given initial domain and boundary value problem, find a sub-domain of a given volume that extremizes the cost functional depending on the solution of a boundary value problem for this subdomain. The most common formulation studied in structural engineering (e.g. [1]) seeks for an optimal shape and topology of an elastic body that minimize the strain energy (compliance) while satisfying the weight constraint and the additional constraints imposed by the boundary value problem. Other examples of such optimization include the number of problems involving Navier, Navier-Stokes, Laplace, Maxwell and Helmholtz equations [2-5]. Non-convex nature of such optimization problems often makes finding globally optimal layouts difficult or impossible.

It is worth noting that recent emergence and development of additive technology [6, 7] wiped out previous technological limitations on admissible design shapes, and led to the situation when (at least on a prototype stage) any shape can be immediately produced at a cost that does not depend on the complexity of the shape and only on the amount of the material spent. This has revived the interest in large-scale three-dimensional techniques of shape and topology optimization.

A number of efficient shape optimization techniques have been developed in the last few decades, including boundary variation algorithms [8], level set methods [9], homogenization [10]. These techniques are usually implemented as extensions to the existing finite element (FEM) modeling packages. Modern FEM parallel optimization frameworks display outstanding performance, enabling simulations with tens and hundreds of millions of degrees of freedom (DOF) [11].

If the problem statement excludes composite layouts, non-homogeneous material properties and body forces, the shape and topology of the domain, as well as the boundary value problem, are fully defined by the domain boundary. Therefore, boundary element method (BEM) and boundary integral approaches appear to be the most appropriate tools for such an analysis. Few recent works [12, 13] have demonstrated that BEM can be a convenient tool for topological-shape optimization. Significant work has been done in the adjacent field of shape sensitivity analysis in elastodynamics using fast multipole BEM [14].

In our recent works [15, 16] we have demonstrated the application of 2D and 3D boundary element codes equipped with fast algebraic solvers to the problem of topological-shape optimization of elastic structures. However, the tools we used did not allow us to reach true scalability, and evaluate the capabilities of our approach on large problems. In our current work, we present topological-shape optimization framework based on extremely scalable implementation of a kernel-independent fast multipole method (KIFMM) [17]. The sin-gle-node performance of our code is comparable to the performance of state-of-the-art FEM solvers, and can be further improved. Fully parallel implementation of KIFMM used in our work [18] guarantees efficient parallelization up to 20k cores. Our present paper discusses only the problems of compliance optimization in elasticity. However, due to kernel-independence, our technique can be straighforwardly employed in many other fields of engineering optimization dealing with corresponding partial differential equations and cost functionals.

The remainder of the paper has the following organization. Section 2 gives the description of our approach. Section 3 presents and briefly discusses the results of benchmark numerical simulations. Section 4 outlines the directions of further work.

One of the key points of our approach is the use of surface integral equations and boundary discretization for the solution of the elasticity boundary value problem at every iteration. We utilize direct boundary integral equation formulation for elasticity [19]:

where $u_i(x)\left(p_i(x)\right)$ is the displacement (traction) at the point on a boundary of the domain, and $U_{i j}(x, \xi)\left(P_{i j}(x, \xi)\right)$ are corresponding fundamental solutions. For the case of a linear isotropic elastic material with the shear modulus $G$ and Poisson's ratio $v$ they are given by:

where $r=|\xi-x|$. Once the solution on the boundary is found, the stress at the point inside the domain can be calculated according to:

where $D_{k i j}(p, \mathrm{x})$ and $S_{k i j}(p, \mathrm{x})$ are given by:

The numerical treatment of this formulation is presented in subsection 2.2. The optimization procedure utilizes the concept of topological derivative (TD) - a cost of making an infinitesimal spherical cavity with a center in a given point of the domain. For the case of strain energy (compliance) cost functional $\psi=\frac{1}{2} \int_{\Omega} \sigma_{i j} \epsilon_{i j} d \Omega$ and 3D linear isotropic elasticity the analytical expression of TD is available [1]:

where $E$ is the material Young's modulus. Optimization strategies employing TDs are usually based on the progressive elimination of the material. In this paper, we present a relatively simple approach based on voxel representation of materials and discrete elimination of voxels. The parallel framework of such optimization is presented in subsection 2.3.

Using the method of collocations and piecewise-constant approximation [19], one can discretize Eq. (1) into the following system of linear equations:

where $\Delta P_{i j}=\int_{S_k} P_{m n}\left(\xi_1, x\right) d S_k, \Delta U_{i j}=\int_{S_k} U_{m n}\left(\xi_1, x\right) d S_k, i=3 k+m, j=m l+n, \xi l$ is $l$-th collocation point and $S_k$ is the area of $k$-th flat triangular boundary element.

After a rearrangement of columns in 8 , we obtain the following system of linear equations, where all unknowns appear in vector $x$, while the tractions or displacements known from boundary conditions appear in vector $y$.

We utilize parallel implementation of the GMRES algorithm [20] available in PETSC library [21, 22] to solve this system of linear equations. The performance of GMRES is conditioned by the performance of generic matrix-vector multiplication procedure, that, for general matrix shape takes $O\left(N^2\right)$ operations. In order to reduce the computational cost of matrix-vector products we employ the recent implementation of the KIFMM - PVFMM [17, 18]. PVFMM provides a "black-box" tool to perform the rapid summation of the kind:

Here $t_{\mathrm{i}}$ is the vector of "target" values being computed, that are positioned at spatial points $x_i ; s_j$ is the vector of known "source" values positioned at points $y_j$. Kernel function $K_{i j}\left(X_i\right.$, $Y_j \mathrm{n}_{\mathrm{j}})$ depends on both source and target coordinates, and, optionally, on the normal $\mathrm{n}_{\mathrm{j}}$ that is specified at a source point. Valid kernel functions belong to a wide class of fundamental solutions of second order partial differential equations. In particular, KIFMM is suitable for fast summation of the matrix components of kernels $2,3,5,6$. PVFMM is highly optimized and extremely scalable implementation of KIFMM, it supports both intranode OMP parallelization and internode MPI parallelization, demonstrating excellent scalability for up to tens of thousands of cores. An in-depth description of PVFMM is beyond the scope of this paper, the reader is referred to [17] for the more detailed description of the method and to [18] for peculiarities of the parallel implementation used.

We utilize a KIFMM summation algorithm for fast calculation of the matrix-vector products in 9 . The summation 10 is performed with "source" coordinates $\mathrm{y}_j$ being the coordinates of quadrature points on triangular boundary elements, "source" values $s_j$ being the constant approximation of the solution on the triangular element, normalized by triangle's area and quadrature point's weight, and "target" values $t_i$ being the values at collocation points, positioned at the centers of triangles. Each triangular element contains 16 quadrature points with corresponding weights. Generic matrix-vector product operation utilized by GMRES is performed in the following way. First, we perform rapid summation $A \cdot x$ using PVFMM. Then follows so-called "local pass" - we subtract inaccurate contributions from the sources corresponding to quadrature points on the same triangle that contains a collocation point, and replace those with the analytical expressions for singular integrals obtained in [19]. This relatively simple scheme is easily parallelizable, since it does not require MPI communications during the "local pass". However, it imposes limitations in terms of achievable model sizes - for large enough models the numerical summation of near-singular integrals over triangles neighboring to a triangle containing collocation point becomes inaccurate and leads to significant artifacts in the surface solution. Since the surface solution is found, the solution in stresses at internal points is found via PVFMM summation of the kernels 5,6 . It is worth noting here, that the sampling of the points inside the domain should not necessarily be uniform - one could successfully use adaptive strategies of sampling points inside the domain, which reduces the computational complexity of the domain computation to $O\left(N_s\right)$, where $N_s$ is the number of surface elements [15].

The optimization procedure performs the following cycle. We start with three-dimensional array of voxels, at the first iteration all the voxels are material. We contour all material voxels with boundary elements, and set up the initial boundary value problem. Since this part is not computationally intensive, it can be done serially. Then the surface and volume meshes that were generated on a master process are scattered over all MPI processes, and parallel solution of boundary integral equation, as well as computation of topological derivatives inside the domain is performed. After that the field of topological derivatives is gathered to a master process, and isosurface of a topological derivative is extracted by thresholding the field with a reasonable threshold value. All the voxels that are below the threshold are assigned to be void. After post-processing procedure excluding isolated voxels and surface irregularities we compute the value of the cost functional $\Psi$ and the ratio $\alpha$ between the current and the initial number of material voxels. If the critical ratio $\alpha_c$ was not reached - proceed to the next iteration.

We test our method on two classical examples of compliance optimization - cantilever supports. We use these benchmarks to measure the performance of our optimization, as well as parallel scalability. We start with a unit cube with one of its sides fixed and two different load configurations applied to the side opposite to a fixed one (Fig. 1a and d). The material properties are $E=1, v=0.3$. The initial model contains $64^3$ voxels, $0.79 \cdot 10^6$ volume DOFs and $4.72 \cdot 10^6$ surface DOFs. Computations were carried out on a cluster machine with up to four nodes, 16 cores each. In both cases optimization iteration took about 15 minutes on 64 cores. This time was mostly conditioned by relatively slow GMRES iterative process, that took about 100-200 iterations per each boundary value problem. However, the volume computation took less than $20 \%$ of that time, the performance of volume computation was about 2,000 points per second at each core, which is close to the performance of state of the art FEM solvers [11]. Although the computationally cheap parts of the code (meshing and post-processing of topological derivatives) were not parallelized, our code displayed nearly linear scalability for up to 64 cores. Figure 1b and e demonstrate the solutions obtained after three iterations. The level of details in the final solution depends on the threshold of the topological derivative and the number of iterations. However, both obtained solutions are in qualitative agreement with 2D and 3D solutions of similar problems obtained earlier [15, 16, 23]. The quality of the surface of the optimal configuration is improved with Laplassian smoothing [24] post-processing step. Figure 1c and f gives the evolution of the cost functional as a function of the current volume ratio for both examples. Figure 2 a and b demonstrate the evolution of the field of topological derivatives in a symmetry plane cross section during the optimization process.

As has been discussed above, our fast BEM technique provides the efficient and universal framework for topology optimization that is capable to solve problems with millions of design DOF. Our technique does not suffer from the problems inherent to FEM homogenization approaches, such as "checkerboard instability" [23], and provides better asymptotic performance, since volume computation can be done adaptively. Also our method inherits such a strong feature of boundary integrals techniques as efficient analysis of infinite (external) problems. Availability of built-in periodic boundary conditions in PVFMM enables the analysis of periodic elastic structures.

One can point out a few directions for further improvement of our method. The highest priority is the improvement of the performance and precision of the boundary integral equation solver. This can be achieved by using higher order approximation, appropriate preconditioners, improvement of the mesh generator and higher quality of the surface integration.

Another important direction is the improvement of the quality of the domain solution. The evaluation of topological derivatives near the boundary suffers from the lack of precision. One approach to the problem is the surface mesh refinement for the points inside the domain that approach the domain boundary. However, this way appears to be computation-ally expensive. The alternative approach is to construct local extrapolation of the field of the topological derivative based on the valid points sufficiently far from the boundary and use it to determine the values of the topological derivative in the points near the domain boundary.

So far we have demonstrated the efficiency of our technique for the case of compliance minimization in linear elasticity. However, the proposed approach has significant flexibility, and due to its kernel independence it is able to solve similar optimization problems in a number of engineering fields. Our further work is aimed at taking advantage of the strong features of our method (ability to solve external infinite problems, ease of incorporation of periodic boundary conditions) in the problems of optimization of periodic cells of metamaterials.