Dielectric breakdown prediction with GPU-accelerated BEM
Cedric Münger Steffen Börm Jörg Ostrowski
math.NA cs.NA 65N38, 35J05, 65D30


Dielectric Breakdown Prediction with GPU-Accelerated BEM





Cedric Münger Seminar for Applied Mathematics, ETH Zürich, Switzerland,
cmuenger@student.ethz.ch
  Steffen Börm Math. Seminar, Christian-Albrechts-Univ. Kiel, Germany,
  boerm@math.uni-kiel.de
  Jörg Ostrowski ABB Corporate Research, Baden, Switzerland,joerg.ostrowski@ch.abb.com






*
    Cedric Münger 
Steffen Börm 
Jörg Ostrowski 

    December 30, 2023
=================================================



The prediction of a dielectric breakdown in a high-voltage device is based on criteria that evaluate the electric field along possible breakdown paths. For this purpose it is necessary to efficiently compute the electric field at arbitrary points in space. A boundary element method (BEM) based on an indirect formulation, realized with MPI-parallel collocation, has proven to cope very well with this requirement. It deploys surface meshes only, which are easy to generate even for complex industrial geometries. The assembly of the large dense BEM-matrix, as well as the iterative solution of the resulting system, and the evaluation along the field-lines all require to carry out the same type of calculation many times. Graphical Processing Units (GPUs) promise to be more efficient than Central Processing Units (CPUs) when it is possible to do the same type of calculations for large blocks of data in parallel. In this paper we therefore investigate if GPU acceleration is a measure to further speed up the established CPU-parallel BEM solver.





§ INTRODUCTION







Every high-voltage device has to pass dielectric type tests, in which a
large voltage is applied to the device.
The test is passed if no dielectric breakdown occurs.
A breakdown usually starts from an electrode-surface with high dielectric
stress, and then propagates through the volume along a field-line of the electric field  towards the opposite
electrode, see Fig. <ref>. The propagation stops if the electric field along this breakdown
path γ is not strong enough. For details see <cit.>.
An inception of a streamer, i.e. the initial state of a breakdown only occurs if the criterion

    ∫_γα_eff(||) ds > K_str

is fulfilled. Here α_eff is the effective ionization function that depends on the strength of the electric field ||, and K_str is the (gas-specific) streamer constant. The prediction of a dielectric breakdown during a type test relies on the evaluation of this criterion along the most probable breakdown paths. It requires the computation of the electric field at all surface points and along field lines in the volume.


Simulation-based dielectric design became a standard procedure because a user-friendly, i.e., fast, robust, reliable, and easy-to-use computational method was developed, see <cit.>. In the following we will first describe this boundary-element-based method and then introduce how general-purpose graphics processing units (GPGPUs) can be used to massively reduce computing times.

0.49


    < g r a p h i c s >

figureThe electric field strength on the surface of a disconnector and possible breakdown paths along the field lines.  


0.49

    < g r a p h i c s >

figureExample of a device that consists of a conductor Ω_0 that is connected to a prescribed electric potential V_0, a floating conductor Ω_1, an insulator Ω_2 , and the unbounded exterior domain Ω_3. 






§ BEM FORMULATION



In this section we derive the BEM-formulation as it is in use since many years at
ABB, see <cit.>.
The device consists of the subdomains
Ω_0,…,Ω_m-1, and the unbounded exterior subdomain is  Ω_m=ℝ^3∖⋃_k={0..m-1}Ω_k, see the example in Fig. <ref>. The electric field =- φ is
calculated by solving the Laplace equation

    ϵφ = 0

for the electric scalar potential φ in each of these subdomains. The permittivity is denoted by ϵ.
We use an indirect formulation with a single-layer potential

    φ() = Ψ_SL[σ]()
      = ∫_∂Ωσ()/4π |-|  dS_y ,

and search for the unknown scalar virtual surface charge density σ that is related to the physical surface charge density σ_s. Each conductor, i.e. each separated conducting part with electrical conductivity σ_el > 0, is on a constant electrical potential. If a conductor is connected to an electric potential V_0, like Ω_0 in Fig. <ref>, then it holds

    φ()=V_0   ∀∈Ω_0.

The electric potential V of floating conductors like Ω_1 in Fig. <ref> is unknown

    φ()=V   ∀∈Ω_1,

