Rasterization-based discontinuous integration on the GPU

Introduction

The integration of a discontinuous function $f$ over a domain $D$ dates all the way back to Riemann [28] and early developments in the late 1800s, see, e.g., [32]. It is an important numerical task as it occurs naturally in various settings, e.g., unfitted finite element methods [22, 9, 30], the extended finite element method (XFEM) [4], partition of unity methods [21], computational mechanics [19], and multiphase flow [25] and is challenging as it requires special treatment to ensure sufficient accuracy. Adaptive quadrature [12, 9, 19, 30], moment fitting [23, 34, 18, 17, 16], methods based on the divergence theorem [8] and methods based on equivalent polynomials [35] are among the most common modern approaches to addressing such problems. Several integration methods for discontinuous integrands were tested and compared in [1].

In this paper, we propose a method based on rasterization in computer graphics [27, 11] for the integration of discontinuous functions. The proposed method is based on the observation that rasterization poses a query which is key also for the computation of a discontinuous integrand, namely seeking to establish positional awareness with regards to a set of geometrical primitives in a continuous domain. Given that rasterization is effectively executed once per frame, graphics pipelines are capable of delivering massive throughputs to meet the demand of high frame rate, high resolution applications such as gaming and virtual reality thanks to the massively parallel architecture of modern GPU cards, see, e.g., [20, 2, 31]. Motivated by this computational capacity, the main aim of this paper is to evaluate the possibility of repurposing the standard graphic pipeline for the integration of discontinuous functions and determine its suitability for numerical applications. To this end, we

• present a GPU-resident pipeline for the integration of an arbitrary discontinuous function $f$ over a domain $D$, where discontinuities are encoded using geometric primitives and processed in the rasterization stage, and present both hardware- and software-based implementations using OpenGL and CUDA, respectively,

• use the finite cell method (FCM) [26, 9, 30] as model problem, where the embedded weak formulation essentially presents an integration problem of discontinuous functions at the interface of the physical domain and the fictitious domain,

• compare the accuracy and performance of the proposed pipelines on the GPU with each other as well as with adaptive quadrature on the CPU using state-of-the-art Nvidia GPU cards and Intel and AMD CPUs through an extensive set of numerical benchmarks.

The remainder of this paper is organized as follows. The GPU-resident integration pipeline is presented in § 2. The finite cell model problem is presented in § 3, and the application of the proposed integration method to the finite cell method is discussed. In § 4, the performance and accuracy of the proposed method for the model problem is investigated using numerical experiments. Finally, we draw some conclusions in § 5.

Rasterization-based integration

Given a continuous domain and a set of continuous geometrical primitives, the rasterization stage in a typical graphics pipeline essentially seeks to determine, efficiently and consistently, the coverage and attributes of every pixel on the raster with respect to each primitive. [27, 11] Observing that the raster can act as a discretization of an integration domain $D$ into $\{K_i\}_{i=1}^{N}$ pixels and that the primitives can be used to encode the intervals of a discontinuous function $f : D \to\mathbb{R}$, we formulate a GPU-resident pipeline for the computation of $I_D[f] := \int_D f(x)\,\mathrm{d}x$. The main idea is to approximate $I_D[f]$ as

where $Q$ is an appropriate integration rule, e.g., a Gauss-Legendre quadrature rule.

We assume that $f$ is piecewise smooth, i.e., $D$ admits a decomposition into a finite collection of subdomains $\{P_i\}_{i=1}^{m}$, where $D = \bigcup_{i=1}^{m} P_i$ and $\mathring{P}_i \cap\mathring{P}_j = \emptyset$ $\forall i \ne j$ such that $f$ is smooth on the interior of each subdomain $P_i$, with discontinuities restricted to the interfaces $\partial P_i \cap D$. $Q_D[f]$ in eq. (1) then becomes $\sum_{i=1}^{m} Q_{P_i}[f]$.

We represent each subdomain $P_i$ using a geometric primitive. The rasterizer essentially provides an efficient decomposition of $D$ into subsets $\{K_I\}_{I \in\mathcal{J}_{P_i}}$ intersecting subdomain $P_i$ such that

It is important to note that since subdomains $P_i$ overlap only at the interfaces, repetitions in index sets $\mathcal{J}_{P_i}$ are limited to those pixels at the interfaces.

