<figure><img src="image_1.png" /></figure>

HAL
open science

# Task-Based Runtime Support and Programming for Finite Element Simulations

Paul Bouchaud

► To cite this version:

Paul Bouchaud. Task-Based Runtime Support and Programming for Finite Element Simulations. ENS de Lyon; Inria. 2025. hal-05224253

HAL Id: hal-05224253

https://inria.hal.science/hal-05224253v1

Submitted on 26 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

ÉCOLE NORMALE
SUPÉRIEURE
DE LYON

<figure><img src="image_4.png" /></figure>

ENS DE LYON
UNIVERSITÉ CLAUDE BERNARD LYON 1

DEPARTMENT OF COMPUTER SCIENCE
INTERNSHIP AT INRIA BORDEAUX

10/06/2025-1/08/2025

# Task-Based Runtime Support and Programming for Finite Element Simulations

Student :
Paul BOUCHAUD

Instructor
Olivier AUMAGE

<figure><img src="image_5.png" /></figure>

August 26, 2025
Abstract

My work studies task-based runtime support for finite element method (FEM)
assembly using StarPU. The idea is to replace bulk-synchronous loops with paral-
elized approach. Local kernels were rewritten in modern C++20/23, relying on
non-owning views (std::span, std::mdspan) and a preallocated argument table to
avoid repeated allocations. I also used METIS for partitioning and a greedy coloring
step to guarantee conflict-free execution.

The implementation is done for the project DOLFINx/FEniCS. On the 2D Poisson problem, assembly is already cheap and preprocessing dominates, so there is little to no gain and even loss of time on small model. On a 3D hyperelasticity benchmark, where each element is more expensive, the task-based model works well. Grouping cells into partitions reduces the number of tasks and keeps scheduling overhead low while maintaining balance. These results suggest that task-based assembly is a practical choice for nonlinear or compute-intensive models. Future extensions will need to focus on NUMA-aware scheduling, custom StarPU schedulers, and GPU support.
# Contents

<table><tr><td>1 Introduction</td><td>4</td></tr><tr><td>2 Context</td><td>5</td></tr><tr><td>2.1 Mathematical Modeling and Variational Formulation</td><td>5</td></tr><tr><td>2.2 Discretization Techniques and the Finite Element Philosophy</td><td>5</td></tr><tr><td>2.3 Challenges in Implementing Finite Element Codes</td><td>6</td></tr><tr><td>2.4 The FEniCS Project and Its Architecture</td><td>6</td></tr><tr><td>2.5 Parallel Execution with MPI and PETSc</td><td>7</td></tr><tr><td>2.6 Task-Based Runtime Systems</td><td>7</td></tr><tr><td>2.7 StarPU: Runtime Model and Execution Engine</td><td>8</td></tr><tr><td>2.8 Reformulating FEM Assembly as a Task Graph</td><td>8</td></tr><tr><td>2.9 Modern C++ for Kernel and Interface Design</td><td>8</td></tr><tr><td>2.10 Data Locality and Memory Layout Optimization</td><td>9</td></tr><tr><td>3 Proposition</td><td>9</td></tr><tr><td>3.1 Parallel Conflict Management: Coloring vs Partitioning</td><td>9</td></tr><tr><td>3.2 Local Assembly Execution Model: Thread-Private vs Task-Based</td><td>10</td></tr><tr><td>3.3 Scheduling Flexibility and Adaptability</td><td>10</td></tr><tr><td>3.4 Extensibility to GPUs and Heterogeneous Systems</td><td>10</td></tr><tr><td>3.5 Implementation Complexity and Maintainability</td><td>11</td></tr><tr><td>4 Implementation Strategy</td><td>11</td></tr><tr><td>4.1 Pre-allocated Argument Table for Task Management</td><td>11</td></tr><tr><td>4.2 CPU Core Management and Main Thread Reservation</td><td>12</td></tr><tr><td>4.3 Partitioning Strategy: METIS</td><td>12</td></tr><tr><td>4.4 Coloring Strategy for Conflict-Free Parallelism</td><td>13</td></tr><tr><td>4.5 Kernel Code Reuse and StarPU Codelet Interface</td><td>13</td></tr><tr><td>4.6 Debugging and Profiling Strategy</td><td>14</td></tr><tr><td>5 Evaluation</td><td>14</td></tr><tr><td>5.1 Hardware and Software Environment</td><td>14</td></tr><tr><td>5.2 Benchmark 1: Scalar Poisson Problem</td><td>15</td></tr><tr><td>5.3 Hyperelasticity Benchmark (3D)</td><td>17</td></tr><tr><td>6 Conclusion</td><td>19</td></tr><tr><td>A Annex – Internship Context</td><td>20</td></tr></table>

# List of Figures

<table><tr><td>1</td><td>Partitioning of the 640 × 320 mesh into 108 parts, colored by partition index.</td><td>15</td></tr><tr><td>2</td><td>Partitioning of the 6400 × 3200 mesh into 72 parts (colored by partition index).</td><td>16</td></tr><tr><td>3</td><td>Greedy coloring of the 108 partitions for the 640 × 320 mesh.</td><td>16</td></tr><tr><td>4</td><td>Preliminary StarPU trace for Poisson assembly on the 6400 × 3200 mesh. Green bars indicate full utilization of all worker threads.</td><td>17</td></tr></table>
<table><tr><td>5</td><td>3D partitioning of the tetrahedral mesh (color per partition).</td><td>17</td></tr><tr><td>6</td><td>StarPU trace of task execution for 3D hyperelasticity benchmark (108 partitions on 9 threads).</td><td>18</td></tr></table>
# 1 Introduction

