IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
1
Accelerating Microswimmer Simulations via a Heterogeneous Pipelined Parallel-in-Time Framework Ruixiang Huang, Weifan Liu
arXiv:2604.12083v1 [cs.DC] 13 Apr 2026
This work has been submitted to IEEE for possible publication. Copyright may be transferred without notice, after which this version may no longer be accessible.
Abstract—Simulating large-scale microswimmer dynamics in viscous fluid poses significant challenges due to the coupled high spatial and temporal complexity. Conventional high-performance computing (HPC) methods often address these two dimensions in isolation, leaving a critical gap for synergistic acceleration. This paper introduces a heterogeneous CPU–GPU computing framework specifically optimized for the long-time simulation of filamentous microswimmers in viscous fluid. We propose a two-level parallelization strategy: (1) high-intensity GPU kernels to resolve the quadratic spatial interactions given by the Method of Regularized Stokeslets (MRS), and (2) a distributed MPI-GPU pipelined Parareal architecture to exploit temporal concurrency. By mapping the asynchronous pipeline onto multiple GPU devices, our framework effectively overlaps coarse and fine propagators, overcoming the serial bottlenecks of traditional Parareal method. Furthermore, we employ a GPU-optimized numerical routine for computing the matrix square root arising in the numerical scheme of the filamentous microswimmer simulations. Theoretical analysis of the efficiency improvement of the pipelined Parareal is presented. Numerical experiments demonstrate that the proposed framework achieves order-of-magnitude speedups over CPU-only methods, providing a scalable pathway for simulating complex emergent behaviors in large-scale biology and physics systems. Index Terms—Heterogeneous computing, CPU–GPU architecture, Pipelined Parareal, Biofluid simulation, Flagellar dynamics, Fluid–structure interaction, Parallel-in-time.
I. I NTRODUCTION The locomotion of microorganisms, such as flagellated bacteria and sperm, is fundamentally governed by low-Reynoldsnumber dynamics, specifically Stokes flow. In computational models, these filamentous structures are typically modeled as thin elastic rods. The dynamic interplay between the rod’s elastic properties and the fluid’s viscous resistance forms the fundamental basis for studying biological propulsion and collective swimming behaviors [1]. The large-scale, long-term simulations of such systems are crucial for understanding the mechanisms of self-organization, emergence, and collective motion. However, the underlying fluid-structure interaction (FSI) problems remain extremely computationally expensive. Initially, the Immersed Boundary (IB) method was developed to address such problems by coupling an Eulerian fluid grid with Lagrangian structural representations. However, this This work is supported by National Natural Science Foundation of China (Grant number 12301547) and the Fundamental Research Funds for the Central Universities (Grant number BLX202245). (Corresponding author: Weifan Liu) Ruixiang Huang is with Beijing Forestry University, Beijing, China, 100083 and University of Washington, Seattle, WA, USA 98195. Weifan Liu is with Beijing Forestry University, Beijing, China, 100083
approach often necessitates substantial grid refinement to accurately resolve thin, filamentous structures, leading to a significant increase in computational overhead [2]. As a mesh-free alternative, the Method of Regularized Stokeslets (MRS) was introduced, replacing singular Green’s functions (Stokeslets) with smoothed, regularized kernels [3]. This regularization eliminates the singularity at the point forces, yielding a numerically stable and mathematically consistent framework for computing velocity fields induced by forces and moments distributed along slender bodies. Other notable techniques include boundary integral equation methods [4] and the Rotne-Prager-Yamakawa (RPY) tensor [5, 6], each offering distinct formulations for handling hydrodynamic interactions in slender structures. Despite the elegance and geometric versatility of these Lagrangian-based frameworks, they often rely on evaluating pairwise interactions. Consequently, the computational complexity scales quadratically with the number of discretization points, posing a primary bottleneck for largescale simulations. A secondary, yet equally daunting challenge lies in the temporal discretization. One approach to model the filamentous structures is to employ a version of the Kirchhoff rod [7, 8, 9]. The FSI problem of the Kirchhoff rod immersed in viscous fluid is inherently stiff. To maintain numerical stability, explicit numerical integrators, such as Runge-Kutta schemes, are often constrained by prohibitively small time steps. Consequently, simulating even a single stroke cycle may necessitate millions of iterations, rendering long-term evolutionary studies or large-scale collective dynamic simulations computationally intractable. To break the sequential barrier of time integration, Parallelin-Time (PinT) algorithms, which are a class of methods designed to parallelize the temporal dimension in the numerical solution of time-dependent differential equations, such as those arising in computational fluid dynamics, complex systems and multiscale physics simulations. Parareal, as one notable PinT method, has been proposed to distribute the temporal workload across multiple processing units [10]. By utilizing a coarsegrained predictor to provide an initial guess and a fine-grained propagator to achieve accuracy through iterative corrections, Parareal allows for the simultaneous calculation of multiple partitioned time intervals. In recent years, PinT methods have further evolved into multi-level time-parallel algorithms such as PFASST [11, 12] and MGRIT [13, 14, 15], and have been systematically investigated in fields including fluid dynamics and multi-physics problems [16, 17, 18, 19, 20]. To overcome the prohibitive computational overhead inherent in Lagrangian frameworks, such as the Method of Regularized
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
2
Stokeslets (MRS), where pairwise interactions impose a sigmaintaining the orthonormal frame of Kirchhoff rods. nificant computational load, GPU acceleration has emerged This routine is explicitly optimized for the SIMT (Single as a pivotal strategy. Leveraging the high-throughput data Instruction, Multiple Threads) architecture, ensuring both parallelism of GPUs, researchers have successfully shifted rigorous numerical stability and high-performance execuintensive tasks, including kernel evaluations and state updates, tion of the long-time simulations of flagellar dynamics. from the CPU to the GPU. Recent efforts in FSI have • Comprehensive performance analysis and validation: further demonstrated the efficacy of CPU—GPU heterogeneous We present a theoretical and empirical evaluation of the architectures, underscoring the potential of heterogeneous efficiency of the proposed framework. By analyzing GPU parallelism in tackling complex FSI challenges [21, 22, 23, 24]. idle times and the dependence of speedup on the coarse-toWhile a wide array of HPC algorithms have been developed fine cost ratio and iteration counts, we validate the superifor various fluid applications, specialized frameworks tailored ority of the pipelined multi-GPU approach over standard for the large-scale, long-time simulation of filament dynamics Parareal methods. Numerical experiments involving both remain remarkably scarce. Existing HPC efforts in this domain single and multiple microswimmers demonstrate order-ofhave predominantly focused on either spatial acceleration via magnitude speedups compared to CPU-only methods. GPU-based kernels [25] or temporal parallelism through CPUThe remainder of this paper is organized as follows. based time-integration schemes [20]. However, these decoupled Section II introduces the Kirchhoff rod model for flagellar approaches often fail to address the synergistic computational dynamics, together with the governing fluid equations and demands of biofluid simulations, which necessitate the simulta- numerical methods. Section III describes the implementation of neous resolution of both high spatial and temporal complexities. the proposed space–time parallel framework on a heterogeneous Consequently, there remains a critical need for a specialized, CPU–GPU architecture. Section IV presents the experimental heterogeneous CPU–GPU pipeline capable of concurrently results and their analysis. Section V concludes the paper with managing spatial complexity and temporal concurrency to discussion. enable high-fidelity simulations of motile microswimmers. In this paper, we address these multifaceted challenges by II. P RELIMINARIES introducing a heterogeneous CPU–GPU computing framework specifically optimized for the filamentous dynamics character- A. Kirchhoff rod model istic of sperm and bacterial simulations in biofluid research. To model the structure of the thin filamentous structure Our approach implements a two-level parallelization strategy. of sperm or bacteria flagella, we employ an unconstrained First, we leverage the massive throughput of GPUs to resolve Kirchhoff rod model. The rod is represented by a centerthe spatial complexity of MRS, utilizing custom-designed line X(s, t) and an associated orthonormal material frame kernels for the high-throughput evaluation of linear and angular {D1 (s, t), D2 (s, t), D3 (s, t)}, where s ∈ [0, L] is the Lavelocities. Concurrently, we deploy the Parareal algorithm grangian parameter, initialized as the arclength in the reference across a CPU-based cluster to exploit temporal parallelism. By state. The vectors D1 and D2 lie in the plane of the rod’s synergizing spatial acceleration with time-parallel integration, cross-section, while D3 = ∂X/∂s remains tangent to the this framework mitigates the prohibitive computational cost of centerline in the absence of shear. A schematic illlustration of sequential time-integration for long-time biofluid simulations the space curve representation of a single filament is shown in of collective phenomena and emergent behaviors. the lower left corner of Figure 1. The balance of linear and Our primary contributions are summarized as follows: angular momentum for a rod element in the viscosity-dominated • A scalable multi-GPU pipelined Parareal framework: regime (neglecting inertia) is given by: We design a novel computing architecture that synergizes ∂F ∂N ∂X spatial and temporal parallelism by integrating a pipelined + f = 0, + × F + n = 0, (1) Parareal structure with a distributed MPI-GPU framework. ∂s ∂s ∂s By leveraging the asynchronous scheduling of the pipeline, where F(s) and N(s) denote the internal contact force and our framework enables the concurrent utilization of moment, while f (s) and n(s) represent the external force and multiple GPU devices, effectively overlapping coarse and torque densities exerted by the fluid onto the rod. fine solvers across clusters. This synergy allows the system The constitutive relations for the internal force and moment to handle the O(N 2 ) spatial complexity of Stokesian are expressed in the material frame as F = P3 F i Di and i=1 P3 dynamics while simultaneously exploiting temporal conN = i=1 N i Di : currency, providing a scalable pathway that leverages abundant GPU resources to achieve higher utilization and ∂Dj ∂X i k i i N = ai · D − Ωi , F = bi · D − δ3i , lower GPU idle time in large-scale simulations. ∂s ∂s • GPU-optimized numerical kernels for filament dy(2) namics: We develop high-intensity parallel kernels for the where (i, j, k) is a cyclic permutation of (1, 2, 3). The coeffievaluation of linear and angular velocities, which represent cients a1 , a2 are bending moduli, a3 is the torsional rigidity, the primary computational bottleneck due to the O(N 2 ) b1 , b2 are shear moduli, and b3 is the stretching stiffness. For an pairwise interactions in Stokesian dynamics. Furthermore, axisymmetric rod with a circular cross-section, we set a1 = a2 we employ a numerical routine for calculating the rotation and b1 = b2 . The vector Ω = (Ω1 , Ω2 , Ω3 ) defines the intrinsic matrix square root—a critical yet sensitive operation for curvature and twist of the rod.
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
3
To solve these equations numerically without encountering the singularities associated with point forces (Stokeslets), we To transform the continuous Kirchhoff model into a computaemploy MRS. In this framework, the singular Dirac delta tionally tractable form, the rod centerline is spatially discretized. distribution is replaced by a smooth regularization kernel ϕε (r), A rod of total length L is divided into M points with a uniform known as a blob function. A common choice for the kernel, segment length ∆s = L/(M − 1). The Lagrangian position which we adopt here, is: of each point is defined as sk = (k − 1)∆s for k = 1, . . . , M . Let Xk (t) denote the position vector of the k-th point, and 15ε4 ϕε (r) = , (10) {D1k , D2k , D3k } represent the orthonormal triad defining the 8π(r2 + ε2 )7/2 local material frame at the point. Using a finite difference method, the discrete internal forces Fk+1/2 and moments where r = |x − x0 | and ε is the regularization parameter Nk+1/2 are defined on the segments connecting nodes k and representing the physical radius of the rod. The velocity u and angular velocity ω at a position x induced k + 1. The balance equations for the discrete force density fk by a force fk and torque nk concentrated at Xk are derived and torque density nk at each node are given by: from the regularized Green’s function Gε and the biharmonic function Bε : Fk+ 21 − Fk− 12 fk = − , (3) µu(x) = fk H1 (r) + [(fk · r)r]H2 (r) + (nk × r)H3 (r), ∆s (11) Nk+ 12 − Nk− 12 1 Xk+1 − Xk nk = − − × Fk+ 12 µω(x) = (f × r)H (r) + n H (r) + [(n · r)r]H (r), ∆s 2 ∆s k 3 k 4 k 5 ! (12) Xk − Xk−1 × Fk− 21 , (4) where the scalar functions H1 through H5 are analytical + ∆s expressions derived from the choice of ϕε . The specific j The internal force and moment components Fk+1/2 and formulations of the H1 and H2 are given in [26]. By the principle of linear superposition, the total velocities j Nk+1/2 are computed based on the discrete strains between at any node i are obtained by a matrix-vector product of the adjacent nodes: form N Xk+1 − Xk 1X j j ui f Fk+ 1 = bj · Dk+ 1 − δ3j , (5) = Mij j , (13) 2 2 ∆s ωi nj µ ! j=1 Djk+1 − Djk i · Dkk+ 1 − Ωi , Nk+ (6) where Mij is the mobility matrix. The O(N 2 ) nature of 1 = ai 2 2 ∆s this summation allows for straightforward parallelization, as where (i, j, k) is a cyclic permutation of (1, 2, 3) and the contribution of each node j to node i can be computed (Ω1 , Ω2 , Ω3 ) is the strain twist vector describing how the rod independently. bends and twists. The mid-segment material frame Dik+1/2 is determined based on the rotation between adjacent nodes. D. Parareal Algorithm Specifically, we define the rotation matrix Ak that maps node The Parareal algorithm is a parallel-in-time algorithm. The k to node k + 1: algorithm is designed to solve time-dependent problems by 3 X parallelizing the time domain of the solution. Consider an 1/2 Ak = Djk+1 (Djk )T , Dik+ 1 = Ak Dik , (7) Initial Value Problem (IVP) of the following form: 2 B. Discretization of the Kirchhoff rod model
j=1 1/2
where Ak represents the rotation through half the angle between the two adjacent frames, ensuring a second-order accurate spatial discretization. C. Mathematical model of Stokes flow The locomotion of microorganisms occurs in a viscositydominated regime where the Reynolds number is nearly zero, rendering inertial effects negligible. The fluid dynamics is governed by the incompressible Stokes equations: µ∆u − ∇p + fext = 0,
(8)
∇ · u = 0,
(9)
where µ is the dynamic viscosity, u is the velocity field, p is the pressure, and fext represents the external force density exerted by the microswimmer on the fluid.
dx = u(t, x) for t ∈ (0, T ] and x(0) = x0 . (14) dt Consider a parallel computing environment with Nc cores. In the Parareal algorithm, the time interval [0, T ] is uniformly split into Nc sub-intervals, each of length ∆T = T /Nc , and the time points are denoted as Tn = n∆T for n = 0, 1, 2, . . . , Nc . The algorithm achieves parallelism through an iterative procedure that alternates between two types of computations: sequential calculation using a low-accuracy coarse integrator and concurrent calculation using a high-accuracy fine integrator. The fine solver typically employs a higher-order numerical scheme with a finer temporal resolution, thus providing greater accuracy than its coarse counterpart. We denote the coarse integrator as G(Tn , Tn−1 , xkn−1 ) and the fine integrator as k F(Tn , Tn−1 , Xn−1 ), both of which are used to advance the solution of the initial value problem (IVP).
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
dx = u(t, x), dt
x(Tn−1 ) = xkn−1 ,
4
(15)
for t ∈ (Tn−1 , Tn ], n = 1, 2, · · · , Nc . The procedure of a standard Parareal algorithm for solving Equation (14) is given below: 1) For n = k, k + 1, k + 2, · · · , Nc , compute in serial 0 Xn0 = G(Tn , Tn−1 , Xn−1 ),
X00 = x0 .
(16)
2) For n = k, k + 1, k + 2, · · · , Nc , compute in parallel ′ k−1 Xnk = F Tn , Tn−1 , Xn−1 .
(17)
3) Let Xnk = Xnk−1 for n = 1, 2, · · · , k − 1 and Xkk = ′ ′ Xkk . For n = k + 1, k + 2, · · · , Nc , correct Xnk by applying the coarse solver sequentially as follows:
1) Cost analysis of the pipelined Parareal: We now quantitatively analyze the impact of different scheduling strategies on the overall runtime. When synchronization barriers exist or the wavefront has not fully developed, some GPUs remain idle. For given computational workload and computing capability of GPUs, earlier participation of GPUs in computation generally leads to a shorter completion time. Therefore, we characterize the difference between the pipeline and regular scheduling strategies by analyzing GPU idle time. To describe the idleness, we define the total GPU wait time W , the sum of the individual waiting periods for each GPU, reflecting the extent of resource underutilization due to synchronization overheads or load imbalance. A lower value indicates better parallel efficiency and minimal resource stalling. W =
m X
(GPUi wait time)
(19)
i=1
Assume we partition the time interval into n sub-intervals, and there are m GPUs are available. Let TF and TG denote −G , the computation time of the fine and coarse solver on one (18) partitioned sub-interval respectively. Let T denote the total k−1 where G Tn , Tn−1 , Xn−1 in Eq. (18) has been computed in computation time of the fine solver run in serial, and r > 1 the the previous iteration. time cost ratio of the fine solver to the coarse solver. Then we have TF = Tn , TG = TrF . Let l denote the number of iterations. For the regular scheduling strategy, the coarse phase of the III. GPU- ACCELERATED HETEROGENEOUS k-th iteration needs to sequentially propagate across n − k IMPLEMENTATION time intervals (k = 0, . . . , l − 1). During this stage, only one GPU is active, while the remaining m − 1 GPUs remain idle. A. Pipelined parareal time parallelization Hence, the idle time generated in the k-th iteration can be In the classic Parareal method, the computational time approximated as domain is partitioned into n subintervals, corresponding to (k) Wreg ≈ (m − 1)(n − k)TG . (20) n parallel workers (e.g., CPU cores or processes). Although the fine solver can be executed fully in parallel over these subinSumming over all iterations yields tervals, the coarse solver must still be advanced sequentially l−1 X at the beginning of each iteration, which constitutes a clear l(l − 1) (n − k). = (m − 1)T ln − W ≈ (m − 1)T . reg G G performance bottleneck. To alleviate this issue, we employ a 2 k=0 pipelined scheduling strategy. Let the global time interval be (21) decomposed into N segments, where the p-th subinterval is For the pipelined scheduling strategy, the idle time mainly assigned to the p-th worker. In the regular Parareal scheme, all occurs during the initial injection stage. The parallelism workers must wait until the coarse solver completes the entire increases gradually from 1 GPU to m GPUs, and for n ≫ m, serial prediction from subinterval 0 to n − 1 before launching the system remains nearly fully utilized after the pipeline is the fine solver. In contrast, under the pipelined strategy, worker filled. Hence, the idle time can be approximated by p can immediately start the fine solver as soon as the coarse m−1 X m(m − 1) solver finishes the prediction on the pth subinterval. Wpipe ≈ (m − q)TG = TG . (22) In other words, the coarse solver no longer needs to finish the 2 q=1 entire serial propagation before the fine solver starts; instead, it streams intermediate results downstream as they become By Equation (21) and (22), the difference in the total wait time available. This allows each worker to begin fine solves as soon between the total scheduling strategies is given by as it receives the required local initial condition, eliminating the l(l − 1) m − . (23) ∆W = W −W ≈ (m−1)T ln− reg pipe G inherent serial waiting. Such a strategy is especially beneficial 2 2 on architectures with abundant CPU and GPU resources, as Substituting T = T /r gives G F it significantly reduces GPU idle time and improves hardware T l(l − 1) m utilization. Figure 2 compares the computation flow of the ∆W ≈ (m − 1) ln − − . (24) rn 2 2 regular and pipelined Parareal approaches. Figure 1 presents an overall schematic of the pipeline combining spatial and Based on Equation (24), we analyze the influence of r, m temporal parallelization on a CPU-GPU architecture, with the and n. bottom panel illustrating the GPU arrangement strategy in the (1) Influence of r. The idle-time difference is proportional to pipelined Parareal. 1/r, i.e. ∆W ∝ 1r . A larger r implies faster coarse propagation ′ Xnk = Xnk +G
k Tn , Tn−1 , Xn−1
k−1 Tn , Tn−1 , Xn−1
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
5
Space curve representation of a filament
Fig. 1: Schematic diagram of the proposed parallel computing pipeline
B. Implementation of the solver
Fig. 2: Schematic comparison of the computation flow of a regular Parareal method and a pipelined Parareal method.
and less synchronization-induced idle time. However, the range of r has an upper bound to ensure numerical stability. An excessively large r may lead to loss of stability for stiff problems. l(l−1) (2) Influence of m. Let A := ln − 2 . Then ∆W ≈ T (m−1) rn A− m 2 . Treating m as a continuous variable yields d∆W dm
∝ A − m + 12 , indicating that ∆W reaches its maximum near m∗ ≈ A + 12 . In practice, the effective time-parallelism is limited by the number of time intervals n. When m ≥ n, adding more GPUs no longer increases parallelism and the performance gain saturates. (3) Influence of n. For fixed T and l, ∆W = (m−1) Tr l − l(l−1)+m . For sufficiently large n, ∆W ∼ (m − 1) Tr l, which 2n implies ∆W becomes relatively insensitive to n and is mainly determined by m, l, and r.
We now present the implementation of the solver, organized into three parts. First, we describe the design of the fine and coarse propagators, each consisting of four stages: initialization, strain-twist vector computation for force calculation, velocity computation, and time integration to update positions and orthonormal triads. To ensure data coherency, we adopt a staged kernel execution model. Second, we present a highly parallelized computational pipeline covering each component of the computation, including force and moment evaluation, linear and angular velocity updates, a GPU-optimized matrix square root routine, and time integration. Finally, we detail a thread-mapping strategy designed for the swimmer simulation. 1) Staged kernel execution and data coherency: Following Eq. (18), the Parareal algorithm consists of two components, G and F , which correspond to the coarse and fine propagators, respectively. The coarse solver is typically a low-order scheme used to compute a coarse but fast approximation of the solution. In this work, we employ the explicit Euler method, which is a single-step scheme. The fine solver is typically a highorder numerical scheme that provides a fine and accurate approximation. Here, we present the case where an nthorder Runge—Kutta method is used as the fine solver. In the numerical experiments, we set n = 2, corresponding to the second-order Runge–Kutta method. In GPU-based rod simulations, CPU—GPU data transfer frequently becomes a bottleneck when managed inefficiently. To minimize this overhead, we structure the solver using a staged execution pattern. Each time step is divided into four stages: initialization, force calculation, velocity evaluation, and configuration update. Each stage is implemented as an independent GPU kernel. To eliminate redundant data
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
6
1) Local computation of forces and moments: The discretization of fk and nk in Eqs. (3)–(4) suggests that the force and moment at the k-th point can be written as fk = f (Xk−1 , Xk , Xk+1 ),
(25)
nk = f (Xk−1 , Xk , Xk+1 , Dk−1 , Dk , Dk+1 ),
(26)
which shows that each point depends only on its immediate neighbors. As a result, the computation can be parallelized across discretization points. Algorithm 1 Local evaluation of internal force fk and moment nk 1: for all k = 0, . . . , M − 1 (in parallel) do 2: if 0 < k < M − 1 then 3: ∆Xk+ 12 ← Xk+1 − Xk 4: ∆Xk− 12 ← Xk − Xk−1 P3 5: Ak+ 12 ← i=1 Di,k+1 ⊗ Di,k q 6: Sk+ 21 ← Ak+ 12 h Di,k ← Sk+ 12 Di,k Compute Fk± 12 and Nk± 12 using Eqs(5) and (6). 9: else Apply boundary conditions at k = 0 or k = M − 1 10: 11: end if 12: end for
7: 8:
Fig. 3: CPU–GPU coupled workflow for the nth-Order Runge– Kutta integration solver. The CPU handles global control and host–device data transfer, while the GPU executes the core computational kernels, including strain–twist evaluation, velocity computation, and explicit time integration. Arrows indicate execution flow and data dependencies.
movement, all kernels share a consistent memory layout and operate on the same global memory region. Figure 3 illustrates the flow of each solver. Arrows indicate data dependencies and execution order, while the layout shows how computation and data movement are organized in practice. By launching these kernels sequentially at each time step, our design keeps data resident on the device, reduces launch overhead, and enhances overall performance. Intermediate data are reused directly rather than being copied between kernels. This approach minimizes memory traffic and prevents frequent host–device synchronization. This strategy is depicted in Figure 4, which will be discussed in more detail later. C. High-intensity parallel computational pipeline We now detail the kernel design, covering force and moment calculation, velocity calculation, and time integration. These kernels are implemented with numba.cuda, offering finegrained control over thread and block configurations while supporting asynchronous execution to overlap computation and communication, thereby minimizing latency. Based on the structure of the problem, the computations are decomposed into components with varying levels of parallelism. This hierarchical mapping ensures each task is efficiently executed on the SIMT architecture, maximizing hardware utilization across all simulation scales.
A key difficulty in this step is the matrix square root in Eq. (7). The standard CPU implementation, scipy.linalg.sqrtm, relies on general-purpose algorithms and is not suitable for GPU execution. In our case, however, the matrices involved are 3 × 3 rotation matrices. We exploit this structure by applying a specialized square-root algorithm for this problem [27]. For a rotation matrix R with det R = 1, there exists a unit vector n = (nx , ny , nz )⊤ and an angle θ such that R = I cos θ + (1 − cos θ) nn⊤ + sin θ K(n),
(27)
where
0 K(n) = nz −ny
−nz 0 nx
ny −nx , 0
K(n)2 = nn⊤ − I. (28)
Equation (27) is the Rodrigues rotation formula. Taking the trace yields cos θ = tr(R)−1 . If sin θ ̸= 0, the rotation axis n 2 can be determined from the skew-symmetric part: R − R23 1 32 R13 − R31 . n= (29) 2 sin θ R21 − R12 Let a rotation matrix S satisfy S 2 = R, and assume that it rotates about the same axis n by the half-angle θ/2. Then θ θ θ S = I cos + 1 − cos nn⊤ + sin K(n). (30) 2 2 2 Using K(n)2 = nn⊤ − I and K(n)n = 0, we obtain S 2 = R, which shows that Eq. (30) provides a closed-form expression for the square root of a rotation matrix. In numerical
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
7
implementations, instabilities may arise near the endpoints; we handle them as follows. When θ ≈ 0, we take S ≈ I. When θ ≈ π, sin θ becomes very small, and n should be estimated from the diagonal entries of R: n2x = 12 (R11 + 1),
n2y = 12 (R22 + 1),
n2z = 12 (R33 + 1), (31)
and then set S = R(π/2, n). Compared to Schur-based methods, which involve iterative QR steps and complex control flow, the proposed approach uses fixed-size matrix operations and a predictable sequence of instructions. This makes it more suitable for GPU execution. In practice, we observe that this implementation is both faster and sufficiently accurate compared with scipy.linalg.sqrtm. 2) Computation of linear and angular velocities: Based on Equation (13), the linear and angular velocities can be written in the form: 1 u= Qij fj + Pij nj (32) µ 1 ω = (Wij fj + Rij nj ) (33) µ where Qij (ε, rij ), Pij (ε, rij ), Wij (ε, rij ) and Rij (ε, rij ) are matrices resulting from the regularized kernel functions associated with MRS for rij = xi − xj . Eq. (13) implies that the linear and angular velocities at all points on the rod are independent during evaluation, meaning the summation order can be arbitrarily permuted. This property allows the velocity computation to be parallelized. The GPU-based parallelization is presented in Algorithm 2. 3) Time integration: In the temporal dimension, we update the rod position using the explicit Euler scheme. Accordingly, the position of the k-th point on the rod at time step n + 1 can be written as Xkn+1 = Xkn + u(Xkn )∆t,
(34)
The angular velocity ω nk = ω(Xkn ) is computed using Eq. (13), and the orthonormal triads Di ’s (i = 1, 2, 3) are updated using Rodrigues’ rotation formula: (Dki )n+1 = cos θ (Dki )n + sin θ · e × (Dki )n (35) + (1 − cos θ) e · (Dki )n e, where θ = |ω nk |∆t, e = ω nk /|ω nk | is the unit vector along the rotation axis, “×” denotes the cross product, and “·” denotes the dot product. 4) Thread mapping strategy: The various kernels within the solver exhibit distinct data dependencies and computational patterns, each requiring a tailored thread-block mapping strategy. To address this, we design such a strategy. For local computations, such as internal force and moment evaluation, we map each discretization point to an individual thread and each rod to a dedicated thread block. Since each point depends only on its immediate neighbors, this mapping preserves spatial data locality and enables efficient reuse of rod-level data within high-speed shared memory. In contrast, nonlocal computations, specifically the velocity evaluation based on MRS, employ a target-based mapping
Fig. 4: Thread–block mapping and data flow in the GPU implementation. In Step 2, each thread p in block i processes a subset of source indices in a strided manner, i.e., j = p, p + T, p + 2T, . . ., where T is the number of threads per block. The lower panel illustrates the interaction pattern at a target point, where threads accumulate contributions from multiple rods. strategy. Each target point is assigned to a thread block, within which threads iterate over the corresponding source points. This decomposition exposes the inherent parallelism in source–target point interactions, allowing the GPU to process multiple interactions concurrently while accumulating partial results in shared memory to minimize global memory access. For time integration and other element-wise updates, each discretization point is handled by a single thread. This stage is embarrassingly parallel and requires no synchronization between threads. Collectively, these mapping strategies are meticulously aligned with the mathematical structure of each computational task, ensuring high parallel efficiency and optimal GPU resource utilization. A schematic representation of the thread mapping strategy is shown in Figure 4. IV. N UMERICAL R ESULTS In this section, we apply the proposed framework to the simulations of filamentous microswimmers in viscous fluid. We first verify numerical convergence and report the solution error. We then evaluate the matrix square root algorithm,
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
8
Algorithm 2 Parallel evaluation of u(xi ) and ω(xi ) using shared-memory reduction
wall. Each rod swimmer is initialized as straight and placed in the domain with random orientation and position. Each N rod is discretized into M points along the rod length, with Require: Targets {xi }N ; sources {X , f , n } ; regularj j j j=1 i=1 ization parameter ε; viscosity µ; (optional) wall and LJ ∆s = L/(M − 1). In the following numerical experiments, we set M = 51 and the regularization parameter of the MRS parameters. method to ϵ = 4∆s. To prevent contact between swimmers, a Ensure: u(xi ) and ω(xi ) for all i. on the Lennard-Jones potential U (r) = h force based 1: for all i = 1, . . . , N in parallel do ▷ one CUDA block repulsive 6 i 12 is applied when the pairwise distance r 4ε σr − σr per target point 2: u ← 0, ω ← 0 stored in shared memory between two points falls below a threshold rc = 21/6 σ, with 3: syncthreads σ = 3ϵ. The time step size used in the fine solver is fixed at 4: for all j = 1, . . . , N distributed among threads do ▷ ∆t = 10−6 . each thread processes a subset of sources 5: Compute rij ← xi − Xj , rij = ∥rij ∥ B. Convergence Verification 6: Evaluate kernels Qij (ε, rij ) and Pij (ε, rij ) ▷ To evaluate the accuracy and verify the convergence of the include wall/image correction if needed solver, we examine two metrics: the true relative error η k , LJ 7: (optional) modify fj ← fj + fij if LJ is enabled which measures the error with respect to the true solution, 8: ∆u ← µ1 Qij fj + Pij nj and the relative increment ηek , which quantifies the difference 9: ∆ω ← µ1 (Wij fj + Ri jnj ) ▷ based on between consecutive iterations obtained at the kth Parareal Eqs. (11)–(12) iteration. 10: atomicAdd(u, ∆u) in shared memory xki − xF i 2 11: atomicAdd(ω, ∆ω) in shared memory (37) ηk = max i=1, 2, ··· , Nσ 12: end for xF i 2 13: syncthreads xki − xk−1 i 2 14: Write U (:, i) ← u and W (:, i) ← ω to global memory ηek = max for k = 1, 2, · · · (38) k i=1, 2, ··· , Nσ x i 2 (one thread) 15: end for Here, xki is the position of the ith point at iteration k, xiF is the fully serial fine-solver solution, and N σ is the number of grid points. Because the serial solution is not known in practice, demonstrating the superiority of the GPU-optimized routine. the stopping criterion instead uses the relative increment, i.e. Finally, we present scaling test results, which also validate the the difference between successive Parareal iterates. asymptotic analysis from Section III-A1. We first verify the convergence of the time-parallel algorithm. The time-parallel solution is compared with both the reference A. Environment Setup serial solution and the previous iteration. We solve for the All experiments are implemented in Python. The GPU filament dynamics for t ∈ (0, T ] with T = 1 with number of k platform is equipped with an NVIDIA A100 PCIe 40GB GPU rods 1, 4, 12, 25. Fig. 5 (a) presents the relative increment ηe k running on an aarch64 architecture with Kylin Linux Advanced and Fig. 5 (b) shows the true relative error η obtained in each Server V10. The CPU is a Kunpeng-920 processor at 3.0 GHz simulation. Both quantities are plotted on a semi-logarithmic with 220 GB of available memory. The CPU-only experiments scale. The results show that the error decay steadily and that k are performed on an AMD 7H12 processor (2.6 GHz) with approximately four iterations suffice for a solution accuracy η −12 128 cores and 512 GB memory. Unless otherwise specified, the of order 10 , indicating good convergence of the algorithm. physical parameters remain identical across all experiments. We model a rod-shaped filament swimmer with the Kirchhoff Rods=1 Rods=1 10 4 10 6 Rods=4 Rods=4 10 7 6 rod formulation, a common approach used in computational 10 Rods=12 Rods=12 k 10 8 k Rods=25 Rods=25 studies of flagellar swimmers [7, 26, 28, 29]. The rod’s 10 8 10 9 10 10 10 centerline, is a space curve that defines its geometry, while 10 10 11 its flagellar motion follows a planar sinusoidal wave, modeled 10 12 12 10 0 2 4 6 8 0 2 4 6 8 based on experimental observations [30, 31]. As in [26], the Sweep Sweep waveform is imposed via a prescribed preferred strain–twist (a) (b) vector: Fig. 5: (a) Relative increment ηek (b) True relative error η k (Ω1 , Ω2 , Ω3 ) = 0, −k 2 A sin(ks + f t), 0 , (36) where s denotes arclength, A the amplitude, f the frequency, and k = 2π/λ with λ the wavelength. The rod is situated in a semi-infinite fluid domain {(x, y, z) ∈ R3 | z ≥ 0}, bounded by a stationary, infinite planar wall at z = 0 where the no-slip condition holds. Initially, at t = 0, the rod is straight and parallel to the z-plane, positioned at a height dz = 1 above the
C. GPU Spatial Parallel Performance First, we examine the performance gain of spatial parallelization on a single GPU. We perform experiments for the rod count 1, 4, 12, 25. In each case, we measure the time cost of the three main components that make up each time step’s calculation:
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
9
TABLE I: Runtime comparison of CPU and GPU implementations for different solver components. Rods 1 4 12 25
CPU 0.010812 0.044339 0.132558 0.276209
Initialization GPU Speedup 0.000683 15.83× 0.000680 65.20× 0.000676 196.05× 0.000687 402.03×
Velocity computation CPU GPU Speedup 0.013908 0.000389 35.79× 0.073934 0.000381 193.97× 0.308140 0.000382 805.78× 0.887801 0.000386 2301.74×
the initialization, velocity computation, and orthonormal triads’ update. To evaluate the relative performance improvement of the proposed GPU-based framework to the CPU implementation, we calculate the relative speedup defined as the ratio between the execution time of the CPU implementation and that of the GPU implementation S1 =
TCPU . TGPU
Orthonormal triads’ update CPU GPU Speedup 0.006234 0.000230 27.15× 0.024397 0.000229 106.56× 0.073840 0.000230 321.33× 0.152188 0.000638 238.60×
CPU 0.030954 0.142669 0.514539 1.316198
Total GPU 0.001301 0.001290 0.001288 0.001711
Tcomm Speedup 23.80× 110.60× 399.52× 769.29×
0.00373 0.00376 0.00401 0.22368
TABLE II: Statistical summary of performance and numerical accuracy. Metric Speedup SciPy error CUDA error SciPy time (s) CUDA time (s)
Mean 2.14 2.33 × 10−15 2.76 × 10−15 7.81 × 10−2 3.64 × 10−2
Median 2.10 2.31 × 10−15 4.71 × 10−16 6.63 × 10−2 3.28 × 10−2
Std / Max 0.53 9.68 × 10−16 1.35 × 10−14 2.46 × 10−1 9.92 × 10−2
(39)
We also measure the communication cost, which is the total data transfer time between the host and the GPU device, i.e. Tcomm = TH2D + TD2H , where TH2D and TD2H denote the host-to-device and device-to-host transfer times illustrated in Figure 3, accumulated over a single solver invocation. In Table I, we report the runtime and corresponding speedup of the CPU and GPU implementations for each component. We observe that the CPU runtime increases approximately linearly with the number of rods, with the velocity computation dominating the overall computational cost. In contrast, the GPU runtime remains nearly constant from 1 to 12 rods, indicating that spatial parallelism effectively absorbs the increased workload. When the number of rods increases to 25, the update of orthonormal triads’ shows a slight increase in runtime, but still remains significantly faster than its CPU counterpart. In this case, the overall speedup reaches 769×, demonstrating the effectiveness of spatial GPU parallelization.
increases the variance of the error but does not affect the convergence behavior of the overall algorithm. 0.25
SciPy time 6 CUDA time S2 5
0.20
4
0.15
Time (s)
3 S2
0.10
2
10 13
SciPy error CUDA error
10 14
S2 R F
10 15
1
0.05 0
1
2
3
0
(a)
0.0
0.5
1.0
1.5
2.0
2.5
3.0
(b)
Fig. 6: Performance and accuracy of the 3 × 3 matrix square root implementation: (a) speedup (b) residual
E. Time-Parallel Performance
In this subsection, we perform filament simulations on multiple GPUs and compare the efficiency of the pipelined Parareal method with that of the regular Parareal method. To evaluate the performance and numerical stability of the Simulations are conducted for rod counts of 1, 4, 12, and CUDA implementation of the 3×3 matrix square root algorithm, 25. Let r = TF /TG . We run the simulations on 2 GPUs with we select a random unit vector as the rotation axis and uniformly r = 2, 5, 8, and on 4 GPUs with r = 4, 10, 16. Table III sample 100 rotation angles θ ∈ [−0.1, π + 0.1] to generate reports the runtime comparison between the regular timerotation matrices. We compute the matrix square root of the parallel scheduling and the pipelined scheduling under different rotation matrices using SciPy’s sqrtm and the half-angle parameter settings. Each entry is presented in the format rotation routine described in Section III-C1. regular / pipeline. Under all tested configurations, the pipelined To demonstrate its improvement over SciPy’s sqrtm, We scheduling consistently outperforms the regular scheduling, and calculate the speedup of the GPU-optimized routine relative this advantage holds across simulations involving different rod TSciPy to sqrtm S2 = TCUDA , where TSciPy and TCUDA denote the counts. execution times of of the two implementations respectively. From the numerical results, it is also clear that the perIn Fig. 6 (a), we present the speedup and the numerical formance gap between the two scheduling strategies is more error of each method. The CUDA implementation achieves pronounced when r is small. For instance, when r = 2 (or approximately 2× computation speed across most angles, with r = 4), the pipelined Parareal reduces runtime by approximately an average speedup of 2.14. We present the error ∥S 2 −R∥F in 25%—30% compared to the regular Parareal. As r increases, Fig. 6 (b), and the associated statistical error and time metrics this gap gradually diminishes. This trend is consistent with in Table II. As summarized in Table II, the proposed method the idle-time analysis presented in Section III-A1: when achieves comparable or higher accuracy than SciPy’s sqrtm. the coarse solver becomes faster (i.e., as r increases), the Numerical instability may appear when θ ≈ 0, π, which results coarse wavefront propagates more quickly, reducing the idle in several outliers in the error distribution. This phenomenon waiting time caused by synchronization barriers in the regular D. Performance and Accuracy of the Matrix Square Root Implementation
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
10
TABLE III: Runtime comparison between the regular and pipelined scheduling. Each entry reports regular / pipeline runtime. Rods 1 4 12 25
r=2 815.71 / 552.87 830.77 / 558.27 1420.24 / 1109.08 3461.85 / 2819.40
2 GPUs r=5 508.06 / 408.64 490.23 / 412.75 1111.41 / 938.49 2756.28 / 2381.44
r=8 413.56 / 381.23 431.24 / 389.06 1031.39 / 898.03 2582.23 / 2277.16
r=4 810.22 / 548.23 827.48 / 552.94 970.64 / 729.38 2333.84 / 1749.87
4 GPUs r = 10 501.01 / 420.98 503.83 / 417.04 664.74 / 572.98 1623.55 / 1294.12
r = 16 424.12 / 381.63 424.42 / 379.09 590.75 / 536.67 1447.81 / 1238.45
scheduling. Consequently, the runtime difference between the 1/r gives a slope of approximately 0.9 × 103 ; when m = 4 two scheduling strategies becomes smaller. (with r = 4, 10, 16), the corresponding slope increases to To further validate the theoretical analysis presented in approximately 2.0 × 103 . The increase in slope with respect Section III-A1, we plot the measured runtime gap Treg − Tpipe to m indicates that a larger number of GPUs amplifies the vs. 1/r in Fig. 7. For all rod counts considered in the synchronization-induced waiting cost in the regular scheduling. experiments (1, 4, 12, and 25), the measured data exhibit Although this growth is not strictly proportional to (m − 1) an approximately linear trend, and the fitted lines agree well due to the −m/2 correction and additional non-idle overheads with the experimental points. This behavior provides direct such as communication and memory access, the overall trend empirical evidence for the theoretical prediction obtained from remains consistent with the theoretical prediction, that is, as the GPU wait-time difference estimate ∆W . Since the regular the number of GPUs increases, the performance advantage of and pipeline scheduling strategies execute exactly the same the pipeline scheduling over the regular scheduling becomes sets of coarse and fine computations, their total computational more pronounced. workload is essentially identical. Hence, the dominant source of runtime difference ∆T = Treg −Tpipe is the difference in GPU idle time. Under this assumption, it is reasonable to approximate ∆T ≈ ∆W . The experimental observation Treg − Tpipe ∝ 1r is therefore consistent with the theoretical prediction ∆W ∝ 1r . This agreement indicates that the measured total runtime gap is indeed primarily determined by the difference in idle waiting time, and thus provides a validation of the analysis of ∆W in F. Weak Scaling GPU idle-time analysis. Based on the analysis, this advantage becomes more pronounced for longer simulations, i.e. larger T in equation (24). To evaluate the scalability of the proposed framework, we conduct weak scaling experiments by increasing the time interval length T alongside the number of GPUs, keeping the workload per GPU approximately constant. We show the weak scaling plot in Figure 8 (a) and the corresponding numerical results in Table IV. As T increases from 0.5 to 4 and the number of GPUs increases from 1 to 8, the total runtime grows modestly, increasing much more slowly than the problem (a) (b) size. This demonstrates solid weak scaling performance of the proposed time-parallel framework. We note that the number of iterations l increases slightly from 3 to 5 as the problem size increases, which is expected since a longer time interval requires more iterations for the Parareal algorithm to converge (for a solution accuracy η k < 10−11 ). In all configurations, the pipelined scheduling consistently achieves shorter runtime (c) (d) than the regular scheduling. As the number of GPUs increases, Fig. 7: Runtime gap Treg − Tpipe versus 1/r for (a) 1 rod (b) synchronization overhead in the regular scheme becomes more 4 rods (c) 12 rods (d) 25 rods. significant, whereas the pipelined scheduling maintains higher resource utilization. In the regular Parareal scheduling, the additional idle time mainly arises from the (m − 1) GPUs waiting during the sequential coarse propagation. Hence, ∆T has a dominant term TABLE IV: Weak scaling results that scales with (m − 1), with an additional term proportional T GPU number m Regular time (s) Pipeline time (s) l to −m/2, which accounts for the finite cost of establishing 0.5 1 4832.14 4283.23 3 parallelism in the pipeline strategy. The experimental results 1 2 5122.68 4539.43 3 2 4 7439.62 6416.24 4 reflect the same trend. For instance, in the case of 25 rods, 4 8 11745.27 9231.66 5 when m = 2 (with r = 2, 5, 8), the linear fit of ∆T versus
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
12000 11000 10000 9000 Time (s) 8000 7000 6000 5000 4000 1
Regular Pipeline
2
3
4
5
6
Number of GPUs
7
8
8 7 6 Sp 5 4 3 2 1
11
Sp Ideal Sp
1
2
(a)
4
Number of GPUs
8
(b)
Fig. 8: Scalability results of the proposed solver: (a) Weak scaling performance as the problem size and GPU number increase proportionally. (b) Strong scaling performance with fixed problem size while increasing the number of GPUs.
G. Strong Scaling To evaluate the strong scalability of the framework, strong scaling experiments are conducted by fixing the problem size and increasing the number of GPUs from 1 to 8. We show the strong scaling plot in Figure 8 (b) and the corresponding numerical results in Table V. The speedup and parallel efficiency are defined as Sp =
T1 , Tp
Ep =
T1 Sp = , p p Tp
(40)
where T1 and Tp denote the execution times using 1 and p GPUs, respectively, and p is the number of GPUs. When the number of GPUs increases from 1 to 4, the algorithm achieves near-linear speedup with parallel efficiency above 96%. When scaling to 8 GPUs, the efficiency decreases to 77.1%, mainly due to increased cross-node communication and synchronization overhead. Overall, the proposed time–space hybrid parallel framework demonstrates good scalability and computational efficiency on both single-node and multi-node GPU systems. TABLE V: Multi-GPU strong scaling test GPU number 1 2 4 8
total time (s) 9597.82 4834.20 2492.66 1555.91
speedup 1.99 3.85 6.17
parallel efficient 100% 99.5% 96.3% 77.1%
V. D ISCUSSION Building upon the experimental results presented above, we now discuss several key factors that determine the performance of the proposed framework, including the efficiency of spatial parallelism, the effectiveness of the pipelined temporal scheme, and the trade-offs associated with algorithmic parameters. The strong performance of GPU-based spatial parallelism is primarily due to the tailored mapping strategy between points and rods, along with the use of shared memory. In the current implementation, a system with 25 rods requires only about 2 GB of GPU memory, demonstrating good memory efficiency. Compared to the regular parallel-in-time method, the proposed pipeline approach significantly reduces GPU idle time. This advantage becomes more pronounced in longer simulations (i.e.,
larger T in equation (24)), where waiting time accumulates in the regular scheme. We find that the parameter r plays a critical role in both performance and convergence. A larger r leads to a faster coarse solver and can improve efficiency, but may degrade convergence as it yields solvers of lower accuracy. In particular, when rods are initially close to one another, the system becomes stiffer, and selecting a large r may result in nonconvergence. This highlights the need to balance computational efficiency with numerical stability in practice. In our framework, the number of time intervals n (determined by the number of CPU cores) has limited impact on overall performance, as most of the computational workload is handled by the GPU. Additionally, the separation of computation and data transfer reduces communication overhead, contributing to strong scalability in practice. Despite the advantages, the proposed framework has certain limitations. As the simulation time T increases and more GPUs are employed, resource underutilization may occur, with some GPUs remaining idle during execution. This issue primarily stems from dependencies among time-parallel tasks and load imbalance. More efficient mapping and scheduling strategies between GPUs and solver tasks are therefore needed to further enhance performance. Additionally, the current Python-based implementation limits peak performance. Moving to lower-level languages could yield further speedups. Overall, the experimental results show that a GPU-dominated framework leveraging both temporal and spatial parallelism can effectively minimize idle time while maximizing efficiency, suggesting a robust and scalable solution for space-time parallel simulations in heterogeneous computing environments. VI. C ONCLUSION In this work, we have presented a space–time parallel framework for simulating a fluid-structure interaction problem of filamentous swimmers on a heterogeneous CPU—GPU architecture. The proposed method integrates spatial parallelism on the GPU with temporal parallelism based on the Parareal algorithm. A key contribution is the use of a pipeline structure to mitigate idle time in parallel-in-time execution, coupled with a GPU-friendly implementation of the underlying solver. Through a tailored mapping strategy and optimized memory usage, the proposed approach achieves high-intensity parallelization and efficient utilization of GPU resources. Numerical results demonstrate that the pipeline method consistently outperforms the regular parallel-in-time scheme, particularly in long-time simulations. The method also exhibits favorable scalability and memory efficiency. These findings suggest that the synergistic combination of temporal and spatial parallelism offers an effective pathway to performance improvement for the long-time microswimmer dynamics simulations in biofluid research. Future work will focus on refining the GPU scheduling strategy and enhancing the robustness of the framework for stiff problems. VII. ACKNOWLEDGEMENTS The 3D effect of the rod swimmer in the schematic illustration in the bottom panel of Figure 4 was enhanced
IEEE TRANSACTIONS ON PARALLEL AND DISTRIBUTED SYSTEMS
using Gemini. R EFERENCES [1] C. S. Peskin, “The immersed boundary method,” Acta Numerica, vol. 11, p. 479–517, Jan. 2002. [2] A.-K. Tornberg and M. J. Shelley, “Simulating the dynamics and interactions of flexible fibers in stokes flows,” Journal of Computational Physics, vol. 196, no. 1, pp. 8–40, 2004. [3] R. Cortez, “The method of regularized Stokeslets,” SIAM. J. Sci. Comput., vol. 23, no. 4, pp. 1204–1225, 2001. [4] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press, Feb. 1992. [5] J. Rotne and S. Prager, “Variational treatment of hydrodynamic interaction in polymers,” The Journal of Chemical Physics, vol. 50, no. 11, pp. 4831–4837, 1969. [6] H. Yamakawa, “Transport properties of polymer chains in dilute solution: Hydrodynamic interaction,” The Journal of Chemical Physics, vol. 53, no. 1, pp. 436–443, 1970. [7] L. Carichino and S. D. Olson, “Emergent threedimensional sperm motility: coupling calcium dynamics and preferred curvature in a kirchhoff rod model,” Mathematical medicine and biology: a journal of the IMA, vol. 36, no. 4, pp. 439–469, 2019. [8] N. Ho, K. Leiderman, and S. Olson, “A three-dimensional model of flagellar swimming in a brinkman fluid,” Journal of Fluid Mechanics, vol. 864, pp. 1088–1124, 2019. [9] S. Lim and C. S. Peskin, “Fluid-mechanical interaction of flexible bacterial flagella by the immersed boundary method,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics, vol. 85, no. 3, p. 036307, 2012. [10] J. Lions, Y. Maday, and G. Turinici, “A parareal in time discretization of pdes. comptes rendus de l’acadé, mie des sciences–series i–mathematics 332 (7), 661–668 (2001).” [11] M. Emmett and M. Minion, “Toward an efficient parallel in time method for partial differential equations,” Communications in Applied Mathematics and Computational Science, vol. 7, no. 1, p. 105–132, Mar. 2012. [12] P. Freese, S. Götschel, T. Lunet, D. Ruprecht, and M. Schreiber, “Parallel performance of shared memory parallel spectral deferred corrections,” arXiv preprint arXiv:2403.20135, 2024. [13] R. D. Falgout, S. Friedhoff, T. V. Kolev, S. P. MacLachlan, and J. B. Schroder, “Parallel time integration with multigrid,” SIAM Journal on Scientific Computing, vol. 36, no. 6, pp. C635–C661, 2014. [14] J. Hahne, B. Southworth, and S. Friedhoff, “Asynchronous truncated multigrid-reduction-in-time (at-mgrit),” arXiv preprint arXiv:2107.09596, 2021. [15] N. Margenberg and T. Richter, “Parallel time-stepping for fluid–structure interaction,” Computer Methods in Applied Mechanics and Engineering, vol. 384, p. 113953, 2021. [16] B. W. Ong, “A review of parallel-in-time algorithms,” 2020. [17] A. L. Blumers, M. Yin, H. Nakajima, Y. Hasegawa, Z. Li, and G. E. Karniadakis, “Multiscale parareal algorithm for
12
long-time mesoscopic simulations of microvascular blood flow in zebrafish,” Computational Mechanics, 2021. [18] M. J. Gander, S.-L. Wu, and T. Zhou, “Time parallelization for hyperbolic and parabolic problems,” Acta Numerica, pp. 1–, 2026, arXiv preprint arXiv:2503.13526. [19] J. G. C. Steinstraesser, P. d. S. Peixoto, and M. Schreiber, “Parallel-in-time integration of the shallow water equations on the rotating sphere using parareal and mgrit,” Journal of Computational Physics, vol. 496, p. 112591, 2024. [20] W. Liu and M. W. Rostami, “Parallel-in-time simulation of biofluids,” Journal of Computational Physics, vol. 464, p. 111366, 2022. [21] A. Eghbal, A. G. Gerber, and E. Aubanel, “Acceleration of unsteady hydrodynamic simulations using the parareal algorithm,” Journal of Computational Science, vol. 19, pp. 57–76, 2017. [22] Y. Zeng, Y. Wang, and H. Yuan, “A stable and efficient semi-implicit coupling method for fluid-structure interaction problems with immersed boundaries in a hybrid cpu-gpu framework,” Journal of Computational Physics, vol. 534, p. 114026, Aug. 2025. [23] W. Xue, H. Wang, and C. J. Roy, “Cpu–gpu heterogeneous code acceleration of a finite volume computational fluid dynamics solver,” Future Generation Computer Systems, vol. 158, pp. 367–377, 2024. [24] Q. Li, R. Li, and Z. Yang, “An incompressible flow solver on a gpu/cpu heterogeneous architecture parallel computing platform,” Theoretical and Applied Mechanics Letters, vol. 13, no. 5, p. 100474, 2023. [25] M. T. Gallagher and D. J. Smith, “Passively parallel regularized stokeslets,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 378, no. 2179, 2020. [26] S. D. Olson, S. Lim, and R. Cortez, “Modeling the dynamics of an elastic rod with intrinsic curvature and twist using a regularized stokes formulation,” Journal of Computational Physics, vol. 238, pp. 169–187, 2013. [27] N. J. Higham, Functions of Matrices: Theory and Computation. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2008. [28] S. Lim, “Dynamics of an open elastic rod with intrinsic curvature and twist in a viscous fluid,” Phys. Fluids, vol. 22, no. 2, p. 024104, 2010. [29] S. D. Olson, “Motion of filaments with planar and helical bending waves in a viscous fluid,” in Biological Fluid Dynamics: Modeling, Computations, and Applications, ser. Contemporary Mathematics, A. T. Layton and S. D. Olson, Eds. AMS, 2014, vol. 628, pp. 109–127. [30] H.-C. Ho and S. S. Suarez, “Hyperactivation of mammalian spermatozoa: function and regulation.” Reproduction, vol. 122 4, pp. 519–26, 2001. [31] D. J. Smith, E. A. Gaffney, H. Gadêlha, N. Kapur, and J. C. Kirkman-Brown, “Bend propagation in the flagella of migrating human sperm, and its modulation by viscosity.” Cell motility and the cytoskeleton, vol. 66 4, pp. 220–36, 2009.