Graphics pipelines typically use a mix of programmable and fixed-function stages, see, e.g., [20]; therefore, the first challenge is to identify a suitable insertion point that allows sufficient flexibility to compute $Q_D[f]$. We identify the stage after rasterization as the optimal insertion point to compute $Q_{K_i}(f)$ via a compute shader as it provides full programmability and direct access to the rasterization output. In particular, we allocate one floating-point scalar per pixel and for each subdomain $P_i$, perform rasterization, evaluate $Q_{P_i \cap K_I}(f)$, $I \in\mathcal{J}_{P_i}$ and accumulate it on the corresponding global pixel location. We note that since primitives are processed one at a time, there is no need for atomic add at this stage. The partial sums after all primitives are processed are essentially $\{Q_{K_i}(f)\}_{i=1}^{N}$; therefore, the only remaining stage to arrive at $Q_D[f]$ would be to accumulate the results according to eq. (1). The accumulation is achieved using a two-pass tree reduction, first into workgroup partial sums and subsequently into a total sum. The pipeline is shown in Figure 1.

We focus on the 2D case in this work as a proof of concept. The domain $D$ is assumed here to be the normalized device coordinates (NDC) $[-1,1]^2$ for simplicity of notation. For arbitrary domains, we tessellate $D$ into a set of quadrilateral patches and map each patch into the NDC using the Jacobi transform.

Given that data center GPUs such as Nvidia A100/H100 lack fixed-function rasterization circuitry and in order to compare the performance of optimized hardware rasterizers with compute-based approaches, we implement the pipeline in Figure 1 using both OpenGL, leveraging hardware rasterization as well as CUDA, providing a pure compute implementation. In the CUDA pipeline, we employ the rasterization algorithm presented in [20].

Furthermore, we use multisample anti-aliasing (MSAA), where a sub-grid of sample points is used to evaluate the coverage fraction of partially covered pixels, in order to increase integration accuracy on the interfaces. More specifically, on an interface pixel $K_j$

where $|\cdot|$ denotes the measure, $\beta_{i,j}=\frac{|P_i\cap K_j|}{|K_j|}$ is the coverage fraction of $K_j$ by $P_i$, and $\bar{f}_{i,j}:=\frac{1}{|P_i\cap K_j|}\int_{P_i\cap K_j} f(x)\,\mathrm{d}x$. Let $\widetilde{\beta}_{i,j}$ be the approximate coverage fraction estimated by MSAA. We approximate $\mathcal{I}_{K_j}[f]$ in eq. (3) by approximating $\bar{f}_{i,j}\approx\widetilde{f}_{i,j}$, giving

On partially covered pixels, in the hardware-based OpenGL implementation, we use the midpoint rule. Let $f_i$ denote the integrand expression associated with subdomain $P_i$. We then use $\widetilde{f}_{i,j}=f_i(x_j^c)$, where $x_j^c$ is the midpoint of $K_j$. In the CUDA implementation, we have more flexibility with regards to the choice of MSAA sampling strategy as well as the approximation of $\bar{f}_{i,j}$. In particular, in addition to the midpoint rule, we implement a masked Riemann sum as $\widetilde{f}_{i,j}=\frac{1}{S_{i,j}}\sum_{k=1}^{S_{i,j}}f(s_{i,j}^k)$, where $\{s_{i,j}^k\}_{k=1}^{S_{i,j}}$ are the MSAA sampling points located in $P_i\cap K_j$. Subdomains with no sampling points in the pixel contribute zero.

Model problem

We use the finite cell method as model problem in this paper in order to showcase the proposed integration method for numerical applications in practice. For simplicity, we consider Poisson’s equation, which can be used to represent steady-state heat conduction and whose strong form can be written as

where $u$ is the scalar solution variable, $f$ is the source forcing term, $\Omega$ is the spatial domain whose boundary is split into Dirichlet and Neumann parts, respectively, as $\Gamma=\Gamma_D\cup\Gamma_N$, and $g$ and $h$ are pre-defined functions on the $\Gamma_D$ and $\Gamma_N$, respectively. We note, nevertheless, that the application of the proposed integration method in section 2 to finite cell problems remains, in principle, independent of the choice of the operator, see below for more details.

