Multi-GPU Acceleration of Interface Problems in FEniCS, with Application to Cardiac Modeling

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS, WITH APPLICATION TO CARDIAC MODELING

arXiv:2611.03181v1 [math.NA] 2 Nov 2026

EDOARDO CENTOFANTI∗ , BENJAMIN PACHEV† , JOSHUA STEYER‡ , JAMES D. TROTTER§ , AND PIETRO BENEDUSI¶ Abstract. GPU acceleration of finite element assembly is well established for standard formulations, but mixed-dimensional problems, in which volumetric subdomains are coupled via lowerdimensional interfaces, remain unsupported in automated frameworks such as FEniCS. We extend cuDOLFINx, a GPU-enabled extension of the DOLFINx finite element framework, to assemble such interface problems on GPUs while preserving the high-level FEniCS programming model. Combined with GPU-resident Krylov solvers and algebraic multigrid preconditioning, this enables end-to-end GPU execution of the EMI (Extracellular–Membrane–Intracellular) model, an elliptic interface problem from cardiac modeling in which many intracellular subdomains are coupled trough membranes. We compare the GPU implementation against a CPU-based one through strong and weak scaling experiments on regular and realistic 3D cardiac geometries with up to 8 000 cells and 22 million degrees of freedom. In terms of assembly, a single GH200 node is roughly equivalent to two 144-core CPU nodes in strong scaling, and to two to three nodes in weak scaling, depending on geometry. For the linear solver, performance depends on geometry, preconditioning, and problem scale. On the regular geometry, GPU advantages emerge mainly at larger problem sizes and node counts. On the realistic geometry, with the same algebraic multigrid preconditioner for both architectures, the GPU is faster at all tested scales. The implementation is available as open-source software and applies to other problems with interface coupling. Key words. finite element method, GPU computing, high-performance computing, algebraic multigrid, EMI, DOLFINx MSC codes. 65Y05, 65Y10, 65N30, 65N55, 65F10

1. Introduction. Modern GPU architectures now dominate many HPC systems, offering substantially higher throughput and memory bandwidth than multicore CPUs, which makes them attractive for large-scale finite element simulations, where both assembly and linear solves are major costs. While GPU-enabled linear algebra is available in several libraries [23, 28, 21], GPU acceleration of finite element assembly remains considerably more challenging, due to complex mesh data structures, code generation, and memory management. FEniCS [3] has become one of the most widely adopted open-source finite element frameworks thanks to its expressive Python interface and automated code generation, enabling the rapid development of finite element applications from high-level mathematical formulations. Recently, the experimental CUDA extension cuDOLFINx [27] demonstrated that matrix and vector assembly can be efficiently offloaded to GPUs while requiring only minimal modifications to existing FEniCS applications. By preserving the high-level programming model of FEniCS, cuDOLFINx makes GPU acceleration accessible without requiring users to rewrite their codes using GPU-specific programming models. Despite these advances, an important class of finite element applications remains beyond the capabilities of current GPU assembly frameworks: mixed-dimensional problems. Many multiphysics models involve interfaces, where multiple subdomains ∗ Università di Pavia, Italy † University of Massachusetts Amherst, USA ‡ Karlsruhe Institute of Technology, Germany § Simula Research Laboratory, Norway ¶ Università della Svizzera Italiana, Switzerland (benedp@usi.ch)

1

2

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

are coupled through lower-dimensional manifolds. Such problems require the simultaneous assembly of volumetric and interface contributions, posing challenges that are fundamentally different from those encountered in standard finite element formulations. Extending GPU assembly to these interface problems is therefore a necessary step towards broadening the applicability of GPU-accelerated finite element software. A relevant example of a mixed-dimensional problem is the EMI (Extracellular–Membrane–Intracellular) or cell-by-cell model [34, 24], used in computational cardiology to describe electrical propagation at cellular resolution [32]. Unlike homogenized monodomain and bidomain models [17], EMI represents intracellular and extracellular domains explicitly, coupled through membrane interfaces. Advances in micrometer-scale 3D reconstruction now make realistic cellular geometries available [31], but even a few thousand cells lead to O(107 ) degrees of freedom, making both assembly and linear solution major bottlenecks. Scalable solvers have been developed, including algebraic multigrid [8] and domain decomposition methods such as BDDC [19] and GDSW [9], but mostly on regular geometries or structured meshes. Realistic unstructured geometries have so far been considered mainly in computational neuroscience [6], where only the extracellular space and a single intracellular compartment are coupled. The cardiac setting is considerably more challenging: neighboring cells are electrically connected through gap junctions, so a large number of intracellular subdomains must be coupled simultaneously, giving rise to highly complex interface problems. To address this challenge, in this work we extend cuDOLFINx to support GPU assembly of interface problems, enabling for the first time GPU-accelerated assembly of both matrices and right-hand sides for the cardiac EMI model on large realistic unstructured meshes. Rather than introducing a problem-specific implementation, our approach generalizes the cuDOLFINx framework to interface formulations involving coupled volumetric and codimension-one domains while preserving the highlevel FEniCS programming model. Combined with GPU-resident iterative solvers and preconditioners, this enables efficient end-to-end GPU execution of large-scale EMI simulations. We assess the proposed implementation through weak and strong scalability studies on state-of-the-art HPC architectures using both idealized regular geometries and realistic cardiac meshes. We stress that this study goes beyond the cardiac application, since many other problems share a similar structure, where elliptic differential operators are coupled via jumps and flux conditions on interfaces, e.g., in fracture modeling, where these conditions are generalized via Robin–Wentzell interface conditions [15]. This paper is organized as follows: in Section 2 we introduce the cardiac EMI model. In Section 3 we describe the implementation, in terms of GPU acceleration and solution strategy. In Section 4, we present the computational geometries and meshes used for benchmarking. Finally, in Section 5, we report the weak and strong scalability results obtained using state-of-the-art HPC systems. 2. The EMI model. In this section, we describe the cardiac EMI model, also known as cell-by-cell model [6, 34], using the same formulation as in [22]. We consider the domain Ω ∈ R3 , as composed of N + 1 subdomains, corresponding to N cells Ωi , with i = 1, . . . , N , and an extracellular medium Ω0 . Ionic exchanges through cell membranes are responsible for the interactions between cells and the extracellular space. Adjacent cells can be connected through gap junctions, which contain protein channels responsible for the flow and exchange of ionic species. From a mathematical perspective, the EMI model is intrinsically mixed-dimensional: it couples volumetric

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

3

PDEs in the intra- and extracellular subdomains with lower-dimensional interface equations (typically ODEs) posed at the membranes. The EMI model describes the time evolution of the electrical potential ui : Ωi × (0, T ] → R within the subdomain Ωi and across the interfaces Γij = Ωi ∩ Ωj , for i, j = 0, . . . , N , and has the following form:   − div(σi ∇ ui ) = 0 in Ωi ,    ∂vij + Iij (vij , wij ) in Γij , i ̸= j, −ni · σi ∇ ui = nj · σj ∇ uj = CM (2.1) ∂t     ∂wij − R(vij , wij ) = 0 in Γij , i ̸= j. ∂t Here, vij = ui − uj is the potential jump between cell i and cell j, across the common boundary Γij , σi > 0 is the conductivity coefficient corresponding to the ith subdomain, ni is the outward normal to Ωi , and CM > 0 is the specific membrane capacitance per unit area of the membrane surface. We notice that here we consider a common capacitance; a more general model can include a set of Cij , i.e. with the option of varying capacitance depending on the specific interface. The Iij term (possibly non-linear in vij ) represents the current flowing through the cellular membrane because of ionic exchanges. Typically  −1 if i ̸= 0 and j ̸= 0, Rg v (2.2) Iij (v, wij ) = Iion (v, wij ) otherwise, where Rg is the gap junction resistance, and Iion is a specific ionic model, which depends on a set of M auxiliary variables wij : Γij → RM (e.g., gating variables), related to the behavior of the ionic channels in the cellular membranes. In this setting, R is usually a non-linear function describing membrane dynamics and recovery. As a result, the second equation in (2.1) represents the physiological conservation of current fluxes between the intra- and the extracellular spaces, as well as an ODE system modeling ion flow dynamics. Finally, we close the EMI system (2.1) with homogeneous Neumann boundary conditions: for i = 0, . . . , N ni · ∇ ui = 0