and is to be determined by a charge neutrality condition, see <cit.>. The total charge Q of the floating conductor can be derived from the Gauss law as

    Q=∫_∂Ω_1σ_s  dS=∫_∂Ω_1·  dS.

The normal component of the displacement field  can be expressed as

    · = ε^+ · = -ε^+ φ· = -ε^+ Ψ_SL [σ] ·.

Here ε^+ denotes the permittivity of the exterior domain.
The Neumann trace of the single layer potential and can be expressed with help of the adjoint double layer K'

    Ψ_SL [σ] · =  1/2σ + K' σ
          K'(σ)() = ∫_∂Ω -/4 π  |  - |^3·() σ() dS_y .

Combining the equations (<ref>)-(<ref>)  with Q=0 due to charge neutrality yields

    ∫_∂Ω_11/2ε^+ σ () + ε^+ (K'σ)() dS_y  = 0.

We model thin floating conductive sheets only by a single surface. Then the electric fields from both sides (± ) need to be considered for charge neutrality, since

    σ_s=·( ^+ -^-) ⟹
          ∫_∂Ω_11/2 (ε^++ε^-) σ () + (ε^+-ε^-) (K'σ)() dS_y  = 0.

There is no surface charge on non-conductors: σ_s=0 on ∂Ω_2

    1/2(ε^+ + ε^-) σ() + (ε^+ - ε^-)(K'σ)() = 0   ∀∈∂Ω_2.

So for our simple but quite general example of Fig. <ref> we have to solve the following set of equations:

    ∫_∂Ωσ()/4π |-|  dS_y    =    V_0   ∀∈∂Ω_0
    ∫_∂Ωσ()/4π |-|  dS_y-V    =    0   ∀∈∂Ω_1
    ∫_∂Ω_11/2ε^+ σ() + ε^+  (K'σ)() dS_y     =    0
    1/2(ε^+ + ε^-) σ() +  (ε^+ - ε^-)(K'σ)()    =    0   ∀∈∂Ω_2

The solution of the system of equations (<ref>) - (<ref>) yields  the virtual surface charge distribution from which the electric field can be compute at any point in space as

    () = ∫_∂Ω -/4 π |  - |^3σ() dS_y   ∀∈ℝ^3.






§ DISCRETIZATION




We use a collocation boundary element approach: the surface ∂Ω is
represented by a collection of triangles τ_1,…,τ_N
with vertices _1,…,_n.
The unknown function σ

    σ_h() = ∑_j=1^n u_j ψ_j(),

is approximated by suitable basis functions ψ_1,…,ψ_n. This approach is to be inserted into (<ref>)-(<ref>) and this set of equations  is only required to hold in the
collocation points _1,…,_n.
This is an n+N_fl-dimensional system of linear equations, with N_fl being the number of floating conductors. The implementation of this approach poses a number of challenges:


  * High-voltage devices have smooth curved surfaces in order to avoid field enhancements. We use piecewise quadratic parametrizations with (curved)
  triangles τ_1,…,τ_N to minimize the geometrical discretization error.

  * The entries of the matrix corresponding to the linear
  system have to be computed.
  In the established CPU-based method we employ an MPI-parallel implementation of suitable
  quadrature rules.

  * The system of linear equations has to be solved.
  We use Krylov subspace methods, since these methods only
  need matrix-vector multiplications, which can easily
  provided in the established method by the MPI-based distributed representation of
  the matrix.







§ GPGPU QUADRATURE




Even on modern processors of parallel computers it is time-consuming to compute the matrix entries v_ij  for complex industrial geometries. It is advantageous to calculate the surface integral for pairs of collocation points and triangles.

    v_ij   = ∫_∂Ωψ_j()/4π |_i-| dS_y        for all  i,j∈{1,…,n}

Matrix assembly is however ideally suited for parallelization on SIMD-type of processors, because
the integrations require us to
carry out mostly identical operations for all matrix entries.
The currently most common processors are general-purpose
graphics processing units (GPGPUs) like the NVidia Tesla™ or
AMD Radeon Instinct™ cards that contain thousands of
floating-point arithmetic units and offer teraflop-level performance.
Porting the quadratures to GPGPUs poses challenges:

  
  * The most powerful GPGPU models are equipped with fast local memory. We have to ensure that geometrical data is efficiently transferred to the local memory.
  
  * We are using piecewise (mapped) linear basis functions and multiple
triangles contribute to the same matrix entry. Thus we have to avoid
collisions between different triangles that may try to simultaneously update the same matrix entry.

For the parallelization it is beneficial that each collocation point corresponds to one row in the matrix. In our implementation we take advantage of this fact by assigning each collocation point to a thread. Then we consecutively iterate through all triangles of the discretization and simultaneously compute the  contribution of a triangle to all collocation points, see Fig. <ref>. This way collisions can be avoided. For good performance it is required that all threads execute exactly the same operations. This can be ensured by grouping the collocation points depending on whether they are on the surface of a conductor or part of a dielectric interface.

For the quadrature we distinguish between regular, near-singular, and singular integration of the pairs of triangles and collocation points. If the collocation point is a node of the triangle then it is a singular pair. When the distance D between the circumcenter of the triangle and the collocation point is larger than a threshold that scales with the circumradius R of the triangle like

    D > η· R

then it is a regular pair. We used the scaling parameter η = 1.2. All other pairs are near-singular. Circumcenters and circumradii of the triangles can be precomputed.

All three types of pairs can be integrated by using the well-known Duffy-transformation, see <cit.>.
 This is straightforward for regular and singular pairs, only the near-singular pairs need to be integrated adaptively for accuracy. They frequently occur, e.g. in cases with narrow gaps in the geometry, or during the postprocessing, when a point near the surface is to be evaluated. In this near-singular case we first compute the point of the triangle that is closest to the collocation point. Next we subdivide the triangle into smaller triangles
  such that this closest point is a corner-node of a subdividing triangle.
  Then we  again employ the Duffy-transformation to integrate over all smaller triangles.
  This yields an adaptive quadrature with increased accuracy around the closest point. The adaptivity can impact the performance of a GPU-computation badly if no attention is paid, because then there is divergence in the control flow on the GPU. In order to minimize this divergence, we first compute the regular and singular pairs, and deal with the near-singular integrals later.

The categorization into regular, near-singular and singular pairs is carried out on the fly during the iteration through the triangles.
 If a pair is marked as near-singular, then it will be marked as not processed, see Fig. <ref>. They are computed in parallel by using the subdivision method after the regular and singular cases have been completed.  This strategy allows that all three types of integrations are carried out in parallel, without the need to mix operations.
 


 


 The full matrix may not fit in the memory of one GPU for larger problems. Therefore, but also to speed up the computation we use multiple GPUs. Due to the independence of matrix rows we split the matrix into multiple blocks of rows that can be computed and stored independently on different GPUs.







































§ NUMERICAL EXPERIMENTS



In this chapter we show some examples that were computed with the novel GPU-implementation that is based on the H2Lib package, see <cit.>. We first validated our implementation for an axial-symmetric case. We compared the results of the H2Lib with the results of the already existing simulation tools Polopt (3D) see <cit.>,  and Elfi (2D) see  <cit.>. Next we compared the performance of the new GPU-parallel H2Lib implementation with the performance of the existing MPI-parallel Polopt tool.






 §.§ Validation



The benchmark problem that is used to validate the GPU-implementation in H2Lib is a bushing, see Fig. <ref>. It consists of an insulator that is wrapped around a conductor on high-voltage. Five thin conducting sheets are embedded in this insulator. They accomplish the field grading. The outermost is grounded, and the potentials of the other (floating) sheets are unknown and are to be determined. The sheets are treated as single surfaces according to equation (<ref>). Their potentials were computed with all three solvers. POLOPT and H2Lib use the same mesh with 3'526 nodes. The results agree very well, see Table <ref>. The small remaining differences are due to the use of different quadratures.

0.5
      



    ELFI(2D)     POLOPT     H2LIB
Vfoil1     70.8kV     70.15kV     70.22kV 

Vfoil2     51.4kV     50.47kV     50.50kV 

Vfoil3     35.0kV     34.00kV     34.02kV 

Vfoil4     18.9kV     18.02kV     18.02kV 
 
tablePotentials of the conducting sheets of the bushing



0.5


    < g r a p h i c s >

figureCross-section of the bushing








 §.§ GPU-Acceleration



After the successful validation of the H2Lib implementation, we compared the computing times for the GPU-parallel H2Lib and the CPU-parallel Polopt. In both cases we assembled the dense BEM matrix and solved the system with an iterative GMRES with diagonal preconditioner. So the expected numerical work  is quadratically depending on the degrees of freedom, i.e. here the number of nodes (# Nodes).

0.49


    < g r a p h i c s >

figureCumulative times for assembly, solving and  surface electric field computation for POLOPT and H2LIB.



0.49


    < g r a p h i c s >

figureCumulative times for assembly, solving and  surface electric field computation for the H2Lib only.





Figure <ref> shows the times that it took for matrix-assembly, iterative solution, and computation of the electric field at the surfaces of a realistic high-voltage device for different mesh-sizes. POLOPT used  180 CPU cores distributed over 5 nodes with two 18-core CPUs each. The CPU cluster was optimized for these calculations because it is in use by ABB product designers. For H2Lib we used a total of 12 NVidia GTX 1080Ti, distributed over multiple nodes. For POLOPT we clearly recognize the quadratic scaling.
The H2Lib also scales quadratically, see Fig. <ref>, the proportionality constant seems however to be much better than the one for Polopt.


0.49

 #Nodes       #GPUs       H2Lib       POLOPT  
 68'218     2		    18s	    33s 

140'183     8         19s       130s 

232'029     20        26s       371s 

330'706     40        34s  	    702s 

432'084     72	    41s       732s^* 
 



^* 360 Cores


tableCalculation times for different meshes. H2Lib used multiple GTX1080Ti, POLOPT is executed on 180 CPU-Cores.



0.49 

    POLOPT         H2Lib  

     180 Cores      4xTesla P100

Assembly     1r|11s        1r|2s  

Solving     1r|11s        1r|2s  

 surface     1r|11s       1r|2s   
 Total     1r|33s        1r|6s  
Evaluation     102s for      67s for   9219

                17 fieldlines       fieldlines 
 


tableTime for computation and fieldline evaluation. Model size: 68'218 Nodes. Polopt in serial for the evaluation.




We computed another larger example see Table <ref>.  We used 180 cores for POLOPT for all meshes except the largest one, where we used 360 cores. The number of GPUs were chosen such that the matrix  fitted in the combined memory of all GPUs in single precision for H2Lib. Again the GPU-parallel implementation clearly outperforms the CPU-parallel version. We also evaluated the breakdown criterion along the most probable breakdown paths (field lines), see Table <ref>. The acceleration seems even higher, however the evaluation is only implemented in serial in Polopt. 



























§ CONCLUSIONS AND OUTLOOK



The usage of GPUs drastically accelerates the computation of electrostatic fields, as well as their evaluation with respect to breakdown inception.
This opens the door not only for the computation of larger problems, but also for the inclusion of additional physical models (e.g. for surface charging), or for optimization.
Most promising is the combination of the GPU-acceleration with a  compression technique, see <cit.>. Another strong acceleration is to be expected in this case. Moreover the asymptotical behavior will no longer depend quadratically  on the degrees of freedom.



99









VHVlab
A. Blaszczyk, J. Ekeberg, S. Pasncheshnyi, and M. Saxegaard
Virtual High Voltage Lab.
Mathematics in Industry,  SCEE 2016, Springer book.


Polopt1
A. Blaszczyk and H. Steinbigler.
Region oriented charge simulation.
IEEE Trans. Magn., 30(5), 1994.



Polopt2
N. De Kock, M. Mendik, Z. Andjelic, and A. Blaszczyk.
Application of 3D boundary element method in the deisgn of the EHV GIS components.
IEEE Magn. Electr. Insul., 14(3), 1998.



Floating
D. Amann, A. Blaszczyk, G. Of, and O. Steinbach.
Simulation of floating potentials in industrial applications by boundary element methods.
J. of Mathematics in Industry, 4(1), 2014.

DUFFY
M. G. Duffy.
Quadrature Over a Pyramid or Cube of Integrands with a Singularity at a Vertex.
SIAM J. Num. Anal. 19(6), 1260-1262 (1982)



GCA
S. Börm and S. Christophersen.
Approximation of integral operators by Green quadrature and
nested cross approximation.
Num. Math., 133(3), 409–552, 2016.


H2Lib
S. Börm, S. Christophersen et al.
H2Lib, a software library for ℋ- and
  ℋ^2-matrices.
Open source, available at <http://www.h2lib.org>


GCA_SCEE
S. Börm and S. Christophersen
GCA-H^2 matrix compression for electrostatic simulations.
SCEE 2018, Taormina, Italy.