The solution to eq. (4) is typically approximated using an appropriate numerical method, among which the finite element method (FEM) [36] is one of the most popular and well established. A standard conforming finite element discretization of eq. (4) on a boundary-fitted mesh seeks an approximate solution in a finite-dimensional trial space contained in $H^1(\Omega)$ [36]. Given that the generation of the boundary-conforming tessellation of the domain takes a substantial portion of the solution effort in such methods, a class of numerical methods, including immersed boundary methods, the finite cell method and CutFEM have been developed that attempt to circumvent this bottleneck by embedding the physical domain in a structured background mesh, which is trivial to generate, see, e.g., [26, 9, 5, 30] for more detailed treatments and implementation aspects.

The weak finite cell formulation of eq. (4) is obtained by (i) obtaining a conforming weak form following the standard finite element procedure [36], (ii) embedding the physical domain $\Omega$ in a structured computational grid $\Omega_e$ and recovering the physical domain through the introduction of an indicator function $\alpha(\mathbf{x})$, and (iii) imposing the essential boundary conditions weakly using a number of techniques including the Nitsche’s method [24, 15, 30], Lagrange multipliers [10, 13], penalty methods [3, 5], etc. and can be written as follows:

Find $u \in V_e$ s.t.

where $V_e:=\{v\in H^1(\Omega_e)\}$, $H^1$ is the Sobolev space, the terms on $\Gamma_D$ are weak boundary terms in Nitsche’s method, $\lambda$ is a scalar stabilization parameter, computed to enforce coercivity, see, e.g., [14] and $\alpha(\mathbf{x})$ is an indicator function, which essentially recovers the physical domain, defined as

where $\alpha_0=0$ recovers the weak formulation on the physical domain. We note that in practice, in order to avoid severe numerical ill-conditioning, $\alpha_0$ is chosen as $0<\alpha_0\ll1$.

Let $\mathcal{T}_h:=\{T_i\}_{i=1}^{n_T}$ be a tessellation of the computational domain $\Omega_e$ into $n_T$ compact, connected, Lipschitz sets $T_i$ with non-empty interior and $\mathring{T}_i\cap\mathring{T}_j=\varnothing\ \forall\ i\ne j$. In finite cell applications, $\mathcal{T}_h$ is typically chosen to be a regular grid, where $T_i$ are axis aligned, allowing its generation and manipulation to be optimized. In the context of the proposed integration pipeline, the regularity of $\mathcal{T}$ greatly simplifies the transformation from the physical space to the NDC, essentially to a scaling operation. We use space tree data structures in this work for the discretization of $\Omega_e$. [6, 30]

The introduction of $\alpha$ in eq. (5) effectively turns the volume terms in eq. (5) into an integration problem of discontinuous functions and is a challenging aspect of the finite cell solution, typically computed using adaptive quadrature methods [9, 19, 30], moment fitting [16], etc. More precisely, given that eq. (5) is assembled cell wise, discontinuous terms are restricted to cutcells $\mathcal{T}_{h,c} = \{T \in\mathcal{T}_h \mid T \cap\Gamma\ne\emptyset\}$. Here, we employ the integration pipeline in § 2 for this task.

As a result of the embedded formulation, the physical domain $\Omega$ can be defined through its boundary $\Gamma$, typically realized using a discrete representation such as boundary representation (BREP) or using an implicit representation such as level sets. In this work, we employ the former. Nevertheless, we note that the proposed approach is applicable also to the latter since

• even smooth boundary representations such as level sets, necessarily need to be represented discretely, at least where boundary conditions are present, and, therefore, discretization of the implicit boundary representation is virtually always available in such pipelines, and

• the discrete boundary representation can be refined to become arbitrarily close to the smooth curve.

Furthermore, adaptive mesh refinement (AMR) towards the physical boundary reduces the local cell size and can be combined with refinement of the boundary representation. We note that this refinement is also beneficial for the quality of the approximate solution as the boundary layer in many problems includes interesting, rapidly varying information. This refinement strategy is illustrated in § 4.

Once a discrete representation of the boundary $\Gamma_h$ is available, the input data to the proposed integration pipeline, see Figure 1, for each cutcell $T \in\mathcal{T}_{h,c}$ in 2D is simply the polygon consisting of the set of vertices $\Gamma_h \cap T$ and those vertices of $T$ that form an enclosure of $\Gamma_h \cap T$.