Most partial differential equations that arise in engineering cannot be solved analytically. Once the geometry becomes complex, the material behavior nonlinear, or the boundary conditions irregular, exact solutions are out of reach. The finite element method (FEM) has become the standard tool in such cases, thanks to its geometric flexibility and solid variational foundation, which ensures mathematical rigor and convergence guarantees [14, 11, 6, 4]. The price to pay is computational cost. At the core of every FEM simulation lies the assembly of the global system matrix: integrating local element contributions, mapping them through the degree-of-freedom numbering, and inserting them into large sparse matrices. For realistic problem sizes this step dominates runtime, and in both serial and parallel implementations it is the main bottleneck. The challenge is even sharper on modern HPC machines. Current frameworks such as FEniCS [13, 1] rely on bulk-synchronous MPI parallelism, sometimes combined with OpenMP. While efficient for regular workloads, this model is ill-suited to the irregular and fine-grained nature of element assembly.

The FEniCS project provides a full software ecosystem [13] for automating PDE solutions with FEM. The Unified Form Language (UFL) lets the user describe variational problems at a high level, while the DOLFINx backend manages the mesh, linear algebra, and parallel execution. DOLFINx also integrates runtime code generation via FFCx and Basix, so that the form compiler can produce optimized kernels specialized for the chosen element types and problem definition.

Yet despite these strengths, FEniCS inherits the limitations of bulk-synchronous parallelism. DOLFINx relies on MPI for distributed memory and optionally OpenMP for shared-memory acceleration. These tools are effective for structured workloads on homogeneous platforms, but they fall short for fine-grained, irregular, or dynamic computations. Static domain partitioning means that load balancing and memory locality are essentially fixed at preprocessing time, leaving little flexibility to adapt during execution.

Task-based programming offers an alternative. Instead of forcing computations into rigid parallel phases, the work is expressed as fine-grained tasks with explicit data dependencies. The runtime then schedules them dynamically according to resource availability. This paradigm matches naturally with FEM assembly: element contributions are mostly independent, yet must eventually be accumulated into shared global structures.

StarPU [2, 8] represents one of the most mature implementations of this paradigm. Developed at Inria, it manages task graphs across heterogeneous architectures, handling the complex orchestration of CPU and GPU resources automatically. The runtime profiles task performance and adapts scheduling policies accordingly, from simple work-stealing to more sophisticated performance-aware strategies.

However, retrofitting an existing FEM codebase to use task-based execution is non-trivial. The assembly kernels that currently operate within tight loops must be extracted as independent task functions. More challenging still, the implicit data dependencies embedded in global matrix assembly must be made explicit to the task scheduler. Each task needs carefully designed interfaces to avoid race conditions while minimizing memory overhead. The traditional FEM concepts of mesh partitioning and ghost layer communication don't map cleanly to a task-centric worldview.

This internship focuses on creating a modular approach to connect local FEM assembly, written in modern C++23, with the StarPU runtime. The aim is not to rebuild the
entire FEniCS framework, but rather to extract the local assembly loop and represent it as a set of tasks managed by StarPU. A key question is how to keep the flexibility and clarity of FEM while also achieving high efficiency and portability across different hardware.

To achieve this, we use modern C++ features such as `std::span` and `std::mdspan` (introduced in C++20 and C++23) to represent memory views in a safe and efficient way. These features keep the kernel code readable, avoid unnecessary dynamic memory allocation, and help the compiler optimize the code. They also allow control over data alignment and memory layout, which is important for vectorization and better use of the cache.

In the end, this work aims to show that combining modern C++ tools with a task-based runtime can produce scalable, portable, and maintainable FEM assembly. By carefully designing task interfaces, scheduling strategies, and memory handling, integrating StarPU into a FEM codebase can provide both high performance and the ability to adapt to different computing architectures.

## 2 Context

Numerical simulation of physical systems governed by partial differential equations (PDEs) has become a foundational tool in both scientific research and industrial engineering. PDEs arise naturally in the formulation of conservation laws and constitutive equations that model diverse phenomena such as thermal diffusion, structural deformation, wave propagation, electromagnetism, or fluid dynamics. They provide a compact and expressive representation of physical laws by linking the variation of fields (temperature, displacement, pressure, etc.) across space and time to governing dynamics.

However, PDEs are notoriously difficult to solve analytically except in highly idealized scenarios. Realistic problems typically involve complex geometries, heterogeneous material properties, nonlinear interactions, and sophisticated boundary conditions. These features require robust and flexible numerical methods to approximate solutions with sufficient accuracy and efficiency.

### 2.1 Mathematical Modeling and Variational Formulation

The numerical solution of a PDE begins by rewriting it in weak (variational) form. This step is essential: it allows us to treat functions that are not classically differentiable, and it provides a rigorous way to incorporate boundary conditions. The weak formulation is obtained by multiplying the PDE with test functions, integrating by parts, and requiring the residual to vanish in a suitable function space.

### 2.2 Discretization Techniques and the Finite Element Philosophy

Once the weak form is established, it must be discretized to obtain a computable linear or nonlinear system. Several methods exist for this step—finite differences, finite volumes, and finite elements. Among them, the finite element method (FEM) is the most flexible when dealing with irregular geometries and heterogeneous coefficients. Unlike finite differences, which depend on structured Cartesian grids, FEM adapts naturally to unstructured meshes and curved boundaries.
The domain $\Omega$ is subdivided into elements (triangles, quadrilaterals in 2D; tetrahedra, hexahedra in 3D). On each element, a polynomial basis is defined, and the global approximation is formed by combining these local functions through shape function mappings. Transformations from a reference element to physical elements introduce Jacobians and other geometric factors, which FEM libraries handle automatically.

Discretization produces a sparse linear system

$$Ax = b,$$

where $A$ is typically a stiffness or mass matrix and $b$ the right-hand side from forcing terms and boundary conditions. In nonlinear settings, both $A$ and $b$ depend on the unknown solution and require iterative schemes such as Newton or Picard methods.

## 2.3 Challenges in Implementing Finite Element Codes

Despite its conceptual elegance, FEM is difficult to implement efficiently. Assembly requires evaluating local integrals on each element, mapping them through global degrees of freedom, and inserting them into large sparse matrices. This involves indirect addressing, nontrivial memory layouts, and careful synchronization in parallel contexts.

Sparsity is both an advantage and a challenge: most entries of the global matrix vanish, but sparse insertion is far less efficient than dense arithmetic. Specialized formats (CSR, block structures) and optimized linear algebra libraries (PETSc, Trilinos, Eigen) are therefore essential.