in

∂Ω ∩ Ωi ,

in and with an initial condition on membrane voltage vij = vij for t = 0, given the in in initial voltage vij : Γij → R and initial auxiliary variables wij . According to (2.1), the solving potential can be determined up to a constant, so that uniqueness has to be enforced, for example with an additional point-wise Dirichlet condition u0 (xD , t) = 0 for xD ∈ ∂Ω0 , or at the solution strategy level, imposing a constant1 . We discretize system (2.1) in space using the finite element method, while we use an implicit-explicit (IMEX) strategy to discretize in time: first, the auxiliary variables wij are computed using the jump vij from the previous time step with an explicit treatment of R; then, the ionic current Iij is computed. Finally, we solve the EMI equations with respect to the electric potential. In the following sections, for completeness, we report the finite element discretization, following the notation in [5].

2.1. Weak formulation. By integrating the first equation in (2.1) by parts over each Ωi and incorporating the second equation (defined on the membrane) into the 1 Experiments show that these two approaches are essentially equivalent

4

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

resulting interface terms, we obtain the weak formulation for the i-th subdomain. Specifically, for t ∈ (0, T ], we seek ui ∈ H 1 (Ωi ) for i = 0, . . . , N , such that: Z σi ∇ ui · ∇ ϕi dx + Ωi

N Z X



Γij

j=0 j̸=i

 ∂vij + Iij ϕi ds = 0, CM ∂t

∀ϕi ∈ H 1 (Ωi ).

Substituting vij = ui − uj , we can express the previous equation as (2.3) Z Z N Z N Z X X ∂ui ∂uj ϕi ds − ϕi ds = − σi ∇ ui · ∇ ϕi dx + CM CM Iij ϕi ds, ∂t ∂t Ωi Γi j=0 Γij j=0 Γij j̸=i

j̸=i

with Γi = ∂Ωi \∂Ω, for all ϕi ∈ H 1 (Ωi ). 2.2. Time discretization. We uniformly discretize the time domain [0, T ] with Nt intervals, so that ∆t = T /Nt . We describe the time integration of a single time step; in the complete time evolution then follows accordingly. Given the initial conditions vij in and wij , we discretize the last equation in (2.1) with explicit Euler, i.e., we obtaining in wij in Γij . Then, we define Iij = Iij (vij , wij ) and equation (2.3) becomes Z

Z σi ∇ ui ·∇ ϕi dx+

(2.4) τ Ωi

ui ϕi ds− Γi

N Z X

uj ϕi ds =

Γij

j=0 j̸=i

N Z X j=0 j̸=i

in (vij −τ Iij )ϕi ds,

Γij

where τ = ∆t/CM and the left hand side potentials ui and uj are evaluated at t = ∆t. For i = 0, . . . , N , we define the bilinear and linear forms N Z X in (2.5) di (ui , ϕi ) := τ ai (ui , ϕi ) + pi (ui , ϕi ), fi (ϕi ) := (vij − τ Iij )ϕi ds. j=0 j̸=i

Γij

where Z

Z σi ∇ ui · ∇ ϕi dx,

ai (ui , ϕi ) :=

ui ϕi ds −

pi (ui , ϕi ) :=

Ωi

Γi

N Z X j=0 j̸=i

uj ϕi ds.

Γij

Given H1 = H 1 (Ω0 ) × H 1 (Ω1 ) × . . . × H 1 (ΩN ), the EMI single-time-step weak formulation reads: find (u0 , u1 , . . . , uN ) ∈ H1 so that N X

di (ui , ϕi ) =

i=0

N X

fi (ϕi ),

∀(ϕ0 , ϕ1 , . . . , ϕN ) ∈ H1 .

i=0

2.3. Space discretization. For the spatial discretization, we use the finite element method with linear elements (P1) on tetrahedral meshes. For i = 0, . . . , N , let V (Ωi ) be the finite element space of piecewise continuous linear functions ui,h in Ωi , encoded by ni degrees of freedom. The continuous problem becomes: find (u0,h , u1,h , . . . , uN,h ) ∈ V (Ω0 ) × V (Ω1 ) × . . . × V (ΩN ) such that (2.6)

N X i=0

di (ui,h , ϕi,h ) =

N X i=0

fi (ϕi,h )

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

5

is satisfied for all (ϕ0,h , ϕ1,h , . . . , ϕN,h ) ∈ V (Ω0 ) × V (Ω1 ) × . . . × V (ΩN ). Alternatively, in matrix form, Eq. (2.6) may be written as: (2.7)

Ku = f ,

with K =

N X

Ki ,

Ki = τ Ai + Mi ,

i=0

PN where n = i ni is the global number of degrees of freedom (DoFs), Ai , Mi ∈ Rn×n are stiffness and mass matrices corresponding to ai and pi in Eq. (2.5), respectively, while u, f ∈ Rn collect the finite element solution coefficients and the right and side in Eq. (2.6). In this setting, the mass matrices contain only interface contributions. We remark that K ∈ Rn×n is symmetric by construction. We refer to [5] for a more detailed description of the block structure of the matrix K, with an analysis of the interplay between ni , n and the number of degrees of freedom of the membrane. Under uniform refinement, K is shown to be spectrally equivalent to N +1 independent Poisson problems, one for each cell. Finally, we remark that a discontinuous Galerkin (DG) discretization over the entire domain can also resolve interface jumps, as in [13]. This avoids restricted assembly in Ωi , but increases the number of degrees of freedom. Restricted assembly instead captures interface discontinuities while using a continuous Galerkin discretization within each cell, reducing the degrees of freedom compared to a global DG scheme. 3. Implementation. For the finite element discretization, we use the FEniCS framework (DOLFINx v0.9.0) [3], which translates weak forms expressed in the Unified Form Language (UFL) into low-level C code via the FEniCS Form Compiler (FFCx). PETSc [28], providing a parallel implementation of linear solvers. We use the multiphenicsx [2] library to assemble mixed-dimensional variational problems coupled across subdomains, as in [16, 18, 7]. The code is openly available in the cardioEMI2 repository. We first describe how mixed-dimensional weak forms are set up and assembled on CPUs, we then motivate the need for GPU-accelerated assembly and describe how it is realized in cuDOLFINx. Finally, we detail the solution strategy adopted for the resulting linear using the CUDA backend for PETSc. 3.1. Assembly of the weak forms. Given a mesh with different volume tags for different cells (so that cells sharing an interface have distinct tags), we construct auxiliary facet tags to mark membrane facets, along with a dictionary encoding the dependency graph between cells, which we use to guide the block assembly. In terms of implementation, disjoint cells can be grouped under the same tag, or index i in Eq. (2.5), to avoid looping over individual cells (which can be thousands). Tags, on the other hand, are heuristically bounded by ≈ 5 in realistic cardiac geometries.3 From the dictionary, interface contributions for index i can be isolated, so that all Γij terms are identified and encoded in the FEniCS weak forms following Eq. (2.4). The detailed pre-processing pipeline is illustrated and described in Figure 1. 3.2. GPU-accelerated assembly with cuDOLFINx. As GPU-accelerated solvers in PETSc can significantly reduce solution time compared to CPU ones [10], the assembly of right-hand side vectors becomes a non-negligible bottleneck unless it is also accelerated. This motivates the use of cuDOLFINx to offload assembly operations to the GPU. We remark that other (non-linear) EMI-based problems may 2 https://github.com/pietrobe/cardioEMI 3 No theoretical bound exists for the corresponding coloring problem in 3D.