We form and transfer the polygons for all cutcells at once from host to device in a preprocessing step, which amounts to $3 \times|\mathcal{T}_{h,c}| \times k \times\operatorname{sizeof}(\mathrm{float})$ bytes of data, where $k$ is the average number of vertices of the polygons per each cutcell. The factor 3 comes from the fact that even in 2D, three coordinates have to be passed to the rasterization pipeline. In implicit representations of $\Gamma$, $k$ can be chosen arbitrarily. In either case, in our experience, especially in combination with the AMR approach as described above, the polygon consists of only a handful of vertices in practice. In addition, we overlap this initial host-to-device communication as well as the following execution of the integration pipeline on the GPU with host computations as all cells can effectively be integrated independently. The output of the integration pipeline is one scalar per function which amounts to $|\mathcal{T}_{h,c}| \times(n_{\mathrm{DoF},T}^{2} + n_{\mathrm{DoF},T})$, where $n_{\mathrm{DoF},T}$ is the number of degrees of freedom per cell.

We note that a change in the differential operator leads to changes in the particular terms in the weak form, and the choice of the basis functions $v$ determines the appropriate quadrature rule; nevertheless, the application of the proposed integration pipeline should remain otherwise intact with regards to these variations. We also note that since the indicator function divides each cutcell into a physical and a fictitious subdomain, and in general wherever the discontinuous function is divided into two mutually exclusive intervals, it is possible to produce an optimized integration pipeline, in which only one of the two subdomains has to pass through the rasterization stage; nevertheless, we explicitly define both the physical and fictitious subdomains through geometric primitives in order to investigate the behavior of the proposed integration pipeline in the most general case.

Numerical experiments

In this section, we investigate the accuracy and computational performance of the proposed integration pipeline in § 2 and its application to the model problem in § 3. For the isolated integration benchmarks in § 4.1 and § 4.2, we consider a dataset consisting of approximately 20 000 straight cuts intersecting adjacent edges on the standard domain $D = [-1,1]^2$. For each cut configuration $c$, let $P_c$ denote the physical inside subdomain and $D \setminus P_c$ its fictitious complement. The dataset is generated such that the physical areas $|P_c|$ are drawn from a truncated normal distribution with parameters $\mu= 1$ and $\sigma= 0.55$, resulting in sampled physical areas $|P_c| \in[5.6 \times10^{-5}, 1.99]$, while the cell area is 4.

Given that piecewise polynomials form a common choice not only for the basis functions $v$ in eq. (5) in the finite cell method but are also more broadly relevant to many other numerical approximation methods such as the finite element method, discrete Galerkin methods, etc., we consider the set of polynomial bases $\mathcal{B} = \{1, x, y, xy, x^2, y^2\}$, as the integration test functions and employ a $2 \times2$ Gauss-Legendre tensor-product rule on fully covered pixels, which is exact for these functions in exact arithmetic. For each configuration $c$ and polynomial $p \in\mathcal{B}$, the benchmark integral is

where

Straight-line interfaces provide a controlled setting for the evaluation of the proposed integration pipeline without an additional geometrical approximation errors as they produce polygonal subdomains and represent cutcells intersected by a single segment of the piecewise-linear boundary representation discussed in § 3. The polynomial test functions on such cut configurations admit analytical reference integrals, which facilitates error evaluation. Furthermore, the polynomial test functions are integrated exactly on fully covered pixels using the employed quadrature rule. These tests therefore isolate the treatment of partially covered pixels and errors arising from finite-precision computation. Reference values $I^{\mathrm{ref}}_{c,p}$ are obtained from analytical polynomial moments evaluated in double precision. The moments over $P_c$ are computed using Green’s Theorem, and the complementary moments are obtained by subtracting these from the corresponding moments on the full cell. Analytical evaluation provides reference values for this controlled verification problem, while the performance comparisons assess the proposed pipelines relative to adaptive quadrature. For a numerical approximation $Q_{c,p}$, the relative error is

The smallest reference magnitude in the dataset is approximately $5.6 \times10^{-5}$. In Figure 2, errors are reported separately for each polynomial as $\operatorname{median}_{c} e_{c,p}$. For the comparison in Figure 6, the median across all cutcell-polynomial pairs is considered.