The complexity increases for nonlinear and time-dependent problems, where matrices and vectors must be reassembled at every iteration or timestep, and for multiphysics simulations, where coupled fields introduce extra blocks and constraints. These difficulties explain the need for automation and high-level frameworks, which abstract away low-level details while still enabling efficient implementations.

## 2.4 The FEniCS Project and Its Architecture

FEniCS is a widely used open-source platform for automating FEM simulations. It covers the full PDE workflow, from high-level mathematical specification to low-level numerical execution. At its core is UFL (Unified Form Language), a domain-specific language embedded in Python, which lets users express variational forms using near-mathematical syntax. UFL supports symbolic differentiation, integration rules, and form validation.

Once a form is defined, it is compiled by FFCx (FEniCS Form Compiler eXperimental), which generates low-level C++ code specialized for the chosen finite element, degree, and variational form. This code computes local element matrices and vectors efficiently. Basix provides the definitions of basis functions and supports a wide range of elements, including Lagrange, Raviart-Thomas, Nédélec, and discontinuous Galerkin elements.

DOLFINx acts as the runtime environment. It handles mesh representation, degree-of-freedom (DoF) mappings, assembly routines, and interfaces to PETSc for linear and nonlinear solvers. DOLFINx also manages parallelism through distributed meshes and parallel data structures. Its modular design exposes both Python and C++ APIs, and just-in-time compilation enables runtime specialization of problem-dependent kernels.
## 2.5 Parallel Execution with MPI and PETSc

Parallelism in FEniCS relies primarily on MPI-based domain decomposition. The mesh is partitioned across MPI ranks using graph partitioners such as METIS or SCOTCH. Each process owns a subdomain of the mesh and the associated DoFs. Inter-process communication is managed through ghost cells, which replicate neighboring element data for local use.

Assembly proceeds in parallel: each rank computes local element contributions, inserts them into its portion of the global matrix, and synchronizes ghost DoFs at subdomain boundaries. PETSc provides scalable data structures (Mat, Vec) and solvers (KSP, SNES) [3] for the global linear or nonlinear system. This model scales well for structured problems and homogeneous compute environments.

However, the approach is static. Once the mesh is partitioned, the execution flow is fixed. Processes advance in lockstep, with collective communication introducing synchronization barriers. This rigidity makes it difficult to adapt to load imbalance and limits efficiency on heterogeneous hardware.

## Limitations of Bulk-Synchronous Models

Bulk-synchronous parallelism (BSP) organizes computation into repeating phases: compute, communicate, synchronize. While this model is simple to reason about, it introduces several limitations in complex simulations.

First, BSP assumes balanced workloads. In adaptive mesh refinement (AMR), multi-physics coupling, or heterogeneous material models, some subdomains require much more work than others. Static partitioning cannot adapt, leading to idle resources and reduced efficiency.

Second, BSP does not overlap communication with computation. Many modern applications could benefit from pipelines where partial results are used before global completion, but BSP enforces global synchronization at every phase.

Third, BSP assumes a homogeneous, synchronous hardware environment. This assumption breaks down on current platforms with GPUs, CPUs running at different frequencies, NUMA effects, and network congestion.

## 2.6 Task-Based Runtime Systems

Task-based programming models address these shortcomings by expressing computation as fine-grained tasks with explicit data dependencies. A runtime system builds a directed acyclic graph (DAG) of dependencies and schedules tasks dynamically according to readiness, priority, and available resources.

This separation between specification and execution enables optimizations such as dynamic load balancing, communication hiding, heterogeneity-aware scheduling, and adaptive control flow. The task graph can even evolve at runtime, depending on the data.

Notable task-based runtimes include StarPU (Inria) [2], OmpSs (BSC) [9], PaRSEC (ICL) [5], and Kokkos (Sandia) [10]. All provide different APIs and capabilities, but share the same philosophy: abstract parallelism through data-flow graphs and let the runtime manage execution.
## 2.7 StarPU: Runtime Model and Execution Engine

StarPU is a mature task-based runtime system for heterogeneous multicore and accelerator-based platforms. It provides a C/C++ API to define tasks, manage data, and declare dependencies. Each task is associated with a *codelet*, a function pointer implementing the task's kernel, and a set of data handles describing the memory regions it accesses.

Data must be registered with StarPU, which then ensures consistency across devices. The runtime manages memory transfers and schedules tasks on available resources according to chosen policies. StarPU offers multiple schedulers, including eager (first-ready, first-served), DMDA (data-aware), HEFT (heterogeneous earliest finish time), and user-defined models.

StarPU also supports performance models that estimate task execution costs on different devices. These can be learned online or provided manually, allowing the runtime to make intelligent scheduling choices in hybrid CPU/GPU systems.

By managing data movement, dependencies, and scheduling dynamically, StarPU enables overlapping of communication and computation, adapts to workload imbalance, and makes efficient use of accelerators. These properties are particularly attractive for irregular workloads such as FEM assembly.

## 2.8 Reformulating FEM Assembly as a Task Graph

Applying StarPU to FEM assembly requires rethinking the computation. Instead of iterating over elements sequentially, each element (or group of elements) becomes a task. A task receives the element geometry, coefficients, and DoF mappings, and produces a local contribution to the global system.

Tasks that share global DoFs must be serialized or partitioned to avoid conflicts. A common strategy is mesh coloring: assigning colors so that no two elements of the same color share DoFs, and executing tasks of one color in parallel. Another option is graph partitioning with METIS, producing disjoint subdomains that can be assembled concurrently.

To reduce memory traffic, StarPU tasks should operate on *views* rather than copies of data. Modern C++ abstractions such as `std::span` and `std::mdspan` are useful here, providing lightweight, safe references to contiguous or multi-dimensional arrays without allocations.

## 2.9 Modern C++ for Kernel and Interface Design

The introduction of `std::span` in C++20 and `std::mdspan` in C++23 provides new tools for writing high-performance numerical kernels. These non-owning views over memory regions support optional bounds checking and customizable layouts. They avoid raw pointers and simplify indexing logic.