6

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

INPUT: mesh with different volumetric tags V1 , . . . , V4

Tagging facets, so that Γij can be uniquely determined by tag F

Ω0 Ω1

Ω2 −→

Ω3

−→

Python dictionary .pickle V1 : {F12 , F13 , F14 } V2 : {F12 , F23 , F24 } V3 : {F13 , F23 , F34 } V4 : {F14 , F24 , F34 }

Ω4 FEniCS weak forms

V1

V2

V3

V4

F12 F13 F14 F23 F24 F34

Fig. 1. The pre-processing starts with an input mesh where M volumetric tags V1 , . . . , VM , with M ≤ N + 1, are assigned to subdomains. Adjacent cells have different tags (left panel). Given this input, facets are tagged, so that each Γij tag is unique (middle panel). We use the following tagging strategy for facets, Fpq = Fqp = M · min(p, q) + max (p, q), with p, q = 1, . . . , M denoting volumetric tags. Finally, we create a Python dictionary in the pickle format, mapping facet tags to volume tags (right panel), so that FEniCS weak forms can be assembled correctly.

require matrix assembly at each time step, so that assembly can become a bottleneck, see e.g. [4]. As we will demonstrate in Sect. 5, CUDA assembly kernels in cuDOLFINx deliver substantial speedups over CPU implementations. Assembly performance is largely memory-bandwidth bound, and modern GPUs provide significantly higher memory bandwidth than CPUs. For example, on the NVIDIA GH200 used in our experiments, the GPU is equipped with HBM3 memory, providing 8x times higher memory bandwidth than the LPDDR memory available to the CPU, as confirmed by a STREAM benchmark run on our hardware. This architectural difference makes GPU acceleration particularly well-suited for assembly routines, even when accounting for multicore CPU parallelism. Moreover, running both assembly and solve on the GPU minimizes host–device data movement, which could otherwise limit overall performance. At assembly time, DOLFINx iterates over all active mesh entities (cells, facets,...), and compute the local element tensor Kec (or local vector fec ), and scatters these local contributions into the global matrix K (or vector f ) through a degree-of-freedom map (dofmap). Here and in the following, we denote the global matrix and vector of a generic variational problem by K and f , respectively, consistently with the notation already introduced for the EMI system in Eq. (2.7). For a bilinear form on a cell c with local DoF indices l, m = 0, . . . , nloc − 1, the scatter reads   (3.1) K π(c, l), π(c, m) += Kec [l, m], where π(c, l) denotes the dofmap, mapping the l-th local DoF of cell c to its global index. Linear forms are handled analogously. This architecture separates numerical integration from global bookkeeping (the dofmap and the scatter). Internally, cuDOLFINx reuses the C element kernels generated by FFCx, embedding them directly within CUDA assembly kernels. Because GPU execution demands a much higher degree of parallelism than CPU assembly loops, cuDOLFINx reorganizes the data structures of DOLFINx, in particular the dofmap, mesh coordinates, and coefficient arrays, into layouts that are efficient for GPU memory access patterns.

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

7

After this one-time data transfer, all element integrations and scatter operations are executed on the GPU and CPU–GPU data movement is limited to updating coefficient and constant values between time steps. On a single GPU, cuDOLFINx assembly kernels can match the throughput of FEniCS running on several many-core CPU nodes [27], and multi-GPU assembly has also been demonstrated [26]. 3.3. Restricted assembly via modified dofmaps. For the monolithic linear system in (2.7) each contribution Ki = τ Ai + Mi ∈ Rn×n has nonzero entries in the rows corresponding to the degrees of freedom of subdomain Ωi . The stiffness matrix Ai , arising from the volumetric integral ai , couples only DoFs within Ωi and thus contributes exclusively to the diagonal block (i, i). The mass matrix Mi , arising from R the interface form pi (cf. (2.5)), contributes to both the diagonal block, through Γi ui ϕi ds, and to the off-diagonal blocks (i, j), through the coupling inteR grals − Γij uj ϕi ds for each neighboring subdomain j. To assemble this system in FEniCS, each block form is defined on a continuous Galerkin (CG) function space that is formally constructed on the entire mesh but is restricted to the cells and facets belonging to the relevant subdomain or interface. The solution is therefore continuous within each Ωi but discontinuous across the interfaces Γij . On CPUs, this assembly is handled by multiphenicsx, which creates restricted function spaces and uses PETSc index sets together with local submatrices, which are accessed through MatGetLocalSubMatrix and released through MatRestoreLocalSubMatrix, to assemble each block into the monolithic system. For each block (i, j), multiphenicsx extracts a local sub-matrix indexed by the block’s row and column index sets, equips it with a custom ISLocalToGlobalMapping that maps unrestricted DoF indices to restricted positions (with −1 for inactive DoFs), and then invokes the standard DOLFINx assembly routine on the sub-matrix; vector assembly follows an analogous extract–assemble–restore pattern. Porting this strategy to GPUs would require extracting and restoring sub-matrices across the host–device boundary for every block, negating much of the GPUs throughput advantage. Rather than replicate the PETSc sub-matrix extraction and local-to-global remapping mechanism, cuDOLFINx implements restricted assembly by modifying the dofmap itself, so that the standard GPU assembly kernel can be reused without any changes. The key observation is that in the scatter operation (3.1), the only component that decides where a local contribution is placed in the global system is the dofmap π. By constructing a modified dofmap π̃i for each subdomain form, one can direct the contributions from form i to the correct rows and columns of the monolithic system, while the element integration kernel remains identical. More precisely, let Vh be the CG finite element space defined on the full mesh, with n global DoFs and standard dofmap π. For subdomain Ωi , let Ri ⊂ {0, . . . , n−1} denote the restriction set: the indices of the global DoFs whose support intersects Ωi , so that ni = |Ri | is the number of DoFs of V (Ωi ) introduced in Section 2.3. Let ri : Ri → {0, . . . , ni − 1} be a renumbering map that assigns a local block index to each restricted DoF. The monolithic offset for subdomain i is (3.2)

oi =

i−1 X

nk ,

o0 = 0.

k=0