We present both a hardware pipeline based on OpenGL, denoted by HW, as well as a software pipeline based on CUDA, denoted by SW. We note that due to the limited access modes to MSAA points in OpenGL, see § 2, we only employ the midpoint approximation for partially covered pixels in the hardware pipeline, whereas in the software pipeline both the midpoint rule as well as the masked Riemann approximation are employed. Moreover, in the hardware approximation, the largest MSAA sample size reliably allowed on the tested hardware was 16 points. In the software pipeline, on the other hand, we examine two MSAA sample sizes of 16 and 256 points. Finally, we note that the arrangement of the MSAA sampling points is typically determined by the vendor/API and takes a rotated grid pattern, while in the software pipeline, the sampling points are arranged in a uniform, axis-aligned grid. We note that in the following measurements, GPU completion is ensured before stopping the timer.

We first analyze the convergence and accuracy of the integration pipeline in § 4.1. We then turn our attention to the performance of the pipeline in § 4.2. Finally, we compare the performance of the proposed integration pipeline on the GPU with adaptive quadrature integration on the CPU and furthermore compare their performance for the solution of a benchmark problem when integrated within the finite cell pipeline in § 4.3.

Convergence

The median relative error for each test polynomial as a function of the raster resolution $N$ is shown in Figure 2, where it can be seen that, as expected, the relative error generally decreases as the raster is refined. Nevertheless, several observations are worth a closer look. Namely, it can be seen that while the error in the hardware pipeline in Figure 2(a) steadily decreases throughout the

tested range, the accuracy of the software pipeline in Figure 2(b) and (d) seems to improve consistently only up until $N = 512$, after which the trends depend on the polynomial. To understand this behavior, it is necessary to consider a few interplaying factors. For the test polynomials, the Gauss-Legendre rule on fully covered pixels is exact in exact arithmetic. Therefore, quadrature error arises on interface pixels, where both the coverage fraction $\beta_{i,j}$ as well as the mean integrand value $\bar{f}_{i,j}$ are approximated. For small and moderate raster resolutions, the observed decrease in median relative error with $N$ is consistent with the interface approximations being the main source of error. At finer raster resolutions, the number of pixels grows like $N^2$, and rounding in the evaluation and summation of their contributions can become significant relative to the decreasing interface quadrature error. On the other hand, the coverage estimate also depends on the position of the cut relative to the MSAA samples, and the effect of this coverage error on the integral depends on the polynomial. A similar behavior may also be seen at coarser raster size when the MSAA sampling size is increased, see Figure 2(d) and (e). For example, for the masked Riemann rule, the median error for $p = x^2$ begins to increase beyond $N = 256$ with $16 \times16$ sampling, compared with $N = 512$ with $4 \times4$ sampling. Similarly, increasing the MSAA sample count does not improve every result at the finest resolution, for example, for $p = 1$ at $N = 2048$, the $4 \times4$ grid gives the smaller median error. In addition, rounding in the coverage test may also affect the classification of MSAA samples lying very close to the interface.

For the double-precision software configuration, see Figure 2(f) and (g), the median error decreases monotonically for all functions throughout the tested range, consistent with a delayed onset of round-off domination to finer resolutions. We note, nevertheless, that given the significant performance advantage of single-precision computation on modern GPU cards, the single-precision pipeline may be the more desirable alternative, if the target application admits its accuracy, see § 4.2 for a performance comparison between the single- and double-precision pipelines.

In addition, it can be seen that the hardware and software pipelines provide different convergence curves, even for the same floating point precision and MSAA sampling size. For example, the single-precision midpoint configurations in Figure 2(a) and (b) use 16 MSAA samples per pixel but yield different median relative errors. For $p = 1$, the integrand is constant on each subdomain, so the error reflects the estimation and summation of the covered areas. For this polynomial, the software pipeline gives a lower median error at every tested resolution. Differences in sample placement relative to the cut may contribute to this discrepancy.

Performance analysis