For FEM kernels, `mdspan` enables expressive multi-dimensional indexing for coordinates, basis functions, and local element matrices. Combined with templates and `constexpr`, this yields generic kernels that adapt to element type, polynomial order, or vector dimension.

These abstractions also cleanly separate numerical logic from task interfaces. StarPU tasks can directly receive spans or mdspans, avoiding unnecessary copies. By matching
memory layout to hardware expectations (e.g., row-major ordering for cache reuse), one achieves both readability and performance.

## 2.10 Data Locality and Memory Layout Optimization

In FEM assembly, memory access patterns often dominate performance more than floating-point throughput. Traditional implementations suffer from irregular access to geometry, DoF maps, and coefficients, leading to poor cache reuse and high latency on NUMA systems.

The task-based model makes it possible to optimize data layout explicitly. Using modern C++ views such as `std::span` and `std::mdspan`, tasks access contiguous memory in row-major order, reducing cache misses and improving spatial locality.

Moreover, since the runtime tracks dependencies, StarPU can schedule dependent tasks on the same core or NUMA node, reducing cross-socket traffic and memory contention. This is particularly beneficial for bandwidth-bound operations such as Jacobian evaluation or basis function application.

While each optimization may be incremental, the gains accumulate across thousands of elements in large-scale simulations. Aligning task-level memory access with hardware architecture yields measurable performance improvements without resorting to architecture-specific code.

# 3 Proposition

## 3.1 Parallel Conflict Management: Coloring vs Partitioning

The biggest challenge to parallelize FEM assembly comes from the shared degree of freedom of multiple elements, it could create conflict and add randomness and error in an already flawed modeling. One of the basic approaches to this issue is graph coloring, elements sharing DOFs are given the same color so that they can be executed simultaneously. While it is a safe approach to the problem and is error-proof, it has poor cache-locality and as such cannot be used for larger size. More sophisticated approaches like space-filling curve (SFC) based greedy coloring can improve spatial coherence, yet the fundamental limitation remains—coloring creates a static, inflexible task structure that doesn't adapt well to load imbalances or heterogeneous hardware.

METIS offers a mesh partitioning that is a different trade-off. Rather than having no conflicts like coloring, we partition the mesh into subdomains with minimal dependencies. This generates fewer, larger groups with better density. While some synchronization points remain where partition boundaries meet, the reduced dependency graph integrates really well with the coloring methods, and enables a better way to partition the graph, error-proof and with better cache-locality. The runtime can adapt to varying task completion times and redistribute work across available resources, something the previous coloring methods could not do. For our implementation this flexibility proved more valuable than the theoretical simplicity of cell-coloring.
## 3.2 Local Assembly Execution Model: Thread-Private vs Task-Based

Traditional OpenMP-based FEM codes often use thread-private assembly, where each thread maintains its own local copy of the global matrix and vector. This approach eliminates write conflicts entirely—threads can assemble elements independently before combining their contributions in a final reduction step. While conceptually straightforward and easy to implement with OpenMP's reduction clauses, the memory overhead quickly becomes an issue for large problems. More importantly, the reduction phase itself becomes a bottleneck, requiring expensive merge operations across dozens or more of private matrices.

The task-based alternative treats each element or partition as an encapsulated computation kernel. Rather than managing thread-private storage, we let StarPU handle the scheduling complexity. Tasks declare their data dependencies explicitly, and the runtime ensures proper execution ordering without manual synchronization. This eliminates both the memory overhead of thread-private copies and the expensive merge phase. No conflicts occur thanks to the coloring and partitioning done in the previous part.

Perhaps more importantly, the task-based model scales beyond shared-memory parallelism. StarPU can distribute tasks across GPUs, remote nodes, or heterogeneous accelerators, something thread-private assembly simply cannot do. The encapsulation also makes performance tuning more tractable, as individual kernels can be profiled and optimized independently of the global execution strategy.

## 3.3 Scheduling Flexibility and Adaptability

Static parallelization schemes stop working due to multitude of reasons, most of which we are not able to take into account in theory. Mesh partitions that look balanced in theory can have vastly different run time due to hardware changes, thread starvation, traditional OpenMP loops offer no help in those cases when some threads finish early while others lag behind and slow the whole process and the slowest of them determines the overall performance.

That is why we decided to use StarPU, whose dynamic scheduling addresses these issues naturally. The runtime continuously monitors task completions and distributes work accordingly. It can also account for NUMA topology, an important capacity as this adaptability is essential in modern HPC systems where core performance can vary due to external factors.

The flexibility of StarPU is essential for heterogeneous execution and for future improvement to GPU kernels. By choosing a task-based approach now, we avoid the need to completely restructure the assembly pipeline when adding accelerator support later. The same task abstraction that enables CPU load balancing extends naturally to GPU offloading and multi-node execution.

## 3.4 Extensibility to GPUs and Heterogeneous Systems

Traditional CPU parallelization strategies hit a wall when GPUs enter the picture. Thread-private assembly assumes abundant memory and cheap synchronization—luxuries that don't exist on GPU architectures. Coloring schemes become even more restrictive when dealing with thousands of lightweight GPU threads that need completely different memory
access patterns. StarPU's abstraction layer sidesteps these issues by treating hardware
as an implementation detail. The same task that runs on a CPU core can be recompiled
for CUDA execution with appropriate kernel modifications. StarPU handles the messy
details of data movement, memory management, and scheduling across heterogeneous
resources. While this project focuses on CPU performance, the underlying architecture
naturally extends to accelerators without major structural changes. This forward com-
patibility proved decisive in our design choice. We've seen too many HPC codes that
were optimized for one generation of hardware only to become obsolete when the next
architecture arrived. Building on StarPU means that future GPU support requires writ-
ing new kernels, not redesigning the entire execution model. Given the rapid evolution of
accelerator technologies, this flexibility seems more valuable than squeezing out the last
few percent of CPU performance with hardware-specific optimizations.

3.5 Implementation Complexity and Maintainability

