Writing Custom Kernel Objects

In scientific computing, particularly in integral equation methods and fast multipole methods, kernel functions play a crucial role. This tutorial will guide you through writing custom PDE kernels using the GenericKernel class template. Predefined kernel functions for Laplace and Stokes in 3D can be found in kernel_functions.hpp. The API for the GenericKernel class is documented in generic-kernel.hpp.

1. Define the Micro-Kernel

   The micro-kernel struct must include the following member functions:

   • Name: Returns the kernel’s unique name.

   • FLOPS: Returns the number of floating-point operations for each scalar kernel evaluation.

   • uKerScaleFactor: Returns the scaling factor for the kernel.

   • uKerMatrix: Computes the kernel matrix given a distance vector (and optionally a normal vector for double-layer kernels).

   Example: Laplace single-layer micro-kernel:

   struct Laplace3D_SL_uKer {
     static const std::string& Name() {
       static const std::string name = "Laplace3D-SL";
       return name;
     }
   
     static constexpr Integer FLOPS() {
       return 6;
     }
   
     template <class Real> static constexpr Real uKerScaleFactor() {
       return 1 / (4 * const_pi<Real>());
     }
   
     template <Integer digits, class VecType> static void uKerMatrix(VecType (&u)[1][1], const VecType (&r)[3], const void* ctx_ptr) {
       VecType r2 = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
       VecType rinv = approx_rsqrt<digits>(r2, r2 > VecType::Zero());
       u[0][0] = rinv;
     }
   };
   

   Example: Laplace double-layer micro-kernel:

   struct Laplace3D_DL_uKer {
     static const std::string& Name() {
       static const std::string name = "Laplace3D-DL";
       return name;
     }
   
     static constexpr Integer FLOPS() {
       return 14;
     }
   
     template <class Real> static constexpr Real uKerScaleFactor() {
       return 1 / (4 * const_pi<Real>());
     }
   
     template <Integer digits, class VecType> static void uKerMatrix(VecType (&u)[1][1], const VecType (&r)[3], const VecType (&n)[3], const void* ctx_ptr) {
       VecType r2 = r[0]*r[0] + r[1]*r[1] + r[2]*r[2];
       VecType rinv = approx_rsqrt<digits>(r2, r2 > VecType::Zero());
       VecType rdotn = r[0]*n[0] + r[1]*n[1] + r[2]*n[2];
       VecType rinv3 = rinv * rinv * rinv;
       u[0][0] = rdotn * rinv3;
     }
   };
   

   The uKerMatrix routine for the double-layer kernel has an additional argument for the normal vector at the source particle. The kernel matrix u[SRC_DOF][TRG_DOF], the distance vector r[DIM], and the normal vector n[DIM] are arrays of SIMD (Single Instruction, Multiple Data) vector type (Vec<Real,DIM>) for enhanced performance. The dimensions SRC_DOF and TRG_DOF are the dimensions of the source density and the target potential; DIM is the dimension of the coordinate space. In the above example, for 3D Laplace, DIM=3 and both the density and the potential are scalars (SRC_DOF=1 and TRG_DOF=1).

2. Define the Kernel Object

   The micro-kernel struct is passed as a template parameter to the GenericKernel class to define the new kernel object:

   using Laplace3D_SL = GenericKernel<Laplace3D_SL_uKer>;
   using Laplace3D_DL = GenericKernel<Laplace3D_DL_uKer>;