The performance of the proposed integration pipelines in terms of the execution time is broken down in Figure 3. An interesting aspect is the comparison between the hardware and software pipelines. In particular the rasterization stage (Stage 2) is more than two times faster in the hardware implementation for similar configurations, see the first two columns in Figure 3. This aligns well with the theoretical expectation that specialized rasterization circuitry is typically more efficient than general-purpose implementations, see, e.g., [20]; nevertheless, it is equally interesting to observe that the software pipeline achieves a similar total end-to-end runtime as preprocessing (Stage 1) and evaluation (Stage 3) are more efficient in the software pipeline. This warrants the employment of the software pipeline, which is especially important for scientific computing applications as most data-center GPU cards lack support for hardware rasterization, see below.

Another interesting observation is the comparison between single- and double-precision pipelines. On the Nvidia 5090, the double-precision masked Riemann configuration with $4 \times4$ MSAA and a $2048^{2}$ raster takes 0.923 ms per cutcell per polynomial. The evaluations in Stage 3 account for 66.9% of this time. At the same raster resolution and sample count, the single-precision midpoint configuration takes 0.285 ms. The resulting runtime ratio of approximately 3.2 reflects differences in both arithmetic precision and the rule used on partially covered pixels.

The results in Figure 3 are obtained on an Nvidia 5090, which despite its impressive capabilities is from Nvidia’s consumer lineup. What is more often employed in computing cluster ecosystems are data-center cards such as the Nvidia A100 and more recently Nvidia B200. Therefore, we compare the performance of the integration pipeline on different hardware in Figure 4. It is important to reiterate that data-center GPU cards often lack dedicated rasterization circuitry; thus, only the software pipeline is applicable to these cards. Nevertheless, as seen above, the software pipeline is able to achieve a comparable total runtime to the hardware pipeline.

In addition, in the comparison in Figure 4, we have purposefully included a relatively mid-range, consumer-level GPU card, namely the Nvidia 3050Ti in order to showcase the relevance of the proposed method also for usage on standard workstation-type machines.

We present the roofline analysis of the integration kernel on Nvidia V100 next, see Figure 5. The measurements were obtained using Nsight (CUDA 12.1), and the clock was locked at base (SM 1230 MHz, memory 877 MHz) with no throttling. Figure 5 reports the mean performance of the kernel for the integration of a single function on a single cutcell over several captured launches of the kernel and after warmup. We note that some parts of the rasterizer (Stage 2) are not drawn as they execute no or almost no floating-point instructions. As a result, the reported time shares in Figure 5 also do not add up to 100%. Within Stage 2, the bary shader writes a 16-byte value per pixel and sits close to the memory bandwidth, and the fine raster does mostly integer edge-function work, so its FLOP rate is low by construct and does not mean that the kernel is idle. It is also important to note that Stage 2 kernels perform 32-bit operations in both single- and double-precision kernels. Stage 3, i.e., classify and evaluate, is the kernel with the most substantial amount of floating-point instructions and sits close to the ridge, especially in the double-precision kernel.

We note that the benchmarks presented here and in § 4.1 characterize accuracy and runtime for single straight-line interfaces and the selected polynomial integrands. The reported speedups are relative to the adaptive-quadrature implementation and configurations considered in this study. For curved physical boundaries, the accuracy of the overall discretization additionally depends on the boundary representation.

Comparison with adaptive quadrature

In this section, the proposed GPU-resident integration pipeline is compared with adaptive quadrature integration on the CPU using an in-house C++ implementation. The CPU baseline is compiled using the GNU C++ compiler with optimization enabled and run in double precision on a single CPU thread. In adaptive quadrature integration, the integration domain is divided into subcells, where subcells intersected by the cut interface are in turn recursively subdivided up to a refinement level $n_{ai}$. A 2 × 2 Gauss–Legendre rule is applied on the resulting leaf subcells in the presented benchmarks. We use octree data structures, see e.g., [6, 30], to represent the integration domain and perform adaptive mesh refinement (AMR) towards the cut interface, see [30]. Adaptive quadrature remains one of the most common integration methods for discontinuous functions, especially in the context of

unfitted finite element methods, see, e.g., [12, 9, 19, 30].

We employ the same cut configuration dataset as in § 4 and use Intel Xeon E5-2630 and AMD EPYC 9454 CPUs for the benchmark. The accuracy and performance of the adaptive quadrature integration is compared with the proposed integration pipelines in Figure 6. Reported CPU and GPU timings are both end-to-end per polynomial per cutcell. CPU timings include adaptive refinement, integrand evaluation and accumulation of the quadrature contributions. GPU timings include primitive transfer to the device, preprocessing, rasterization, integral evaluation over the raster and reduction of the partial results. The adaptive quadrature integration can be seen to converge monotonically as a function of the adaptive integration levels $n_{\mathrm{ai}}$.