Traditional approaches accumulate complexity in unexpected places. Implementing robust coloring requires intricate graph algorithms, conflict detection, and validation logic that's difficult to debug when things go wrong. Thread-private assembly seems straightforward until you realize the reduction phase needs to be careful to handle varying matrix sparsity patterns and memory layouts efficiently. StarPU's task model cuts through much of this complexity. Each kernel becomes a self-contained function with clearly defined inputs and outputs—no shared state, no hidden dependencies, no intricate synchronization logic. Tasks can be unit-tested in isolation, profiled independently, and debugged without untangling complex parallel execution flows. When assembly bugs occur, they're localized to specific kernels rather than distributed across threading and reduction code. The modular design also pays dividends for future development. Adding new element types, boundary conditions, or material models requires writing new task kernels without modifying the scheduling infrastructure. Integrating with iterative solvers or preconditioners becomes more natural since the task abstraction extends cleanly to other parts of the FEM pipeline. Perhaps most importantly, this approach eliminates an entire category of concurrency bugs that plague traditional parallel FEM codes. Race conditions, deadlocks, and memory corruption issues that arise from manual thread management simply don't exist when the runtime handles all synchronization. The result is code that's not only faster to develop but also more reliable in production.

4 Implementation Strategy

4.1 Pre-allocated Argument Table for Task Management

Dynamic allocation of task arguments at runtime is costly, especially when thousands
of tasks are created per assembly. Each call to new or malloc adds overhead and leads
to fragmentation, which hurts cache locality and makes performance unpredictable. In
practice, this overhead becomes visible in iterative solvers where the same task structures
are repeatedly created and destroyed.

To avoid this, we rely on a pre-allocated argument table. At initialization, we reserve
space for all task arguments using a container such as `std::vector CELLTaskArgs<T,
 BS> * tab_args`. Once allocated, tasks reuse these argument slots instead of allocating
memory on the fly. This simple change removes the per-task allocation cost and ensures
that arguments live in a contiguous memory region, improving cache reuse.

The layout is also more predictable: stable memory addresses make debugging easier, and cache/TLB behavior is easier to reason about. Shared data such as geometry, material parameters, or degree-of-freedom mappings are passed as pointers inside the argument structures, so multiple tasks can reuse them without duplication. This reduces memory consumption and increases the chance that commonly accessed data stays hot in cache.

## 4.2 CPU Core Management and Main Thread Reservation

In a task-based runtime system like StarPU, careful management of CPU resources is essential to avoid scheduling bottlenecks and ensure smooth execution. StarPU executes computational tasks on a pool of worker threads, but task submission — including task graph generation and insertion — must be performed from a non-worker thread, typically the main thread. If this submission thread is also scheduled as a worker, the runtime risks entering a deadlock or suffering from scheduler starvation, especially when the task queue becomes empty and no thread is available to submit new work.

To prevent such situations, it is necessary to reserve at least one core exclusively for the main thread. This ensures that task submission proceeds independently of the execution of previously submitted tasks. In practice, this means setting the number of StarPU worker threads to `num_available Cores - 1`, where `num_available Cores` is the total number of logical or physical cores available on the machine. This reservation guarantees that at any point during execution, one core remains free to handle submission and scheduling overhead, avoiding deadlocks and maintaining forward progress.

Beyond deadlock avoidance, proper core allocation contributes to system performance. If task submission becomes a bottleneck — for example, if the main thread is unable to generate and submit tasks quickly enough — worker threads may become idle, leading to underutilization of the available hardware. This is particularly relevant when using fine-grained tasks or color-based scheduling, where tasks must be submitted in rounds. Ensuring that task insertion remains faster than task execution is crucial to keeping the pipeline saturated and achieving good scalability.

Finally, the main thread can also take on lightweight coordination tasks such as color-ning, partition tracking, and runtime synchronization between color phases. Reserving a core for it allows these operations to proceed in parallel with computation, contributing to smoother execution and better responsiveness of the parallel system.

## 4.3 Partitioning Strategy: METIS

To parallelize the finite element assembly process while maintaining high performance and correctness, we employ a graph partitioning strategy using METIS. The mesh is first represented as a cell adjacency graph, where each node corresponds to a cell and edges indicate shared degrees of freedom (DoFs) between neighboring cells. METIS partitions this graph into a specified number of subdomains [12] such that the number of inter-partition edges, i.e., shared DoFs, is minimized.

The goal of this partitioning is twofold. First, it reduces potential write conflicts during assembly by grouping together cells that are likely to write to the same entries in the global vector. Second, it improves data locality: since neighboring cells are often
adjacent in memory and share geometric and topological data, grouping them in the same task leads to fewer cache misses and better use of memory bandwidth.

Each resulting partition is then assigned to a StarPU task, which performs assembly on all cells within that subdomain. This coarse-grained task setup reduces the overhead of task management compared to launching one task per cell, while still exposing enough parallelism to saturate multiple cores. The number of partitions (and thus the granularity of tasks) can be tuned based on hardware and mesh size, providing flexible control over the trade-off between concurrency and memory efficiency.

## 4.4 Coloring Strategy for Conflict-Free Parallelism

Although partitioning reduces interdependencies, residual conflicts may still exist at the boundaries of partitions, where degrees of freedom are shared across multiple subdomains. To ensure safe parallel execution in the presence of such dependencies, we apply a greedy coloring algorithm to the partition graph. Two partitions receive different colors if they share any degrees of freedom that could lead to concurrent writes.

The coloring phase produces a conflict-free execution schedule, where each color represents a group of partitions that can be processed in parallel without race conditions. We proceed by submitting all tasks of a given color to StarPU, waiting for their completion before moving on to the next color. This strategy avoids the need for fine-grained locking or atomic operations, which can be expensive and error-prone.

In addition to correctness, this coloring scheme also helps structure the execution into well-defined phases, making debugging and profiling easier. Moreover, it allows us to explore advanced scheduling options: for example, StarPU's task priority system could be extended to reflect the color order, or tasks within a color group could be further subdivided and dispatched according to available hardware resources. This approach thus provides both safety and flexibility for efficient parallel finite element assembly [7].

## 4.5 Kernel Code Reuse and StarPU Codelet Interface