The modified dofmap for form i is then (
oi + ri π(c, l) , (3.3) π̃i (c, l) = −1,

if π(c, l) ∈ Ri , otherwise,

8

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

where −1 is a sentinel value indicating an inactive DoF whose contributions are discarded during the scatter. With π̃i in place, the assembly of form i on the GPU proceeds exactly as in (3.1). For bilinear forms of arity two (i.e., matrix assembly), the row and column function spaces may belong to different subdomains; accordingly, independent modified dofmaps π̃irow and π̃jcol are constructed from the row and column restriction sets, with their respective offsets. The scatter then reads   (3.4) K π̃irow (c, l), π̃jcol (c, m) += Kec [l, m], c∈T, where T denotes the set of mesh entities (cells or facets) active for the form at hand. For vector assembly (arity one), only the test-space dofmap is modified, and the right-hand side is assembled analogously:   (3.5) f π̃i (c, l) += fec [l], c ∈ Ti . No special scatter logic, sub-matrix extraction, or host–device copies are needed: the standard assembly kernel writes directly into the monolithic system at the correct locations, governed entirely by the pre-computed modified dofmap. 3.4. Software design. At the C++ level, the class CUDADofMap stores the modified dofmap as a flat integer array on the GPU. Its constructor receives the original DOLFINx DofMap, an offset oi , a separate ghost offset (for parallel ghost DoFs, which are stored after all owned DoFs in the monolithic vector), and an optional restriction map ri . It “unrolls” the dofmap, expanding any block structure, and applies the offset and restriction in a single pass, producing the device-side array that the assembly kernel indexes into (cf. (3.3)). The class CUDAForm holds a reference to the underlying DOLFINx form and, when set restriction is called, constructs restricted CUDADofMap objects for each function space of the form. At the Python level, the BlockCUDAForm class orchestrates the block assembly: it collects all subdomain forms, computes the cumulative offsets (3.2) by iterating over the function spaces, builds the restriction sets from user-supplied index arrays, and calls set restriction on each CUDAForm. Boundary conditions are similarly remapped through CUDADirichletBC, which applies the same offsets and restrictions to the constrained DoF indices. This design offers several practical advantages. First, the element integration kernels generated by FFCx are reused verbatim: no kernel rewriting or recompilation is needed for restricted assembly. Second, the dofmap modification is a one-time setup cost; once the modified dofmap arrays reside in GPU memory, every subsequent assembly call (including the per-time-step right-hand side assembly) proceeds at full GPU throughput with no additional overhead. Third, no intermediate sub-matrices or sub-vectors are created, and no PETSc scatter or copy operations are required; the monolithic system is populated directly. As a consequence, the GPU-accelerated version of the code has no multiphenicsx dependency; the restricted-assembly functionality is self-contained within cuDOLFINx. 3.5. Solution strategy. To assess the performance and robustness of iterative solvers, we performed a systematic exploration of Krylov methods and preconditioning strategies. Throughout this work, we compare four solver configurations: Conjugate Gradient (CG) and Generalized Minimal Residual (GMRES), each combined with either Algebraic Multigrid (AMG) or Jacobi preconditioning. AMG, implemented through hypre BoomerAMG, was selected as the primary preconditioning strategy because of its optimal complexity for elliptic operators. Its performance was compared

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

9

against a simple Jacobi preconditioner, which is highly parallel and therefore well suited to GPU architectures. To ensure a fair comparison between different solver and preconditioner combinations, all simulations employed a relative tolerance of 10−7 and monitored convergence using the unpreconditioned residual norm. We remark that matrix-free multigrid variants, which are likely better suited to GPUs, have recently been studied for interface problems on regular geometries [20]; however, no matrix-free AMG formulation currently exists, so this approach sacrifices the geometric flexibility needed for complex geometries. A substantial part of this study was devoted to identifying robust AMG parameters for CPU and GPU runs. The main difference between the two architectures lies in the setup phase of AMG, namely in the algorithm that splits the unknowns into coarse (C) and fine (F) points. Unless prescribed otherwise, BoomerAMG uses HMIS coarsening with extended+i interpolation on CPUs, and PMIS coarsening with a reformulation of extended+i interpolation in terms of sparse matrix–matrix products on GPUs [23, 25]. Both algorithms reduce the operator complexity of classical parallel coarsening schemes [12]. PMIS selects the C-points as a maximal independent set of the strength-of-connection graph, breaking ties randomly; it requires only local operations and is currently the only coarsening available on GPUs [14, 23]. HMIS first applies the inherently sequential first pass of the classical Ruge–Stüben algorithm within each process subdomain, and then PMIS across process boundaries [12], so that the resulting hierarchy depends on the parallel partition. Compared with HMIS, PMIS tends to produce coarser hierarchies with lower complexity, at the price of a less effective coarse-grid correction and, hence, of a slower AMG convergence [12, 14, 25]. Depending on the geometry type (regular or realistic), we use different smoothers and coarsening types, motivated by an extensive iper-parameter search. For regular geometries, we used the default coarsening and interpolation of each architecture and the hybrid L1 -scaled SOR/Jacobi relaxation scheme4 . In this case, the difference between the CPU and GPU preconditioners is negligible, and AMG-CG requires essentially the same number of iterations for both architectures. For realistic geometries, whose meshes are unstructured and highly heterogeneous, the quality of the C/F splitting matters more, and the dependence of the AMG preconditioner on the parallel partition becomes relevant. Moreover, the L1 scaled SOR/Jacobi smoother may yield an AMG preconditioner that is not positive definite, which can cause AMG-CG to fail. To compare the two architectures fairly, we therefore prescribe the same, partition-independent AMG configuration on CPUs and GPUs, with PMIS coarsening and extended+i interpolation, a Chebyshev polynomial smoother on all levels, including the coarsest one, and no C/F-ordered relaxation. With these settings, CPU and GPU runs build essentially the same hierarchy , resulting in the same iteration counts across architectures and node counts, cf. Tables 8 and 9. We chose the Chebyshev smoother over the L1 -scaled Jacobi one, which is also independent of the partition, because it was slightly faster per time step. In all cases, we use the default strong threshold of 0.25, and no aggressive coarsening. The complete set of PETSc and hypre options required to reproduce the experiments is reported in Appendix A. 4. Benchmark models. In this section, we outline the model scenarios used to evaluate the performance of the EMI framework in the cardiac setting. These models are designed to bridge the gap between idealized computational benchmarks and 4 This smoother applies Gauss–Seidel within each process subdomain and a L -scaled Jacobi 1

across subdomains. It is equivalent to an inexact block Jacobi using SOR within the blocks.

10

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

Mesh name

Myocytes (N )

Elements

DoFs

A1 A2 A3 A4

20×10×5 20×10×10 20×20×10 20×20×20

8 290 738 16 585 884 33 185 104 66 381 339

2 122 086 4 207 204 8 377 104 16 683 928

Table 1 Model A: number of elements and degrees of freedom (DoFs) for each configuration.

biologically relevant simulations, directly addressing our research questions regarding solver scalability and robustness. We employ two distinct configurations: Model A, which utilizes a modular geometry designed for scalability tests, and Model B, which represents a transition toward real-world use cases. For clarity, we denote the number of degrees of freedom n as “DoFs” throughout the subsequent tables. 4.1. Model A. To evaluate scalability and performance in a controlled environment, we first employ a modular, synthetic geometry. In this model, cardiac cells are represented as regular blocks with smaller cubes attached to each face to represent gap junctions. This highly structured configuration is particularly suitable for performing weak and strong scalability tests; its inherent symmetry allows for an even distribution of the computational load among processors or GPUs, minimizing load imbalance and providing a clear baseline for hardware performance. We show an example of this geometry in Figure 2, while the number of elements and DoFs are reported in Table 1 for each configuration. 4.2. Model B. To assess performance in a biologically relevant scenario, we consider a second realistic geometry representing a 1 mm3 volume of cardiac tissue. Figure 3 illustrates the subdomain tagging, while Figure 4 shows the corresponding solution. This geometry was created with the anni algorithm [30], which creates unstructured meshes starting from imaging data of myocyte alignments. In short, the algorithm first generates a point cloud of myocyte centers and links neighboring centers, splitting each link halfway. The union of the half-link domains associated with a center defines a myocyte. Heuristic constraints, such as a minimum distance between centers or a maximum myocyte length, govern the creation of this network. The network then defines a level-set function describing the cell membranes, which is discretized with the remesher mmg [11]. To facilitate weak scalability tests and investigate solver robustness, the mesh is progressively refined, again with mmg, across four levels, ensuring that surfaces are preserved without changing the topology of the myocytes alignment. The statistics of the different refinement levels are given in Table 2. We remark that this geometry is highly complex, and the meshing process can produce nearly degenerate tetrahedra (here removed), typically where the extracellular space between neighbouring cells has almost vanishing thickness. We therefore stress the importance of careful mesh generation and inspection for physiological simulations. 4.3. Ionic models. For both model A and B, we consider two separate models for the interface current Iij , depending on the type of membrane, following equation (2.2). Among the wide range of available ionic models, we select two models with a favorable trade-off between computational cost and physiological accuracy,

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS Mesh name

B1

DoFs nelems Myocytes (N ) h ρ ρmin Γmyo VECS /VICS

3 391 119 10 799 476 7 843 11.7 ± 5.6 0.420 ± 0.215 1.0 · 10−6 13 800 ± 2 700 0.888

B2

B3

11

B4

6 601 594 10 999 350 22 122 183 25 447 781 44 738 748 94 367 416 7 845 7 845 7 845 8.37 ± 3.65 6.74 ± 2.33 5.15 ± 1.56 0.436 ± 0.263 0.562 ± 0.262 0.664 ± 0.221 1.0 · 10−6 1.0 · 10−6 1.0 · 10−6 15 800 ± 3 200 16 000 ± 3 200 16 100 ± 3 200 0.650 0.632 0.617 Table 2 Model B statistics. If applicable, values are shown as the mean and the standard deviation over all the entities of the mesh. All spatial quantities are given in µm. The variable nelems denotes the number of tetrahedra in the meshes, N the number of myocytes, h the edge length, and ρ the mesh quality, defined as the ratio of inscribed to circumscribed sphere radius (normalized to 1 for a regular tetrahedron); ρmin is its minimum over the mesh. Γmyo denotes the membrane area per myocyte (including cell–cell junctions) and VECS /VICS the ratio between extracellular and intracellular spaces.

Fig. 2. Model A: an example of 3×3×3 myocytes regular geometry. Different colors indicate the subdomain tags. This case presents a total of M = 3 tags (see also Figure 1).

capturing intercellular conduction while remaining inexpensive to evaluate. For gap junctions (i.e., interfaces between cells), we adopt a Plonsey passive linear model [29] with no additional state variables, corresponding to purely ohmic conduction, so that only Rg in (2.2) is set. For the cellular membrane (i.e., interfaces involving the extracellular media), we employ the Aliev–Panfilov (AP) ionic model [1], where Iion is a nonlinear function of the transmembrane potential and a single recovery variable w = w is present: (4.1)

Iion (v, w) = kv(v − a)(v − 1) + vw,

where k > 0 is an excitability parameter and a > 0 a threshold parameter. The evolution of the gating variable w is governed by   µ1 w (4.2) R(v, w) = ϵ + (−w − kv(v − a − 1)) , v + µ2 where ϵ, µ1 , µ2 > 0 are recovery parameters. The AP model is expressed in dimensionless form. In particular, following [1], the normalized variable vij is related to the

12

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

Fig. 3. Model B geometry, with the coarsest setting (B1). Different colors indicate subdomain tags. Six different tags (M = 6) are needed to avoid neighboring myocytes with the same tag. The myocyte network is surrounded by an extracellular bath (transparent box). The magnified circle highlights the tetrahedral elements employed for the intra- and the extracellular space discretization.

Fig. 4. Model B example solution. Left panel: intracellular potential (here denoted by ui ) and extracellular one (ue ) at t = 100 ms, for a traveling wave starting from the left part of the domain. Let us note that the potential wavefront is not homogeneous, because of the non trivial connectivity network between myocardiocytes. Right panel: corresponding membrane potential vij = v01 and recovery variable w evolution at a reference spatial point in Γij .

measured transmembrane potential E (in mV) through the affine transformation (4.3)

vij =

E + 80 , 100

so that vij = 0 corresponds to the resting potential Erest = −80 mV and vij = 1 corresponds to a depolarized state with amplitude 20 mV. Similarly, time is rescaled according to tphys [ms] = αt t, where t denotes the dimensionless time variable used in the model. This scaling ensures that the action potential duration matches physiological values (the action potential duration is ≈ 25.58 units, while in the canine physiological case is ≈ 330 ms [1]). We provide the list of the two models parameters in Table 3.

13

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

Parameter

Symbol

Value

Unit

Common parameters

Membrane capacitance Timestep

CM ∆t

1 0.005

µF/cm2 ms

Passive (gap junction)

Gap junction resistance

Rg

0.003

kΩcm2

Excitability parameter Threshold parameter Recovery parameter Recovery dynamics Recovery dynamics Minimum potential (scaling) Maximum potential (scaling) Time scaling factor

k a ϵ µ1 µ2 vij (min) vij (max) αt

8.0 0.15 0.002 0.2 0.3 -80 20 12.9

– – – – – mV mV ms−1

Aliev–Panfilov

Table 3 Physical and model parameters. The AP model is solved in dimensionless form; vij (min) and vij (max) define the affine transformation between physical voltage (mV) and normalized variables, while αt accounts for time rescaling between model units and milliseconds.

5. Experimental results. In this section, we evaluate the performance of the GPU-accelerated implementation and compare it to the CPU version. We begin with an evaluation of the right-hand side assembly routines, and proceed to investigate the performance of the linear solver, presenting both strong and weak scaling experiments. The experiments for Model A, as well as the assembly experiments for Model B, are run on the Vista system at the Texas Advanced Computing Center (TACC), and the linear solver experiments for Model B on the Alps system “Daint” at the Swiss National Supercomputing Centre (CSCS). Vista provides both GPU and CPU nodes: each GPU node hosts one NVIDIA GH200 Grace Hopper Superchip, while each CPU node contains two 72-core NVIDIA Grace CPUs, for a total of 144 cores. On Daint, each node is equipped with four GH200 superchips5 . To compare CPU and GPU performance consistently across the two systems, we use the following convention throughout the paper: a GPU node denotes a single GH200 driven by one MPI rank, and a CPU node denotes 144 Grace cores, each driven by one MPI rank. We stress that a GPU node and a CPU node are not directly comparable in cost, power consumption, or architecture; the one-to-one node comparison used below is adopted purely for convenience in organizing the results, not as a claim of fairness. In all the experiments, we consider one hundred time steps (Nt = 100), and report the average results per time step. In terms of initial conditions, for all i, j, we consider the initial in membrane potential vij = −80 · (1 − χS ) mV, with S ⊂ Ω a initially activated region in and χS its characteristic function and wij = 0. 5.1. Right-hand side assembly. We perform both strong and weak scaling tests of the right-hand side vector assembly. We remark that matrix assembly is time-independent, so it has to be performed only once, thus it becomes negligible for large Nt . The results are summarized in Table 4. The strong scaling results show that, for assembly, a single GH200 matches the performance of at least two 144-core CPU nodes. For the weak scaling test, one node was used for A1 and B1 meshes, and the number of nodes was doubled for each refined 5 Daint nodes are connected through an HPE Slingshot-11 network.

14

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

Nodes 1 2 4 8

A1 GPU

A1 CPU

B1 GPU

B1 CPU

0.05 0.03 0.02 0.01

0.19 0.09 0.04 0.02

0.28 0.15 0.08 0.04

0.61 0.31 0.15 0.08

Table 4 Strong scaling of the EMI vector assembly routine for the meshes A1 and B1. Performance numbers are reported in seconds per timestep.

Nodes 1 2 4 8

A GPU

A CPU

B GPU

B CPU

40.51 37.63 36.45 33.80

11.02 10.20 9.90 9.31

11.96 9.43 8.87 8.19

5.54 6.31 6.32 6.44

Table 5 Weak scaling of the EMI vector assembly routines. Performance numbers are reported in millions of DoFs per second, normalized by the number of nodes. Larger numbers indicate better performance. Meshes A1 to A4 and B1 to B4 are used when increasing the number of nodes.

mesh. This led to a roughly equal ratio of DoFs to nodes across runs. To facilitate comparison between runs, performance is reported in millions of DoFs per second and normalized by the number of nodes. The results are summarized in Table 5. Unlike the strong scaling test, the weak scaling test shows different performance ratios between the GPU and CPU for the Model A and Model B meshes. For the regular meshes (A1-A4), each GPU node is the equivalent of three CPU nodes, while for the realistic meshes (B1-B4), one GPU is approximately 30% to 100% more efficient than a single CPU node. 5.2. Linear solver. For the regular geometry (Model A), we used AMG preconditioned CG both in the GPU and CPU cases. For GPUs, we also included results using Jacobi preconditioning, since this strategy is more effective than AMG in some cases. For the realistic geometry (Model B), AMG preconditioning is used in all cases, with the same configuration on CPUs and GPUs (further details in Section 3.5). We perform strong and weak scaling tests, with the same mesh progression used in Section 5.1 for experiment on assembly times. 5.2.1. Regular geometry. For Model A, Tables 6 and 7 contain strong and weak scaling results. In terms of strong scaling, we note that Jacobi preconditioning is more effective than AMG for a lower GPU count, despite the much larger number of Krylov iterations. This trend reverses as the number of GPUs nodes increases, since communication scales linearly with the number of iterations. The CPU solver does not scale past two nodes, which can also be explained by communication overhead. Each CPU node contains 144 cores, so the CPU solver utilizes a much larger number of MPI processes than the GPU solver, which requires only one MPI process per node. Single-node GPU performance is worse than single-node CPU performance, but the GPU solver scales better. The weak scaling study in Table 7 confirms the same trend. The communication overhead and the increasing iteration count explains the performance deterioration

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

Nodes 1 2 4 8

GPU Jacobi-CG

GPU AMG-CG

CPU AMG-CG

0.90 (1803) 0.65 (1802) 0.52 (1802) 0.47 (1803)

1.41 (9) 0.98 (10) 0.55 (10) 0.39 (10)

0.39 (10) 0.25 (10) 0.25 (9) 0.44 (10)

15

Table 6 Strong scaling of the EMI linear solve routines for mesh A1. Performance numbers are reported in seconds per timestep. The number of Krylov iterations is reported in parentheses.

Nodes 1 2 4 8

GPU Jacobi-CG

GPU AMG-CG

CPU AMG-CG

2.36 (1803) 1.98 (1854) 1.85 (1917) 1.51 (2265)

1.51 (9) 1.79 (10) 2.06 (10) 2.09 (9)

5.38 (10) 4.04 (10) 2.48 (9) 1.67 (9)

Table 7 Weak scaling study of linear solve routines for a sequence of regular geometry meshes (from A1 to A4 in Table 1). Performance numbers are reported in millions of DoFs per second, normalized by the number of nodes. Larger numbers indicate better performance. Krylov iterations are reported in parentheses for the linear solvers. The runtimes in seconds are plotted in Figure 5.

with more nodes. In terms of weak scaling, the CPU solve is sub-optimal so that the GPU implementation outperforms the CPU one on a per-node basis for the largest run (mesh A4 on 8 nodes). As expected, AMG preconditioning is robust to mesh refinement in all cases. 5.2.2. Realistic geometry. Tables 8 and 9 contain strong and weak scaling results for Model B. Strong scaling saturates early on B1, for both CPU and GPU. On B2 the runtimes keep decreasing up to 8 nodes, although the GPU runs out of memory on a single node. The number of iterations stays essentially constant across node counts in all cases and GPU runs are consistently faster than the corresponding CPU runs. From two to eight nodes, B2 speeds up by a factor of 1.7 on GPUs and 1.5 on CPUs. On the GPU, the cost per Krylov iteration stays between about 9 and 12 ms in all the B1 runs, including the single-GPU one, where no MPI communication is involved, and placing up to four GPUs on the same node made the runs slightly slower rather than faster. This floor is therefore set by the latency of the many small kernels and global reductions of each AMG V-cycle rather than by the network, and strong scaling pays off only while the work per GPU stays well above it. B2 does not fit on a single GPU: with the current implementation, the memory required by the cuDOLFINx forms limits a GH200 (96 GB) to about five million DoFs. Table 9 reports the weak scaling of the linear solvers. The number of Krylov iterations generally decreases as the mesh is refined. This is consistent with the theory in [5]: uniform refinement shifts the problem, spectrally, closer to a pure Laplacian, for which robust or even faster convergence is expected. The iteration counts therefore do not explain the loss of per-node throughput observed at larger sizes. On the GPU the AMG-CG throughput drops from 4.81 to 1.94 MDoFs/s per node between B1 and B4, while on the CPU it stays roughly constant up to B3 and then falls to 1.75. This degradation is likely due to parallel overheads at higher node counts. The GPU is faster than the CPU on B1–B3, up to about a factor two on B1, while the two are comparable on

16

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

Weak Scaling - Structured Mesh

Runtime [s]

Assembly 0.2

1.4 1.2

0.15

1 0.8

GPU CPU

0.1

CG solver

0.6

GPU Jacobi GPU AMG CPU AMG

0.4 5

10

15 #106

DoFs

5

10

DoFs

15 #106

Fig. 5. Average runtimes per timestep for the EMI linear solvers and assembly operations for a weak scaling test on regular meshes (A1-A4). Corresponds to Table 7.

B1 Nodes 1 2 4 8

B2

GPU

CPU

GPU

CPU

0.70 (68) 0.75 (68) 0.67 (66) 0.62 (68)

1.57 (66) 0.99 (66) 0.92 (66) 0.96 (64)

OOM 0.97 (56) 0.67 (54) 0.58 (56)

2.33 (57) 1.33 (56) 0.95 (55) 0.87 (55)

Table 8 Solution strategy (AMG-CG) strong scaling for the realistic B1 and B2 meshes. Performance numbers are reported in seconds per timestep. The number of Krylov iterations is reported in parentheses. OOM denotes an out-of-memory failure.

B4. AMG-CG is generally as fast as or faster than AMG-GMRES, with very similar iteration counts. We remark that, with the configuration described in Section 3.5 and the meshes cleaned of nearly degenerate tetrahedra, a single set of solver options converged for all the entries of Tables 8 and 9, with no per-entry tuning. Weak Scaling - Unstructured Mesh

Assembly

Runtime [s]

0.6

AMG-CG solver GPU CPU

0.5

1.4 1.2

0.4

1 0.8

0.3 0.5

1

DoFs

1.5

2 #107

0.5

1

DoFs

1.5

2 #107

Fig. 6. Average runtimes per timestep for the EMI AMG-CG linear solver and assembly operations for a weak scaling test on realistic geometry meshes. Corresponds to Table 9.

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

Model B1 B2 B3 B4

17

GPU AMG-CG

GPU AMG-GMRES

CPU AMG-CG

CPU AMG-GMRES

4.81 (68) 3.40 (56) 3.80 (48) 1.94 (50)

4.32 (72) 3.23 (56) 2.91 (60) 1.76 (52)

2.17 (66) 2.48 (56) 2.54 (48) 1.75 (58)

1.89 (71) 2.57 (51) 2.01 (58) 1.42 (69)

Table 9 Linear solver weak scaling for a sequence of realistic geometry meshes. Performance numbers are reported in millions of DoFs per second, normalized by the number of nodes. Larger numbers indicate better performance. Krylov iteration counts are reported in parentheses for the linear solvers. The corresponding runtimes in seconds are depicted in Figure 6.

6. Conclusions. This work presents the first performance evaluation and comparison of GPU-accelerated linear algebra routines within cuDOLFINx. We consider the EMI (Extracellular-Membrane-Intracellular) model as a mixed-dimensional benchmark on unstructured meshes. EMI is an elliptic model posed on a domain partitioned into subdomains coupled through interface conditions, and its numerical solution is of significant practical interest [34, 6, 8, 33]. The presence of interfaces makes the EMI problem more challenging to assemble and solve than a standard elliptic problem, making it a meaningful stress test for the cuDOLFINx framework. We presented a GPU-accelerated implementation of the EMI model and evaluated its performance against a CPU-based implementation across both regular and realistic geometries. We assessed right-hand side assembly and linear solver performance using strong and weak scaling experiments on multiple super-computing systems. For the assembly routines, strong scaling results showed that a single GH200 GPU node is roughly equivalent to two 144-core CPU nodes. Weak scaling results revealed a similar but geometry-dependent trend: each GPU node was equivalent to approximately three CPU nodes for the regular meshes, and 30% to 100% more efficient than a single CPU node for the realistic meshes. In all cases, the weak scaling of assembly routines was close-to-optimal. For linear solvers, the results were more nuanced and depended on geometry, preconditioning strategy, and problem scale. For the regular geometry, Jacobi preconditioning outperformed AMG on GPUs at low node counts despite requiring substantially more Krylov iterations, though this trend reversed as the number of nodes increased and communication costs began to dominate. The CPU solver failed to scale past two nodes, a consequence of the much larger number of MPI processes per node compared to the single-process-per-node GPU configuration. In the weak scaling regime, the GPU implementation exhibits near-optimal scaling and eventually overtook the CPU solver on a per-node basis at the largest problem size tested. For realistic meshes, both architectures are robust in terms of iteration counts and show similar strong scaling, with the runtimes close to saturation after a few nodes. A GPU node is consistently faster than a CPU node, by a factor between 1.3 and 2.2. In the weak scaling tests, the CPU solver retains a close-to-optimal parallel efficiency, whereas the GPU solver scales less well but remains overall faster, with one GPU node being between 1.1 and 2.2 times faster than one CPU node. The iteration counts are robust also in this case. We emphasize that the one-to-one node comparisons used throughout this study are adopted for convenience and do not reflect equivalence in cost, power consumption, or hardware architecture between GPU and CPU nodes. Overall, the experiments show, for both assembly and solution, that one GH200 GPU node is equivalent to roughly 160 to 310 CPU cores, depending on the

18

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

setting. Finally, we remark that robust convergence on the realistic geometries requires care. Hybrid smoothers, which apply Gauss–Seidel within each process subdomain and Jacobi across subdomains, depend on the parallel partition and may yield AMG preconditioners that are not positive definite. Removing degenerate mesh elements after careful inspection, and adopting the same partition-independent AMG configuration, based on Chebyshev smoothing, on both architectures, made a single solver setting converge in all the tests, with iteration counts essentially independent of the architecture and of the number of nodes. On the GPU, the strong scalability on these meshes is limited by the latency of the AMG cycle and by the memory footprint of the current implementation (about five million DoFs per GH200), which are natural targets for future work.

Appendix A. Detailed solver parameters. This appendix lists the specific configuration parameters used for the PETSc solvers and hypre BoomerAMG preconditioners. For a detailed explanation of these options, the reader is referred to [23, 28]. We used PETSc 3.23.7, and hypre 2.32.0.

A.1. Regular geometry configurations (Model A). The solver configuration for regular geometries employed the following baseline AMG settings both for CPU and GPU:

# -ksp_type here is cg or gmres -ksp_rtol: 1e-7 -ksp_norm_type: unpreconditioned -pc_hypre_boomeramg_relax_type_all: l1scaled-SOR/Jacobi -pc_hypre_boomeramg_agg_nl: 0 -pc_hypre_boomeramg_strong_threshold: 0.25

Coarsening and interpolation were left to the defaults of each architecture (HMIS on CPUs, PMIS on GPUs, see Section 3.5). For comparison, the Jacobi preconditioner was configured as follows:

-ksp_type: cg -ksp_rtol: 1e-7 -ksp_norm_type: unpreconditioned -pc_type: jacobi

A.2. Realistic geometry configurations (Model B). For the realistic models, CPU and GPU runs used exactly the same options. CG, with homogeneous Neumann conditions and the constant nullspace, and GMRES, with one Dirichlet point, share the following configuration:

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

19

# -ksp_type here is cg or gmres -ksp_rtol: 1e-7 -ksp_norm_type: unpreconditioned -pc_type: hypre -pc_hypre_boomeramg_coarsen_type: PMIS -pc_hypre_boomeramg_interp_type: ext+i -pc_hypre_boomeramg_relax_type_all: Chebyshev -pc_hypre_boomeramg_relax_type_coarse: Chebyshev -pc_hypre_boomeramg_no_CF: true -pc_hypre_boomeramg_agg_nl: 0 -pc_hypre_boomeramg_strong_threshold: 0.25 -pc_hypre_boomeramg_min_coarse_size: 2 PMIS coarsening and extended+i interpolation are the GPU defaults; they are set explicitly so that the CPU runs use them as well. C/F-ordered relaxation, enabled by default on CPUs only, is disabled. Aggressive coarsening is disabled because PETSc fixes different interpolation operators on the aggressively coarsened levels on CPUs (multipass) and GPUs (two-stage extended+e). The minimum coarse-grid size of two prevents a 1 × 1 coarsest matrix, which vanishes for the singular Neumann operator. The hypre statistics (number of levels and operator complexity) confirmed that CPU and GPU runs build the same hierarchy. Acknowledgments. BP gratefully acknowledges the use of the Vista system at the Texas Advanced Computing Center under the “ADCIRC” allocation. EC acknowledges the CINECA award under the ISCRA initiative, for the availability of high-performance computing resources and support (project DDO2CARD). EC, JS and JT have been supported by the European High-Performance Computing Joint Undertaking (JU) EuroHPC under grants agreement No 955495 (MICROCARD) and No 101172576 (MICROCARD-2). The JU receives support from the European Union’s Horizon 2020 research and innovation program and France, Italy, Germany, Austria, Norway, Switzerland. PB and EC acknowledge the support by grants from the Swiss National Supercomputing Centre (CSCS) under project IDs sm111 and lp157 on Alps. The authors would like to thank Joshua Steyer’s PI, Axel Loewe, for helpful discussions and for the freedom to pursue this independent study. The authors also wish to thank Mark Potse for his help with the mesh utilities. AI tools were used for language and grammar editing; the authors assume responsibility for all content. REFERENCES [1] R. R. Aliev and A. V. Panfilov, A simple two-variable model of cardiac excitation, Chaos, Solitons & Fractals, 7 (1996), pp. 293–301. [2] F. Ballarin, multiphenicsx - easy prototyping of multiphysics problems in FEniCSx, 2020, https://multiphenics.github.io/. [3] I. A. Baratta, J. P. Dean, J. S. Dokken, M. Habera, J. S. Hale, C. N. Richardson, M. E. Rognes, M. W. Scroggs, N. Sime, and G. N. Wells, DOLFINx: The next generation FEniCS problem solving environment , Dec. 2023, https://doi.org/10.5281/ zenodo.10447666. [4] P. Benedusi, A. J. Ellingsrud, H. Herlyng, and M. E. Rognes, Scalable approximation and solvers for ionic electrodiffusion in cellular geometries, SIAM Journal on Scientific Computing, 46 (2024), pp. B725–B751, https://doi.org/10.1137/24M1644717. [5] P. Benedusi, P. Ferrari, M. Causemann, and S. Serra-Capizzano, Dense cell-by-cell systems of PDEs: approximation, spectral analysis, and preconditioning, BIT Numerical

20

E. CENTOFANTI, B. PACHEV, J. STEYER, J.D. TROTTER, P. BENEDUSI

Mathematics, 66 (2026), https://doi.org/10.1007/s10543-026-01136-6. [6] P. Benedusi, P. Ferrari, M. E. Rognes, and S. Serra-Capizzano, Modeling excitable cells with the EMI equations: spectral analysis and iterative solution strategy, Journal of Scientific Computing, 98 (2024), p. 58, https://doi.org/10.1007/s10915-023-02449-2. [7] M. Causemann, V. Vinje, and M. E. Rognes, The pulsating brain: An interface-coupled fluid-poroelastic interaction model of the cranial, Mathematical Modelling of the Human Brain II: From Glymphatics to Deep Learning, (2026), p. 47. [8] E. Centofanti, N. M. M. Huynh, L. F. Pavarino, and S. Scacchi, Parallel algebraic multigrid solvers for composite discontinuous galerkin discretization of the cardiac EMI model in heterogeneous media, Computer Methods in Applied Mechanics and Engineering, 442 (2025), p. 118001. [9] E. Centofanti, N. M. M. Huynh, L. F. Pavarino, and S. Scacchi, GDSW Preconditioners for the EMI Model in Cardiac Electrophysiology, in Proceedings of Domain Decomposition Methods in Science and Engineering XXIX, Springer, 2026, pp. 249–256. [10] E. Centofanti and S. Scacchi, A comparison of Algebraic Multigrid Bidomain solvers on hybrid CPU-GPU architectures, Computer Methods in Applied Mechanics and Engineering, 423 (2024), p. 116875. [11] C. Dapogny, C. Dobrzynski, and P. Frey, Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems, Journal of Computational Physics, 262 (2014), pp. 358–378, https://doi.org/https://doi.org/10. 1016/j.jcp.2014.01.005. [12] H. De Sterck, U. M. Yang, and J. J. Heys, Reducing complexity in parallel algebraic multigrid preconditioners, SIAM Journal on Matrix Analysis and Applications, 27 (2006), pp. 1019–1039, https://doi.org/10.1137/040615729. [13] A. J. Ellingsrud, P. Benedusi, and M. Kuchta, A splitting, discontinuous Galerkin solver for the cell-by-cell electroneutral Nernst–Planck framework, SIAM Journal on Scientific Computing, 47 (2025), pp. B477–B504, https://doi.org/10.1137/24M1653367. [14] R. D. Falgout, R. Li, B. Sjögreen, L. Wang, and U. M. Yang, Porting hypre to heterogeneous computer architectures: Strategies and experiences, Parallel Computing, 108 (2021), p. 102840, https://doi.org/10.1016/j.parco.2021.102840. [15] M. Favino, Fractures and thin heterogeneities as Robin-Wentzell interface conditions, Applied Mathematical Modelling, (2025), p. 116355. [16] E. Fevola, F. Ballarin, L. Jimenez-Juan, P. Triverio, S. Grivet-Talocia, and G. Rozza, Automated estimation of patient-specific boundary conditions for cardiovascular simulations: an optimal control approach, in 2021 Summer Biomechanics, Bioengineering, & Biotransport Conference, SB3C Foundation, Inc., 2021. [17] P. C. Franzone, L. F. Pavarino, and S. Scacchi, Mathematical cardiac electrophysiology, Springer, 2014. [18] I. Fumagalli, Discontinuous Galerkin method for a three-dimensional coupled fluid-poroelastic model with applications to brain fluid mechanics, Mathematics in Engineering, 7 (2025), pp. 130–161. [19] F. Göbel, N. M. M. Huynh, F. Chegini, L. F. Pavarino, M. Weiser, S. Scacchi, and H. Anzt, A BDDC preconditioner for the cardiac EMI model in three dimensions, SIAM Journal on Scientific Computing, 48 (2026), pp. A646–A667. [20] Z. Guo, C. Cui, and Q. Zhang, GPU-based matrix-free multigrid methods of stable generalized finite element methods for interface problems, Computer Methods in Applied Mechanics and Engineering, 461 (2026), p. 119275, https://doi.org/https://doi.org/10.1016/j.cma. 2026.119275. [21] A. Heinlein, A. Klawonn, S. Rajamanickam, and O. Rheinbach, FROSch: A fast and robust overlapping Schwarz domain decomposition preconditioner based on Xpetra in Trilinos, in International Conference on Domain Decomposition Methods, Springer, 2018, pp. 176–184. [22] N. M. M. Huynh, F. Chegini, L. F. Pavarino, M. Weiser, and S. Scacchi, Convergence analysis of BDDC preconditioners for composite DG discretizations of the cardiac cell-bycell model, SIAM Journal on Scientific Computing, 45 (2023), pp. A2836–A2857. [23] Hypre online documentation. https://hypre.readthedocs.io/en/latest/solvers-boomeramg. html. Accessed: 2026-03-27. [24] K. H. Jæger, K. G. Hustad, X. Cai, and A. Tveito, Efficient numerical solution of the EMI model representing the extracellular space (E), cell membrane (M) and intracellular space (I) of a collection of cardiac cells, Frontiers in Physics, 8 (2021), p. 579461. [25] R. Li, B. Sjögreen, and U. M. Yang, A new class of AMG interpolation methods based on matrix-matrix multiplications, SIAM Journal on Scientific Computing, 43 (2021), pp. S540–

MULTI-GPU ACCELERATION OF INTERFACE PROBLEMS IN FENICS

21

S564, https://doi.org/10.1137/20M134931X. [26] B. Pachev, J. Trotter, and I. Baratta, Multi-GPU acceleration of FEniCSx codes with cuDOLFINxx, in FEniCS 2025, 2025. [27] B. A. Pachev, J. D. Trotter, and I. A. Baratta, cuDOLFINx: A CUDA extension for FEniCSx, in The FEniCS Project: The FEniCS 2024 Conference, Springer, 2026, pp. 27– 38. [28] PETSc online documentation. https://petsc.org/release/manualpages/. Accessed: 2026-03-27. [29] R. Plonsey, R. C. Barr, and A. Bioelectricity, Quantitative approach, Cham, Switzerland: Springer, (2007). [30] M. Potse, L. Cirrottola, and A. Froehly, A practical algorithm to build geometric models of cardiac muscle tissue, in 8th European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2022), 01 2022, https://doi.org/10.23967/eccomas. 2022.027. [31] A. Shapson-Coe, M. Januszewski, D. R. Berger, A. Pope, Y. Wu, T. Blakely, R. L. Schalek, P. H. Li, S. Wang, J. Maitin-Shepard, et al., A petavoxel fragment of human cerebral cortex reconstructed at nanoscale resolution, Science, 384 (2024), p. eadk4858. [32] J. Steyer, F. Chegini, M. Potse, A. Loewe, and M. Weiser, Continuity of microscopic cardiac conduction in a computational cell-by-cell model, in 2023 Computing in Cardiology Conference (CinC), vol. 50, 2023, https://doi.org/10.22489/CinC.2023.385. [33] J. Steyer, F. Chegini, T. Stary, M. Potse, M. Weiser, and A. Loewe, Electrograms in a cardiac cell-by-cell model, in Workshop Biosignals, Göttingen, Germany, Feb. 2024, https://doi.org/10.47952/gro-publ-194, https://doi.org/10.47952/gro-publ-194. [34] A. Tveito, K.-A. Mardal, and M. E. Rognes, Modeling excitable tissue: The EMI framework, Springer, 2021.

Related documents

• Warp-Geo: Differentiable Geometry Representation for Dynamic-Boundary Simulation and Shape Optimization#1034936
• GLASS: Architecture-Tuned, Composable, Device-Side Linear Algebra for Edge Robotics and Beyond#1064723
• Kamino: GPU-based Simulation of Rigid Multi-Body Physics with Kinematic Loops#1411766
• Automating Cardiac Segmentation for Complex Congenital Heart Disease 3D Modeling Using nnU-Net#651657
• Agentic Kernel Optimization: Generating State-of-the-Art GPU Kernels Without Hand-Written CUDA#610063