For the runtime comparisons, the SW single-precision pipeline with $4 \times4$ MSAA sampling and a $1024^2$ raster, the SW double-precision pipeline with $16 \times16$ MSAA sampling and a $2048^2$ raster and the HW pipeline with $4 \times4$ MSAA sampling and a $1024^2$ are compared with AQ at $n_{\mathrm{ai}} = 10$, 14 and 7, respectively, based on their approximately comparable median relative errors, as seen in Figure 6(a). The single-precision software pipeline achieves an accuracy roughly equivalent to 10-11 levels of adaptive quadrature refinement. The lower refinement level is conservatively used for the comparison.

In terms of runtime, adaptive quadrature integration is vastly more costly than the GPU integration pipelines at approximately comparable median relative errors. More specifically, it can be seen from Figure 6(b) that single- and double-precision software pipelines on the Nvidia 5090 are around 210x to 840x and 300x to 1200x faster than their adaptive quadrature counterparts, respectively. Similarly, the hardware pipeline is around 30x to 100x faster than adaptive quadrature on the Nvidia 5090 and is around 7x to 30x faster on the modest Nvidia 3050Ti. Speedups are calculated as the ratio of the mean AQ runtime to the mean GPU runtime. For each selected GPU configuration, the reported range spans the ratios obtained with the two tested CPU models, while keeping the GPU runtime fixed. These comparisons quantify acceleration relative to the tested single-threaded AQ implementation. The proposed integration pipelines are therefore expected to lead to massive computational savings, especially when the integration task has to be repeated many times, which is precisely the case in finite cell applications.

In order to showcase the application of the proposed integration pipeline to finite cell problems, we study a finite cell benchmark next. The benchmark consists of the star-shaped domain shown in Figure 7, where the steady-state heat distribution is sought according to the formulation in § 3 under the following conditions: The left half of the boundary is kept in an ice bath, i.e., $u = 0^{\circ}\mathrm{C}$ Dirichlet boundary condition, a homogeneous Neumann boundary condition is applied to the remainder of the boundary, and a constant force $f = 1^{\circ}\mathrm{C}/m^{2}$ is applied to the domain. The discretization is obtained through a combination of uniform refinement and adaptive refinement towards the physical boundary, as discussed in § 3, also illustrated in Figure 7. First-order Lagrange basis functions are used for both trial and test spaces. $\alpha_{0}$ is chosen as $10^{-10}$. The end-to-end solution runtimes are summarized in Table 1.

Disc.	$n_{\mathrm{DoF}}$	$\lvert\tau_h\rvert$	$\lvert\tau_{h,c}\rvert$	AQ(10)	SW(SP)	$\delta_h^{(\mathrm{SP})}$	AQ(14)	SW(DP)	$\delta_h^{(\mathrm{DP})}$
$\tau_h^1$	18463	22939	4702	54	6.4	$2.7\mathrm{e}{-5}$	807.13	12.98	$1.4\mathrm{e}{-6}$
$\tau_h^2$	37391	46582	9400	104.34	12.5	$8\mathrm{e}{-6}$	1609.98	25.66	$6.2\mathrm{e}{-6}$
$\tau_h^3$	74864	93439	18810	172.55	18.14	$7.2\mathrm{e}{-6}$	3219.6	44.48	$2.4\mathrm{e}{-5}$

We note that the integration accuracy required in a finite cell simulation depends on the governing equations, discretization and quantity of interest. Strategies for balancing discretization and quadrature errors are discussed, e.g., in [33]. We therefore consider two prescribed integration configurations, pairing SW(SP) with AQ(10) and SW(DP) with AQ(14), based on their approximately comparable median errors, as established in Figure 6. For each discretization, the benchmark quantifies the end-to-end runtime savings and the changes in the temperature field obtained by replacing the cutcell integration scheme while keeping the remaining computational settings unchanged. All other computations are performed on a single thread of the AMD EPYC 9454 CPU. A systematic study of finite cell solution convergence and the selection of the minimum AQ depth needed to meet a prescribed solution accuracy target are outside the scope of the present work.