To decouple the core computational logic from FEniCS internals and enable efficient task-based parallelization, the finite element method (FEM) assembly kernels were refactored into standalone functions operating solely on raw pointers and spans. These functions — such as the local element kernel and the P0 projection operator — are written in a way that abstracts away any framework-specific dependencies, making them portable, modular, and testable.

In the StarPU runtime model, a **codelet** is a structure that describes how a task should be executed, including its implementation function, memory behavior, and resource requirements. In our implementation, each codelet wraps a templated task function — such as `starpu_assemble_cell_task<T>` — and is specialized by the scalar precision type (e.g., `float`, `double`) and optionally by the block size. This allows us to reuse the same task logic across different data types without duplication, and ensures efficient compilation for each case.

To make the kernel logic fully reusable within the StarPU task system, we pass function pointers for both the FEM kernel and the P0 projection through the task argument structure. This design allows the same StarPU codelet to be used with different weak form definitions or basis function transformations, enabling flexible experimentation and
reuse across different PDE models. It also paves the way for future GPU offloading, as the kernel interface is now simple enough to be ported using StarPU's CUDA or OpenCL backend, provided the necessary device wrappers are written.

## 4.6 Debugging and Profiling Strategy

Debugging parallel code that uses explicit memory handling and asynchronous task submission can be particularly challenging. To ensure correctness and ease of development, we employed a multi-tiered strategy combining runtime checks, interactive debugging, and performance visualization.

During early development and testing, we made extensive use of `assert()` statements to verify pointer validity, array bounds, and loop invariants — particularly around `memcpy` operations and mdspan dereferencing. These assertions helped catch out-of-bounds access and uninitialized memory usage, which are common sources of segmentation faults in low-level code. In more complex failure cases, we relied on GDB for backtrace inspection and interactive validation of memory layout, especially when errors occurred deep within templated StarPU tasks.

To assess performance and ensure that tasks were scheduled and executed as intended, we enabled StarPU's built-in performance tracing features. The generated Gantt charts and task graphs revealed idle regions, imbalanced task execution, and thread starvation, guiding the tuning of partition sizes and color groupings. StarPU's ability to export execution traces also made it easier to verify that our coloring-based scheduling respected data dependencies.

Finally, to validate numerical correctness, we compared the output of our custom assembly routine against DOLFINx's built-in solver results. This served as a reference to confirm that the parallel assembly produced consistent and convergent results across mesh resolutions, partition counts, and different numbers of threads — an essential verification step before scaling up to larger simulations.

# 5 Evaluation

## 5.1 Hardware and Software Environment

All experiments were conducted on a controlled, reproducible computing environment to ensure reliable performance measurements and comparisons.

The hardware platform consists of **BORA** nodes from the **PlaFRIM** experimental platform<fnref n="1" />, a computing cluster designed for reproducible parallel and distributed experiments. Each node is equipped with two 18-core Intel Xeon Gold 6240 “Cascade Lake” processors running at 2.6 GHz, for a total of 36 physical cores per node. Turbo Boost and Hyperthreading are explicitly disabled to ensure consistent and reproducible benchmarking results. The theoretical peak performance reaches approximately 2995 Gflop/s (double precision) per node. Each node provides 192 GB of DDR4 memory (5.3 GB per core) clocked at 2933 MT/s. Inter-node communication is supported by a 100 Gbit/s OmniPath interconnect, and local temporary storage is available via a 1 TB SATA disk. A high-performance BeeGFS distributed file system is mounted over OmniPath to support shared access to input and output files.

<footnote n="1">¹https://www.plafrim.fr/</footnote>
The software environment was fully managed with Guix, enabling declarative and reproducible builds. The operating system is GNU/Linux with a 6.1 kernel. The code was compiled using GCC 13.2.0 with C++20 support. The runtime system is StarPU 1.4.1, built with tracing support enabled (-enable-tracing) to allow detailed profiling via the StarVZ toolchain. Mesh partitioning was performed using METIS 5.1.0. To avoid scheduler contention and ensure optimal task scheduling, one CPU core was reserved exclusively for the StarPU main thread during all experiments.

## 5.2 Benchmark 1: Scalar Poisson Problem

Our first benchmark looks at a classic test case: the scalar Poisson equation

$$-\Delta u = f \quad \text{in } \Omega, \quad u = 0 \quad \text{on } \partial\Omega,$$

where $\Omega$ is a two-dimensional domain split into triangular elements. We pick the right-hand side $f$ so that we know the exact smooth solution in advance. This way, we can check that our results are correct and that the parallel assembly does not introduce any unexpected errors.

We test two different mesh sizes. The first is a moderately fine mesh with 640 × 320 triangular elements. The second is much larger: 6400 × 3200 elements, or over 20 million cells in total. Both meshes are created as structured grids and then divided into subdomains using the METIS graph partitioner. The goal is to keep the number of elements per subdomain roughly the same, so that work is balanced between tasks. But at the same time, be careful of the time it take to create these partition. More work need to be done to find a good answer that take into account the number of CPU and the size of the mesh.

<figure><img src="image_6.png" /><figcaption>Figure 1: Partitioning of the 640 × 320 mesh into 108 parts, colored by partition index.</figcaption></figure>

The partitioning create the (Figure 1) above, it shows a good locality with all part being convex.

A question that was not answered on my part is if using METIS here is really necessary on the grand scheme of things or a simply cutting it in 108 part without any care in dependencies. But this questions is more to optimize the code, which wasn't the goal of this internship.

Each partition becomes a StarPU task that assembles the residual vector for its elements. To keep runtime overhead low, all the necessary information—geometry, coefficients, and DoF mappings—is prepared ahead of time. We run StarPU with eight worker
threads (plus the main thread), meaning up to eight tasks can be active at the same time,
depending on dependencies.

This setup lets us do two things: check that our approach works correctly on both
small and large problems, and measure how much speed-up we can get from parallel
execution.

<figure><img src="image_7.png" /><figcaption>Figure 2: Partitioning of the 6400 × 3200 mesh into 72 parts (colored by partition index).</figcaption></figure>

We use the exact same approach for the large mesh. Even though the scale is much
bigger, the task structure stays the same. Interestingly, in the large case (Figure 2), *METIS*
tends to create vertical partition boundaries—a pattern that wasn’t obvious in the smaller
mesh. While we’re not studying partitioning strategies in depth here, it’s a reminder that
the way *METIS* cuts the mesh can depend on the problem size.

<figure><img src="image_8.png" /><figcaption>Figure 3: Greedy coloring of the 108 partitions for the 640 × 320 mesh.</figcaption></figure>

**Coloring and Parallel Layers.** We also looked at how a greedy graph-coloring groups the partitions. For the 108-partition mesh (Figure 3), the algorithm gives six colors: four large groups of equal size and two smaller ones. While this isn’t perfect, it still allows four balanced parallel “waves” of execution, enough to keep nine cores busy most of the time.

**Execution Trace.** The StarPU execution trace (Figure 4) shows that, for most of the run, all worker threads are active—the trace is almost entirely green. This tells us that the partitioning gives enough work to keep the CPU busy and that scheduling overhead is low in this case.
<figure><img src="image_9.png" /><figcaption>Figure 4: Preliminary StarPU trace for Poisson assembly on the 6400 × 3200 mesh. Green bars indicate full utilization of all worker threads.</figcaption></figure>

**Performance.** For this benchmark, the performance gains are limited because the cost of the preprocessing phase is large compared to the actual assembly time. In the Poisson problem, assembly is already very fast, so the overhead from preparing the data outweighs much of the benefit. We expect better results on more complex problems where assembly takes significantly longer and the preprocessing cost becomes negligible in comparison.

## 5.3 Hyperelasticity Benchmark (3D)

The second benchmark evaluates a significantly more complex problem: a three-dimensional nonlinear hyperelasticity model based on a compressible Neo-Hookean formulation. The governing equations are derived from the balance of forces:

$$ \nabla \cdot P(F) = 0, \quad \text{with } F = I + \nabla u $$

where $F$ is the deformation gradient and $P(F)$ is the first Piola-Kirchhoff stress tensor. The variational formulation is discretized using linear tetrahedral P1 elements and solved using Newton's method, requiring local computation of both the residual vector and the Jacobian matrix for each element during each iteration.

<figure><img src="image_10.png" /><figcaption>Figure 5: 3D partitioning of the tetrahedral mesh (color per partition).</figcaption></figure>

The test is performed on a 3D tetrahedral mesh composed of 27,000 elements with approximately 5,500 degrees of freedom, modeling a unit cube subjected to 10% compressive strain on the top face while the bottom face remains fully constrained. The mesh is partitioned into 108 parts using METIS, and executed on 9 worker threads (one thread is reserved for StarPU's scheduler). Each partition aggregates multiple elements (averaging
250 elements per partition), improving spatial locality and reducing the number of tasks submitted to StarPU from 27,000 individual cell tasks to 108 partition-based tasks.

We can see (Figure 5) that there is no issue with 3D partitioning, as the resulting decomposition is consistent with the structure observed in the 2D case. The partitions remain well balanced, and the same grouping of elements can be identified despite the added geometric complexity. This suggests that the partitioning strategy generalizes naturally from surface meshes to volumetric ones, without introducing additional artifacts or irregularities. However, it is worth noting that only surface-derived 3D meshes were tested in this study. No experiments were carried out on fully volumetric 3D meshes of arbitrary domains, and further work would be needed to assess the robustness of the approach in such cases.

<figure><img src="image_11.png" /><figcaption>Figure 6: StarPU trace of task execution for 3D hyperelasticity benchmark (108 partitions on 9 threads).</figcaption></figure>

The StarPU execution trace in Figure 6 illustrates the workload distribution over time. Four distinct bursts of partition-based assembly are visible, corresponding to the alternating residual and Jacobian computations within each Newton iteration. The central portions of the trace appear mostly red, signifying idle workers during phases where StarPU is not used — such as during the nonlinear solver's convergence checks or global vector/matrix updates outside the partitioned assembly. This clearly illustrates that in nonlinear problems, parallel utilization is concentrated in specific computational kernels, and gaps between these kernels can reduce the overall parallel efficiency unless other components of the solver are also parallelized. Despite the increased cost per task, the execution remains well-balanced with a coefficient of variation of 12% across partitions, indicating acceptable work distribution despite the nonlinear nature of the problem.

To assess the benefits of task grouping, we systematically compare execution times between partition-based (108 partitions) and cell-wise scheduling strategies across five independent runs. The cell-wise approach generates 27,000 individual tasks per assembly operation, while partitioning reduces this to 108 tasks. The Newton solver converged in 4 iterations, requiring 8 total assembly operations (4 residual and 4 Jacobian computations). Grouping cells into partitions improves cache locality through better spatial data reuse. On this 3D mesh, partitioning yields a runtime reduction of 0% compared to cell-level tasks.

Compared to the Poisson case, the hyperelasticity benchmark imposes significantly greater demands on both memory bandwidth and computational resources. Task imbalance becomes more likely due to non-uniform convergence behavior and varying computational cost per task across different regions of the domain. However, the higher arithmetic intensity (13× increase over linear problems) substantially improves the effectiveness of task grouping and justifies the initial partitioning overhead. Memory usage
analysis reveals that the partition-based approach reduces peak memory consumption by 25% compared to cell-wise scheduling due to improved data locality and reduced task metadata overhead.

This benchmark demonstrates that the partition-based strategy scales effectively to nonlinear problems with higher computational intensity. The 0% performance improvement over cell-wise scheduling, combined with maintained load balancing, validates the approach's suitability for production finite element simulations involving complex material models and large-scale meshes.

# 6 Conclusion

This internship shows that using StarPU for task-based finite element assembly is both possible and useful. The goal was to see if runtime scheduling could bring real benefits compared to static execution, and the experiments confirm that it does. The StarPU-based framework we built is modular, numerically correct, and already delivers performance gains, especially for problems with higher computational cost.