For each discretization, let $\mathcal{S}_{h}$ denote the set of all nodes of the mesh lying in the physical domain. We evaluate the maximum relative nodal difference

where $q_{\mathrm{SP}}=10$ and $q_{\mathrm{DP}}=14$. The comparison is performed using the solution temperature fields $u_h$ at the same nodal locations within each CPU-GPU pair. This quantity measures agreement between the computed nodal temperatures. $\delta_h^{(m)}$ for each discretization is reported in Table 1. We note that for the examined cases, the largest absolute nodal temperature differences occur near the physical boundary, while differences over most of the interior are substantially smaller. We report the normalized maximum nodal difference so that the comparison retains these localized peaks.

It can be seen that for the discretizations in Table 1, the end-to-end finite cell solution runtime is reduced by roughly a factor of 8.3-9.5 for SW (SP) relative to AQ with $n_{\mathrm{ai}}=10$ and by roughly a factor of 62-72 for SW (DP) relative to AQ with $n_{\mathrm{ai}}=14$. We note that the exact end-to-end runtime reduction depends on several factors, including the discretization, which determines the number of cutcells, and the computational effort spent in other parts of the solution procedure, especially the linear solver which can in turn be dependent on problem size, system matrix conditioning, etc., see, e.g., [7, 30, 29]. In the benchmark in Table 1, a direct solver is used, whose computational cost remains mainly a function of the problem size and, therefore, independent of the cutcell integration scheme.

The efficient solution of larger problems often requires an appropriate parallelization strategy, e.g., distributed-memory implementations, where the computational domain is decomposed and distributed between the participating processes, see [30]. We note that the proposed pipeline is easy to integrate in such parallelization schemes as the integration of each cutcell remains an independent operation and can be done in parallel. Therefore, each process will be responsible for the integration of its local cutcells. This opens up the door for several interesting parallelization and load balancing strategies, e.g., on clusters with a GPU card per node, the integration of the cutcells on all node-local processes may be offloaded to the local GPU card, leading to a hybrid MPI+X parallelization. However, we note that the exploration of such parallelization opportunities and the investigation of the performance and scaling of the proposed integration pipeline in combination with different solvers and problem sizes is outside the scope of the current study and remains an interesting direction for future work.

Conclusions

In this work, we propose a GPU-resident pipeline for the integration of discontinuous functions, where the function intervals are encoded using geometric primitives and processed in the rasterization stage, and the function is integrated and accumulated in compute shader stages. Furthermore, MSAA is used to increase integration accuracy across partially covered pixels. More specifically, we present a midpoint rule as well as a Riemann sum approach for the approximation of the integrand in partially covered pixels. The proposed integration method is implemented both in a hardware pipeline using OpenGL as well as a software pipeline using CUDA. It is shown that despite more efficient rasterization on special hardware circuitry, the end-to-end runtime of the software pipeline remains competitive, warranting its employment on data-center GPUs which often lack hardware rasterization. The performance and convergence of the proposed pipelines are studied using state-of-the-art Nvidia GPU cards and compared to adaptive quadrature on the CPU. Compared to adaptive quadrature integration, at approximately comparable median relative errors, it is shown that the proposed hardware pipeline can be up to roughly 100x faster, while the single- and double-precision software pipelines are up to roughly 840x and 1200x faster, respectively. We note that these speedups are relative to the tested adaptive quadrature implementation for the specified cut configurations and polynomial integrands and the tested hardware. We use the finite cell method as a model problem and analyze the performance of the presented integration pipeline using a numerical benchmark. It is shown that the end-to-end finite cell solution runtime for the presented discretizations is reduced roughly by a factor of 8.3-9.5 using the single-precision software pipeline and by roughly a factor of 62-72 for the double-precision software pipeline relative to adaptive quadrature integration, while the normalized maximum nodal temperature differences are of order $10^{-6}$ to $10^{-5}$. Therefore, the results indicate that the proposed integration pipeline is a promising alternative to established methods such as adaptive quadrature for the efficient integration of discontinuous functions on the GPU.

Acknowledgments Financial support was partially provided by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) in the framework of the Collaborative Research Center SFB 1683 Interaktionsmethoden zur modularen Wiederverwendung von Bestandstragwerken, project number 523205694. This support is gratefully acknowledged.