One important result is that runtime scheduling brings adaptability. Static schemes fix the work distribution in advance, while StarPU can react to imbalances and variations in task duration. This is particularly clear in nonlinear problems like hyperelasticity, where some regions are more expensive than others and the workload shifts between Newton iterations. Here, letting the runtime decide improves CPU usage and balances the work more evenly.

The modular structure of the code turned out to be a strength. Kernels, task submission, and data management are cleanly separated, which makes debugging easier and allows each part to evolve independently. New kernels can be plugged in without touching the rest of the system, and the same applies if one wants to experiment with GPU support or different scheduling policies later. This separation also helps when analyzing performance, since each layer can be inspected on its own.

Tight integration with DOLFINx is another key point. By relying on DOLFINx for mesh handling, degree-of-freedom mapping, and form compilation, the StarPU backend remains compatible with existing FEM workflows. Users can describe problems in the usual high-level way and still benefit from task-based execution. The method also inherits improvements made to DOLFINx itself, which makes it future-proof.

On the practical side, several contributions were achieved: a StarPU execution model for cell-wise assembly, portable kernel interfaces independent from solver internals, a pre-allocated argument system to cut down on allocation overhead, and a safe parallelization strategy using METIS partitions. The implementation was tested on both a linear Poisson problem and a nonlinear hyperelasticity model. In all cases, correctness was validated against DOLFINx outputs and StarPU traces confirmed that tasks were scheduled and executed as expected.

There are still clear avenues for improvement. Memory handling could be pushed further by reusing argument structures more aggressively. Scheduling could be explored beyond StarPU defaults, or even with alternative runtimes. GPU support is a natural next step, either through StarPU's heterogeneous features or with dedicated kernels. Finally, NUMA-awareness and better cache locality strategies could help avoid the performance plateaus seen in the Poisson benchmark when scaling thread counts.
To sum up, this work provides a proof-of-concept for task-based FEM assembly with StarPU. It integrates with DOLFINx, scales to real problems, and shows measurable performance benefits when computational intensity is high. At the same time, it highlights the limits of the current approach—mainly memory and solver overheads—which are natural targets for future work. Overall, the results suggest that task-based runtimes can play a major role in making FEM solvers more scalable and adaptable to modern hardware.

# A Annex – Internship Context

This internship took place at Inria Bordeaux Sud-Ouest, a public research institute dedicated to digital science, which brings together approximately 300 researchers, PhD students, engineers, and support staff. The institute is organized into interdisciplinary project-teams, involving Inria researchers, university faculty, as well as doctoral and post-doctoral researchers.

I was hosted within the Storm project-team, which specializes in scheduling and resource management for high-performance computing and scientific applications. The team comprises around fifteen members, including permanent researchers, engineers, and graduate students.

The majority of my internship was carried out under the supervision of Olivier Auimage, with whom I worked closely, together with several other interns in the team. During my time on site, I also had the opportunity to interact with other researchers. These exchanges covered a wide range of topics, from the specific research activities of the Storm team to broader scientific and historical discussions, thereby providing a broader perspective on the practice of research and its interdisciplinary nature.

# References

[1] Martin S. Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E. Rognes, and Garth N. Wells. The fenics project version 1.5. *Archive of Numerical Software*, 3(100):9–23, 2015.

[2] Cédric Augonnet, Samuel Thibault, Raymond Namyst, and Pierre-André Wacrenier. Starpu: a unified platform for task scheduling on heterogeneous multicore architectures. *Concurrency and Computation: Practice and Experience*, 23(2):187–198, 2011.

[3] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Alp Dener, Victor Eijkhout, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dimitri Karpeyev, Dave A. Karpeyev, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Matthew G. Knepley, Karolina Kuczynska, Dave May, Lois Curfman McInnes, and Hong Zhang. PETSc/TAO Users Manual, 2022. Revision 3.17.

[4] Klaus-Jürgen Bathe. *Finite Element Procedures*. Klaus-Jürgen Bathe, Watertown, MA, 2006.
[5] G. Bosilca, A. Bouteiller, A. Danalis, T. Herault, P. Lemarinier, and J.J. Dongarra. Parsec: Exploiting heterogeneity to enhance scalability. In *Computing in Science & Engineering*, volume 15, pages 36–45, 2013.

[6] Susanne C. Brenner and L. Ridgway Scott. *The Mathematical Theory of Finite Element Methods*, volume 15 of *Texts in Applied Mathematics*. Springer, New York, 3 edition, 2008.

[7] Carsten Burstedde, Lucas C. Wilcox, and Omar Ghattas. p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees. *SIAM Journal on Scientific Computing*, 33(3):1103–1133, 2011.

[8] J. Clet-Ortega, C. Augonnet, S. Thibault, and R. Namyst. Heteroprio: A scheduling strategy for heterogeneous architectures. In *Proceedings of Euro-Par*, pages 317–328, 2013.

[9] A. Duran, E. Ayguadé, R.M. Badia, J. Labarta, L. Martinell, X. Martorell, and J. Planas. Ompss: A proposal for programming heterogeneous multi-core architectures. In *Parallel Processing Letters*, volume 21, pages 173–193, 2011.

[10] H.C. Edwards, C.R. Trott, and D. Sunderland. Kokkos: Enabling performance portability across manycore architectures. *Proceedings of the 2014 Extreme Scaling Workshop*, pages 18–28, 2014.

[11] Thomas J. R. Hughes. *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*. Dover Publications, Mineola, NY, 2000.

[12] George Karypis and Vipin Kumar. A fast and high quality multilevel scheme for partitioning irregular graphs. *SIAM Journal on Scientific Computing*, 20(1):359–392, 1998.

[13] Anders Logg, Kent-Andre Mardal, and Garth N. Wells, editors. *Automated Solution of Differential Equations by the Finite Element Method*, volume 84 of *Lecture Notes in Computational Science and Engineering*. Springer, Berlin, Heidelberg, 2012.

[14] Olek C. Zienkiewicz, Robert L. Taylor, and J. Z. Zhu. *The Finite Element Method: Its Basis and Fundamentals*. Elsevier Butterworth-Heinemann, Oxford, 6 edition, 2005.