# High-performance cell-based simulations of cardiac electrophysiology using task parallelism

Giorgio Bettonte

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

## To cite this version:

Giorgio Bettonte. High-performance cell-based simulations of cardiac electrophysiology using task parallelism. Distributed, Parallel, and Cluster Computing [cs.DC]. 2025. hal-05224224

HAL Id: hal-05224224

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

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>

UNIVERSITÉ DU
LUXEMBOURG

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

HIGH-PERFORMANCE CELL-BASED SIMULATIONS OF
CARDIAC ELECTROPHYSIOLOGY USING TASK
PARALLELISM

Author

Giorgio Bettonte

Supervised by

Dr. Olivier Aumage, Inria researcher

Reviewed by

Dr. Pierre Talbot, University of Luxembourg

Examined by

Dr. Georgios Kafanas, University of Luxembourg

A Thesis submitted in fulfillment of requirements for the degree of

**Master in High Performance Computing**

Department of Computer Science
University of Luxembourg
2025
<figure><img src="image_5.png" /></figure>
<figure><img src="image_6.png" /></figure>

UNIVERSITÉ DU
LUXEMBOURG

# Declaration of authorship

For single-authored work:

I hereby declare that I am the sole author of the work entitled

HIGH-PERFORMANCE CELL-BASED SIMULATIONS OF CARDIAC ELECTROPHYSIOLOGY USING TASK PARALLELISM

and here enclosed, and that I have compiled it in my own words, that I have not used any other than the cited
sources and aids, and that all parts of this work, which I have adopted from other sources, are acknowledged
and designated as such. I also confirm that this work has not been submitted previously or elsewhere.

17/08/2025

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

Date

Signature

Giorgio Bettonte

First name and surname
<figure><img src="image_8.png" /></figure>
**Abstract**

Cardiac electrophysiology is the study of electrical activity in the heart; disorders of the heart's electrical synchronisation system, named cardiac arrhythmias, are responsible for half of the deaths related to cardiovascular disease. Cardiac electrophysiology relies on computer models that simulate the behaviour of the heart and how it is affected by cardiac arrhythmias.

The existing software are already very sophisticated and valid tools, but not yet powerful to take into account all the individual cardiac cells in a human heart, which are around 2 billion; instead, they group them and assume the same behaviour for all the cells belonging to the same group. This limitation does not allow one to describe precisely the interaction between neighbouring cells, which is fundamental for a better understanding of the cardiac diseases.

Finer models, called cell-based models would be able to provide more accurate simulations, but only with a relevant backward in terms of computational cost. Due to their complexity, cell-based softwares need to use large computer architectures, called supercomputers. Although, using these platforms is not trivial: specific competences are required to achieve efficient performance and the variability between platforms leads to portability issues. To deal with these issues it can be convenient to rely on runtime systems, libraries handling the executions of programs; if well designed, a runtime system can improve the performance of a software and extend its portability, while requiring moderate efforts by the programmer.

The challenges represented by cell-based models are the focus of many research efforts; this work is part of them, and specifically it investigated the possible benefits for cardiac cell models derived by the use of a runtime system. This document describes the introductory work made to combine a cardiac cell-based model application, DEMI software, with a runtime system, StarPU. The used approach will be illustrated, and a perspective on how it could be possible to extend the work will be provided as well.
<figure><img src="image_9.png" /></figure>
## Plain Language Summary

Having a good understanding of the heart and related diseases is necessary for the medicine. Today, computers are involved in the scientific research dedicated to cardiology, and programs have led to many advances in this field. Unfortunately, many questions are still open and there is the need to adapt our softwares. In this context, I worked on a cardiac application and tried to adapt it to tackle large computer infrastructures, called supercomputers. Supercomputers are complex computer systems, and they can be very different between them. This results in higher difficulties when the same application is needed on different supercomputers. This problem is defined *portability issue* and was one of the main concerns of the work.
<figure><img src="image_10.png" /></figure>
## Acknowledgments

This work was supported by the MICROCARD-2 project (project ID 101172576). The project is supported by the EuroHPC Joint Undertaking and its members (including top-up funding by ANR, BMBF, and Ministero dello sviluppo economico). Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or EuroHPC. Neither the European Union nor EuroHPC can be held responsible for them. The activity was conducted at the INRIA laboratory of Bordeaux [Bora], as part of the on-going collaboration between the team STORM [INR25] and the SIMULA laboratory center [cen25b].

I would like to thank James Trotter, developer of DEMI, for his guidance and availability; I also thank Professors Pierre Talbot and Gianluca Palermo for accepting to be my academic supervisors and for their valuable advice. Finally, a special thanks goes to my company supervisor, Olivier Aumage, for his thoughtful mentorship throughout my internship and for his constant support.

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

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

EuroHPC
Joint Undertaking

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

Co-funded by
the European Union

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

STORM
STATIC OPTIMIZATIONS - RUNTIME METHODS

simula
<figure><img src="image_15.png" /></figure>
# Contents

<table><tr><td>Declaration of Authorship</td><td>iii</td></tr><tr><td>Abstract</td><td>iv</td></tr><tr><td>Plain Language Summary</td><td>vi</td></tr><tr><td>Acknowledgements</td><td>viii</td></tr><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Background</td><td>3</td></tr><tr><td>2.1 Computer models and current state of the art for heart simulation</td><td>3</td></tr><tr><td>2.2 Cell based approach</td><td>5</td></tr><tr><td>2.2.1 EMI model</td><td>5</td></tr><tr><td>2.2.2 DEMI software</td><td>5</td></tr><tr><td>2.2.3 High computational costs</td><td>5</td></tr><tr><td>2.3 HPC</td><td>6</td></tr><tr><td>2.3.1 Programmability Enhancements Issue</td><td>7</td></tr><tr><td>2.3.2 Performance Portability Issue</td><td>7</td></tr><tr><td>2.4 Runtime systems</td><td>7</td></tr><tr><td>2.4.1 Task-based runtime systems</td><td>8</td></tr><tr><td>2.4.2 STARPU</td><td>8</td></tr><tr><td>2.5 The project: recap and aim</td><td>13</td></tr><tr><td>3 Contribution</td><td>14</td></tr><tr><td>3.1 FHN and TP06 models</td><td>15</td></tr><tr><td>3.1.1 Involved structures and registration</td><td>16</td></tr><tr><td>3.1.2 Partitioning of handles, Euler method codelet and insertion of tasks</td><td>17</td></tr><tr><td>3.1.3 DAG for FHN</td><td>18</td></tr><tr><td>3.1.4 ViTE traces FHN model</td><td>18</td></tr><tr><td>3.1.5 Time performance analysis FHN model</td><td>20</td></tr><tr><td>3.2 DIFFUSION model</td><td>24</td></tr><tr><td>3.2.1 Datatype definitions</td><td>25</td></tr><tr><td>3.2.2 CG on single node</td><td>28</td></tr><tr><td>3.2.3 Data partitioning and distribution</td><td>31</td></tr><tr><td>3.2.4 CG multiple nodes</td><td>32</td></tr><tr><td>3.2.5 ViTE analysis for CG multiple node</td><td>34</td></tr><tr><td>3.2.6 Time performance analysis DIFFUSION model</td><td>35</td></tr><tr><td>4 Conclusions</td><td>39</td></tr><tr><td>References</td><td>41</td></tr><tr><td>5 Appendices</td><td>43</td></tr><tr><td>5.1 DEMI structures</td><td>43</td></tr></table>
# 1 Introduction

Cardiac electrophysiology is the study of electrical activity in the heart. It is of crucial importance, as disorders of the electrical synchronisation system of the heart, called cardiac arrhythmias, are responsible for half of deaths related to cardiovascular disease.

Today, electrophysiology uses computer models to simulate and describe how the electric potential is propagated in the human heart. The existing software are already very sophisticated and valid tools, but not yet powerful to take into account all the individual cardiac cells, which are around 2 billions; instead, they group them and assume the same behaviour for all the cells belonging to the same group. This aspect represents a limitation: current models can not provide extensive replication of some events determined by reduced electrical coupling between cells and related to cardiac disease. Therefore, there is a need for models in which individual cells and their connections are taken into account.

These kinds of models are called cell-based models. Thanks to a finer modelling, they are able to provide a very accurate description at the cellular level, and they have been considered to be more suitable than classical models for some specific cases.

However, due to the high number of cells present in the heart, simulating with a cell-based approach results in a extreme amount of computation, representing a considerable obstacle. Supercomputers can represent a solution; supercomputers are special machines designed to perform massive computations with high efficiency, and they are largely used for scientific applications. The related computer science field is called High Performance Computing (HPC). However, these architectures come with new challenges. HPC infrastructures are large complex environments, where often different types of hardware live together; code development is harder and requires specific competences. Moreover, among HPC platforms, there is not uniformity, and they change constantly: therefore, necessary continuous adjustments to ensure portability are necessary.

Taking into account all these aspects, the use of HPC resources is not trivial. To address all these HPC challenge, there exist means and techniques; runtime systems are among them: these libraries handle the execution of a program, ensuring efficient performance and portability. There are different possible runtime paradigms, *task based* was adopted for this work. Tasks are basically an amount of work to do and can be considered as functions in the parallel context. Tasks are the basis concept for task-based runtime systems.

Task-Runner time systems have already shown their capabilities when used to manage scientific applications, and they can be beneficial for cell models as well. This work aims to explore this possibility and has specifically focused on the integration of a cell-based software, DEMI software, and a task-based runtime system, StarPU. Specifically, my work consisted of an initial combination with StarPU, and it represents the first attempt in this sense for DEMI. Due to the size of this scientific application, trying to do this integration in one big run would be impracticable; instead, an incremental approach was adopted. DEMI is composed of many components, each of them being a model on its own, some simpler than others. These submodels cooperate, and in some cases are built on top of each other. For this work, starting from the simplest ones, I handled successively three submodels of DEMI; this means, think using a "task" logic approach and redesign the code accordingly. The long-term goal will be to obtain a unified code version for DEMI, able to execute on different heterogeneous architectures efficiently.

In Section 2, I will give more information about the context of the field, the state of the art, and the instruments mentioned above. In Section 3, I will present what I specifically did during the project. Finally, in Section 4, I will discuss what it has been done and how the work can continue.
**Contributions.** The contributions of this work are summarized as follows:

1. Partial porting of the transmembrane models of the DEMI software (Section 3.1);

2. Partial porting of the diffusion model of the DEMI software (Section 3.2).
## 2 Background

Cardiac electrophysiology is a branch of medicine that focusses on the electrical activities of the heart. The conduction of electricity in the human heart is fundamental: the pumping function depends on the contraction of the heart, which is initiated by an electrical signal spreading through the cardiac muscle. It is well established that abnormalities in cardiac conduction are associated with fatal consequences[Jæg+19]; disorders of the heart's electrical synchronisation system, called cardiac arrhythmias, are responsible for half of deaths related to cardiovascular disease [pro]. The knowledge about the electric conduction in the human heart has evolved in the last decades and intense research, based on experiments, mathematical models, and computer simulation has led to important results. Still, many questions are still open, and current computers models seem to be unable to address them. Many computer models, due to simplistic assumptions, lack in accuracy; therefore, other cardiac models, named cell-based cardiac models, were proposed to fill this gap. Unfortunately, the higher provided precision comes with a high incremental of computational need. Special infrastructures, named supercomputers, can provide the resources necessary for cell-based models. Supercomputers, or High Performance Computing (HPC) platforms, are complex and varied systems. To deal with them, libraries called runtime systems are exploited. The combination of cell-based models, HPC power and runtime systems techniques could contribute to the progress of cardiac electrophysiology, and it is the focus of many research efforts; this work is one of them. This synergy was investigated working in the integration of a scientific application, named DEMI, implementing a cell-based approach, with a runtime system StarPU, with particular attention to portability needs. DEMI is developed by the HPC department of the Simula laboratory of Oslo [cen25b]; StarPU is developed by the team STORM [INR25], at the INRIA laboratory centre of Bordeaux [Bora] [Borf] [LaB25]

In this background section, we will discuss the fundamental background concepts of this work; what are computer models and the state of the art 2.1.; why the cell-based approach is interesting 2.2, a possible implementation for it 2.2.2 and the related difficulties 2.2.3. We will also talk about the HPC means used: why they are used 2.3 and what their challenges are 2.3.2 2.3.1; how particular libraries, called runtime systems, can be useful 2.4, and what the basis concepts of one of them, StarPU 2.4.2. Finally, in 2.5 what is the goal of work and its motivations.

### 2.1 Computer models and current state of the art for heart simulation

Science investigation has always based on real life experiments to validate hypothesis and formulate conclusions. This is not changed; observations and measurements are basic block for every real scientific study. Although, design and conduct experiments is not trivial: every real life experiment requires time, appropriate setup, physical resources and reproducibility. This can be challenging, extremely expensive and consequently a obstacle for the research. Fortunately, nowadays it is possible to rely on another important resource: computer models. Software simulations can complement and improve experiments; they allow to replicate virtually phenomena, study their behaviour and design further real-life experiments more consciousness. Computer models are therefore largely used in many scientific domains and they are the results of competencies coming from mathematics, physics, computer sciences and domain knowledge.

A phenomena needs to be studied and carefully observed; after that, mathematics is deployed to formalize its description. The resulting continuous mathematical model, is then dis-
cretized through *numerical modelling*; indeed, to be actually solved my machine, problems need to be finite representable. At this point *scientific computing*, through the use of a calculator, is in charge of finding an approximate solution for the original problem. Finally, the solution is evaluated, compared with the observations and the following activities are repeated. Figure 1 illustrates this circular workflow.

<figure><img src="image_16.png" /><figcaption>Figure 1: Circular workflow for the development of scientific code. Image taken from [Sun+07]</figcaption></figure>

This development procedure is adopted mainly for every scientific applications, included those ones that are used to simulate the behaviour of the human heart. In particular, we are interested in software useful to cardiac electrophysiology, which a branch of the medicine dedicated to the study of the electrical activity in the heart. The conduction of electricity in the human heart is fundamental: the pumping function depends on the contraction of the heart, which is initiated by an electrical signal spreading through the cardiac muscle. It is well established that abnormalities in cardiac conduction are associated with fatal consequences[Jæg+19]; disorders of the heart's electrical synchronisation system, named cardiac arrhythmias, are responsible for half of the deaths related to cardiovascular disease [pro].

Thanks to an intense research, cardiac electrophysiology has achieved relevant results. Cardiac computer models has played a crucial role and they have lead to a better comprehension of heart disease. Although, many questions are still open and unfortunately current cardiac models seems to not be sufficient to address some of them. Classical cardiac computer models don't take into account all the individual cardiac cells, instead, they group them and assume the same behaviour for all the cells belonging to the same group. This simplistic assumption work well in many cases, but in others has shown its limited capabilities; this modelling can not allow an extensive replication of some events related to cardiac disease and determined by reduced electrical coupling at cell-level.
## 2.2 Cell based approach

To overcome the inner limitations of the classic approach, we need models in which the individual cells and their connections are taken into account: cell-based models. [pro] There exist various cell-based models, but for this project the EMI model was considered.

### 2.2.1 EMI model

We refer to the model as the EMI model because it includes an explicit representation of the Extracellular space, the cell Membrane and the Intracellular space [Jæg+19], while in more classical approaches, there is an homogenization of the tissue and the resulting models therefore assume that the three spaces exist everywhere [Tve+17] and they are not distinguished. Figure 2, shows an EMI scenario for a simple case involving two cardiac cells.

<figure><img src="image_17.png" /><figcaption>Figure 2: EMI configuration for two cardiac cells. Image taken from [CT24]</figcaption></figure>

The feasibility of the approach has been investigated [Jæg+24] [Tve+17] and the model has been applied in different research articles [Rei+23] [Jæg+19], where it has shown its capabilities.

### 2.2.2 DEMI software

Again, even for the EMI model different implementations are possible; for this project, I worked with DEMI software, a free distributed-memory parallel solver [cen25a], written in C language, developed by SIMULA research centre. [cen25b]

The solver employs a temporal operator splitting and an iterative, non-overlapping domain decomposition. First, the temporal splitting is used to decouple the nonlinear cell model and membrane dynamics from the rest of the equations. Second, the domain decomposition is applied to solve the equations for the intra- and extracellular domains separately, [cen25a] using the Conjugate Gradient (CG) method.

DEMI offered working different and separate code versions including programming parallel paradigm such as OpenMP, MPI and CUDA. Moreover, DEMI can be seen as a composition of subparts, actual models on their own, built one on top of each other.

The large variety of DEMI represents a strength, but can as well lead to portability issues; we will address more in details what does this terms mean in 2.3.2.

### 2.2.3 High computational costs

Although the EMI model emerges as the only viable choice in some specific cases [Jæg+24], the greater accuracy offered by the model comes at higher computational cost. In [Jæg+24], different cardiac models software implementations were compared and evaluated.

The experiment setup was the following: 13x65 cardiac cells problem for a Dell Precision 3640 Tower with an Intel Core processor (i9-10900K, 3.7 GHz/5.4 GHz) with ten cores with two threads each.
<figure><img src="image_18.png" /><figcaption>Figure 3: Comparison between different cardiac models. CPU time for EMI is much larger compared to the others; this reflects the higher complexity of the model. Image taken from [Jæg+24]</figcaption></figure>

Figure 3, taken from the article, shows the result of this experiment.

Leaving out the medical considerations and looking at the CPU time, it is clear that the EMI model is much more resources consuming than the others. This result is even more alarming if we consider that the reported example was conducted on few hundreds of cells, while in the human hearts there are about 2 billions of cells. To be a feasible approach and therefore actually useful, EMI needs large computational resources.

## 2.3 HPC

Due to its huge complexity, EMI needs to use powerful and large architectures, and the field of computer science that deals with this topic is called High Performance Computing (HPC). The purpose of the HPC is to design supercomputer hardware, software tools, and programming techniques to solve relevant scientific problems of large scale.

HPC is a crucial resource in the world of scientific computing, since it allows to simulate complex phenomena. Indeed, simulations can complement, prepare, and in some cases even entirely supplant real life experiments —with substantial economical benefits as the most immediate advantage— while enabling new studies for which real experiments would be impractical. [Aum20]. However, HPC capabilities come with new challenges.
### 2.3.1 Programmability Enhancements Issue

HPC platforms, or supercomputers, are complex computing systems, and they are often referred to as "heterogeneous architectures"; indeed, they are the result of the combination of different machines with different structure and purpose. In a supercomputer, we can find general purpose CPUs, designed to address a large variety of needs, and accelerators, special units optimised for specific workloads, for example GPUs.

This integration can lead to rather complex systems, and the process of writing applicative codes for HPC inherently requires a specific expertise often beyond the skills of a scientific simulation programmer in a given application domain[Aum20]. Indeed, in a supercomputer, all the components are designed to offer a high level of parallelism to achieve better performance, but this means also combine and synchronise many different works: the code design needs to interleave executive code with control execution code and has to be appropriate for the target hardware.

### 2.3.2 Performance Portability Issue

Many of the choices made to adapt an applicative code to a platform may be changed again when the code has to be ported to another platform.

The differences among supermachines, are easily observable looking at the top500 List, which ranks and details the 500 most powerful non-distributed computer systems in the world. [Pro25]

Already in the top10, there are three powerful machine which present peculiar and different characteristics: HPC6 [HPC25], JUPITER Booster [Boo25] and Supercomputer Fugaku. [Fug25]: the first two are both engaged with GPUs, but from different vendors, while Fugaku does not have any GPU. This all means, for an applicative code: not uniformity and not compatibility.

Moreover, it is important to consider that, further the mutual divergencies, platforms evolve and change rapidly: the TOP500 is updated every 6 months and that is not surprising to see significant changes from a ranking to the next one.

Variability among platforms and constant evolution of the technologies require from HPC software to have a property called performance portability, that is, the ability for the resulting program to run efficiently on a variety of hardware technologies with only moderate adaptation redevelopment costs. To enforce this property, compilers can offers static performance portability: they are able to offer to the programmer an high level of abstraction, while maintaining optimized performance for the target architecture. [Aum20]

However, a static approach can not be not sufficient for modern HPC environment, due to their complexity. [Aum20] This limitation derives mostly by their lack of adaptivity at runtime. Indeed, a priori evaluation can not take into account aspects like variability of computer resources' performance, the combinatorial explosion caused by the number of available computing resources/memory-hierarchies/devices, some inherent irregularity in the application algorithm, or even some system noise. [Aum20]

For all these reasons, runtime systems can be interesting and useful.

## 2.4 Runtime systems

Runtime systems are software layers linked with an application and/or with numerical libraries, whose responsibility is to drive the processing of work requests expressed by the ap-
plication or libraries into some hardware resources [Aum20]. Even if they represent an addi-
tional layer for the computation, runtime systems are a valid helper for the programmer, who
can exploit a defined general interface avoiding taking direct decisions about the hardware
management. There are different types of runtime systems and again various implementa-
tions for each of them; we will see in detail the task-based runtime systems and the StarPU
runtime library.

**2.4.1 Task-based runtime systems**

A task-based runtime system is based on the notion of *tasks*, *task data dependencies* and *DAGs*

**Concept of task** A task can be considered as a function which operates in the context of parallel programming and is in charge of some work.

**Concept of data dependences between tasks** Different tasks can work on the same data; for this reason, they are linked and related by some data dependencies. This aspect is crucial: a task cannot be executed until all its dependences are satisfied. There are various ways to encode dependencies, for the purpose of our work we will illustrate only the *Sequential Task Flow* (STF). STF programming model consists in fully relying on sequential consistency, exclusively using implicit dependencies. [Agu+17] With this approach, tasks are submitted by non-blocking functions to the runtime system, which is in charge of the asynchronous scheduling and parallel execution.

**Concept of data DAG** A direct acyclic graph (DAG) is created. In the DAG, tasks are represented by vertices, while edges express dependencies. The graph is a logical representation of the code and its constraints and is used to implement correct scheduling.

**2.4.2 STARPU**

StarPU is a task programming library for hybrid architectures [Bore] [Bord], developed by the LaBRI laboratory [LaB25] and the Inria Research Centre [Bora] at the University of Bordeaux [Borf]. The StarPU runtime system has been designed from the start to provide a task oriented parallel programming model targeting multicore nodes of HPC platforms, optionally equipped with accelerator boards or general purpose GPU boards [Aum20]. StarPU is still an object of research and its related literature can count on several publications [Borc]. Moreover, it has already shown its capabilities, and there are some software that uses it to tackle heterogeneous architectures. Some of them are reported in [Borb].

As we said, StarPU is a task-based runtime system; therefore, it is based on the concepts addressed in the previous section and uses the STF paradigm.

In addition to the general task properties, StarPU has its own API. In the following, we discuss the main blocks of StarPU, giving some simple but concrete code examples in Listings 1 2 3 4 5. Specifically, all the combined examples show how to submit a task that sets all the elements of a vector to a given integer.

**StarPU handles** The first basic step, when using StarPU, is the registration of the interested data: the application performs this operation in order to get an handle for a data buffer and to "notify" StarPU about its existence; without this step, StarPU can not actually work on that data. An example of registration for a vector is reported in Listing 1.
<figure><img src="image_19.png" /><figcaption>Listing 1: Vector registration on StarPU</figcaption></figure>

**StarPU codelets** Then, another fundamental block is *codelet*. A codelet is a data structure associated with a task. The codelet provides to the runtime system instructions about the task to execute: for which architectures there exists a version of the kernel, how many data buffers are used and how they are accessed, etc.

All this information allows StarPU to: correctly detect all the task data dependencies, thanks to the specified access mode for a buffer, better chose between different available machine at a given time; indeed, if, for the same purpose, we have two different function implementations, but for two different architectures, specifying both in the codelet, it is possible to postpone the decision of what hardware to use at runtime, when a more conscious decision can be made by StarPU. An example of a codelet is reported in Listing 2.

<figure><img src="image_20.png" /><figcaption>Listing 2: StarPU codelet</figcaption></figure>

**StarPU tasks submission** To be executed, a task needs to be submitted to StarPU and the inputs must be coherent with what is specified in the codelet. An example of task submission is reported in Listing 3; note that the insertion example is consistent with what is shown in Listings 1 and 2.

<figure><img src="image_21.png" /><figcaption>Listing 3: StarPU task submission</figcaption></figure>

**StarPU tasks kernel** Finally, we have the task kernel implementation. An example of a task kernel is reported in Listing 4; at the internal of the tasks, it is necessary to retrieve the pointers for the handles and for the inline arguments.
<caption>Listing 4: StarPU task kernel</caption>

```c
void scal_cpu_func(void *buffers[], void *cl_args) {
    /*retrieve pointer and size for the vector*/
    struct starpu_vector_interface *vector_interface = buffers[0];
    int *vector = STARPU VECTOR GET_PTR (vector_interface);
    int size_vector = STARPU VECTOR GET_NX (vector_interface);

    /*retrieve pointer for inline argument*/
    int value;
    starpu_codelet_unpack_args (cl_arg, &value);

    /*set all the entries of the vector*/
    for (int i = 0; i < size_vector; i++)
        vector[i] = value;
}
```

**Handle partitioning** An existing piece of data can be partitioned into parts to be used by different tasks. The number of partitions can be decided statically or dynamically. Listing 5 shows a possible data partition for our example. This partitioning actually introduces real parallelism in our ongoing example.

<caption>Listing 5: Synchronous partitioning</caption>

```c
/*Declaration of vectors and variable*/
...

/*Registration on StarPU*/
...

int PARTS;
/*static or dynamic assignment for PARTS*/
...

/* Partition the vector in PARTS sub-vectors */
struct starpu_data_filter f =
{
    .filter_func = starpu_vector_filter_block,
    .nchildren = PARTS  //number of partitions
};
starpu_data_partition (vector_handle, &f);

/*Insertion of a task for each subpartition*/
for (i=0; i< starpu_data_get_nb_children (vector_handle); i++) {

    /*retrieve subpartition and insert a task working on it*/
    starpu_data_handle_t sub_handle = starpu_data_get_sub_data (handle, 1, i);
    starpu_task_insert (&set_value_cl,
                       STARPU_RW, sub_handle,
                       STARPU_VALUE, value, sizeof (value),
                       0,
                       );
}

starpu_data_unpartition ();
...
```
What is shown in Listing 5 is called synchronous partitioning: `starpu_data_partition()` and `starpu_data_unpartition()` both wait for all the tasks currently working on the data. This can be a bottleneck for the application, but it is possible to obtain an asynchronous partitioning using the functions `starpu_data_partition_plan()` and `starpu_data_unpartition_plan()` instead of the ones just mentioned.

**Task granularity** StarPU, performing its operation, introduces some overhead, which is not always negligible. This overhead is covered by efficient task scheduling and workload; therefore, the duration of tasks plays a role in the counterbalance. The good duration can differ among applications, but in general the minimum acceptable time length is around 100 µs.

To check the length of the tasks and the general behaviour of the program it is possible to use some tools that give insides on the performance. There exist various possibilities, but for this work I used ViTE, an open tool, developed in Bordeaux, which allows the visualisation of pajé traces generated by StarPU [Bor25].

**Real example with Cholesky factorization** Once we have defined these basic concepts, we are ready to generate concrete data flow submission of and understand the relative StarPU DAG. Tasks are submitted with the STF programming model and executed in the background by the StarPU scheduler. Figure 4 shows how a Cholesky factorisation for a $N \times N$ matrix A, can be implemented in StarPU (Listing 6), and the corresponding DAG for a tiled matrix $3 \times 3$. Cholesky factorisation is often used as an example for task-based runtime systems. The example is taken from [Agu+17]; for conciseness, we omit here the definitions of codelets and functions and the registration of matrix A.
<figure><img src="image_22.png" /><figcaption>Figure 5: StarPU-MPI for Cholesky factorization. Yellow regions underline the different 'views' for the MPI processes: what tasks are executing and on which nodes. A different point of perspective is underlined by the arrows in red: the same data dependency, requires for MPI rank 0 a send message, while for rank 1 a receive operation. StarPU detects the dependencies and handles the needed data transfers.</figcaption></figure>

## 2.5 The project: recap and aim

In 2.2.1, we saw how the EMI model, a cell-based cardiac model, has been appreciated for its accuracy in describing phenomena, at the cellular and even subcellular level, acting in the human heart. EMI model might allow for more complete knowledge about the transmission of electrical potential in the heart and help research against arrhythmia.

Unfortunately, even if EMI is considered promising, in many cases it is not considered a feasible approach due to its huge complexity. A cell approach requires a finer discretization and consequently higher computational costs, which would become extremely significant if we had to tackle all cells of a human heart.

HPC platforms can solve this problem, but they are complex systems, not uniform, and constantly changing. Therefore, using these resources efficiently is not trivial and it can be convenient to use runtime systems, tools designed for this purpose.

**Problem statement** Many complex scientific applications benefit from the use of runtime systems for their execution and portability, and it can also be the case for cell-based cardiac software, by unifying the multiple per-programming model codebases (sequential, OpenMP, MPI, GPU) into a single one. The main aim of this work project was to investigate this possibility; in particular, I worked on an initial integration of DEMI, a cell cardiac software, with StarPU, a task runtime system. In Section 3 I explain in detail how I conducted my work.
# 3 Contribution

This work investigates the portability of a cardiac cell-based software, DEMI, on top of a run-time system, StarPU. It represents the first attempt ever for the combination of these two software, and for this reason an approach from scratch was adopted, and only tasks executing on CPUs were addressed due to time constraints.

The three following paragraphs mention the main key aspects of the project and briefly anticipate what we will discuss more in detail in 3.1 and 3.2

**Challenges** The conversion in a series of tasks for a program not initially designed for that is not trivial. First, it is fundamental to have a deep understanding of the program workflow, which is necessary to understand how the data are managed, how they evolve during the execution, and when dependencies occur. Once we have outlined a task logical abstraction, it is time to deal with the technical development issue, represented by the classical programming challenges and by the specific runtime systems needs. In StarPU, handles (2.4.2) are of crucial importance; data need registration, intelligent partitioning, and correct distribution, therefore, are necessary. In DEMI, some of the data structures are not suitable to be managed directly with StarPU; therefore, especially for the most complex model addressed, it is necessary to define new datatypes for them.

**DEMI code structure** DEMI can be seen as a combination of different models cooperating or built one on top of each other. Even if the models present different complexity, they are all complete models on their own and, in turn, they can be composed by submodels. Therefore, in the following discussion we will interchangeably use the terms model and submodel. Figure 6, shows how the models mentioned and discussed in this document are related in DEMI.

All models present the same logic and the code structure follows the steps, and it is illustrated in Figure 7.

**Used approach and goals** The work carried out mainly focused on the solution step, which implements numerical methods. For three submodels of DEMI, these numerical methods have been redesigned to delineate StarPU tasks and express data dependencies.

Due to the size of the application of DEMI, an approach aiming to port completely all the models simultaneously would have been impractical, and therefore I decided to proceed step by step, starting with the simplest models, improving them incrementally, test their functionalities and eventually moving to other models. Since the work aimed to be the first working integration of DEMI with a runtime system, the main concern was related to ensure correctness and portability. All the results obtained with the StarPU-MPI were validated and compared with those obtainable by the original version of DEMI. Although it was not the main concern of the work, some performance analyses were performed.

In the following section, we will discuss the three submodels addressed: *demiFHN 3.1*, *demiTP06 3.1* and *demiDIFFUSION 3.2*; the first two will be presented together: they are very similar and roughly only differ by their computational weight.

**Used means** The platform used for the experiments was PlaFRIM (Plateforme Fédérative pour la Recherche en Informatique et Mathématiques) [web]. PlaFRIM is a scientific instrument located in Bordeaux (France) and designed to support experiment-driven research in all areas of applied mathematics related to modelling and high performance computing. On the platform, it is possible to load and use provided modules, while SLURM is the used job scheduler. The nodes used to execute the code are shown in Table 1.
<figure><img src="image_23.png" /><figcaption>Figure 6: Some models of DEMI and their relation. The demiMONODOMAIN model is the combination of the demiDIFFUSION model with a chosen transmembrane model.</figcaption></figure>

<table><caption>Table 1: Used nodes on Plafrim</caption><thead><tr><th colspan="2">Used nodes on Plafrim</th></tr><tr><th>Node</th><th>CPU</th></tr></thead><tbody><tr><td>bora001-004</td><td>2x18-core Intel CascadeLake</td></tr><tr><td>sirocco07-13</td><td>2x16-core Intel Broadwell</td></tr><tr><td>sirocco14-16</td><td>2x16-core Intel Skylake</td></tr><tr><td>sirocco17</td><td>2x20-core Intel Skylake</td></tr><tr><td>sirocco18-20</td><td>2x 20-core Intel CascadeLake</td></tr><tr><td>sirocco21</td><td>2x 24-core AMD Zen2</td></tr><tr><td>sirocco22-25</td><td>2x 32-core AMD Zen3</td></tr></tbody></table>

## 3.1 FHN and TP06 models

DEMI includes in its internal subsections demiTP06 and demiFHN. demiTP06 is used to implement the Tusscher-Panfilov (2006) (TP06) model [TP06]; demiFHN to implement the FitzHugh-Nagumo (FHN) model [Sun+07]. They are nonlinear system of ordinary differential equations, which describe transmembrane potential induced by an ionic current across a cell membrane. Depending on the chosen $\Delta t$, to solve the linear system, demiFHN is based on one of the Euler methods, while demiTP06 is based on the forward Euler method or rushlarsen method. Another important difference: demiFHN works with only one state variable, while demiTP06 works with 18 different state variables [cen25a].

To avoid redundancy and without loss of generality, I will discuss the work done on these
<figure><img src="image_24.png" /><figcaption>Figure 7: Steps performed by the models in DEMI</figcaption></figure>

two models only referring to `demiFHN`, due to its higher simplicity; the reasoning made for it can be straightforwardly applied to `demiTP06`. Moreover, for the explanation I will show the technical implementation details of only one numerical method: the forward Euler method.

### 3.1.1 Involved structures and registration

As we anticipated above, FHN computes a quantity *v* (transmembrane potential), using a support variable *w*, for a cell. Then, these variables can be easily extended to vectors of length *n*, and we can consider a collection of *n* cells, using two vectors *v* and *w*: each entry of them is associated with one cell. In `demiFHN`, *v* and *w* are represented using two DEMI data structures, called *demivector* and *demifinstate*. In Appendices 5 their definitions are completely reported. Taking a quick look at the definition of the structure *demivector*, it is easy to understand how this structure is complex: it is a special kind of vector capable of supporting partitioning and ghost entries. This structure is fundamental in DEMI and is used for many purposes. To fully handle it with StarPU, it would be necessary to define a new datatype for it, but for the needs of `demiFHN`, this was not strictly necessary and I limited myself to handle only two
normal vectors: demivector→x, (for v) and demifhnstate→w; since they are normal vectors, their registration is naturally supported in StarPU. Listing 7 shows their registrations, which are completely analogous to what we show in 1.

<figure><img src="image_25.png" /><figcaption>Listing 7: registration of v and w on StarPU</figcaption></figure>

### 3.1.2 Partitioning of handles, Euler method codelet and insertion of tasks

Iterative methods are used to solve the system of equations associated with a time moment. The iterative procedure implies time dependencies: each iteration depends on the solution computed for the previous timestep; although, this does not completely prevent us to leverage parallelism. The transmembrane potential of a cell depends only on its previous solution, therefore many cells can be computed independently at a given timestep, and the computation can easily be "vectorised".

At this point, the partitioning comes naturally: each task can work on a different subhandle independently and safely perform the forward Euler operation for each cell of the subpartition. Listing 8 shows the codelet for the task implementing the Forward Euler method, while Listing 9 shows the partitioning and the insertion of tasks for each subhandle. For conciseness, I reported only the version with the dynamic definition for the number of partitions and the asynchronous approach, but during the development process the static definition and the synchronous approach were addressed and validated as well.

<figure><img src="image_26.png" /><figcaption>Listing 8: Codelet forward euler method</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Listing 9: Partitioning and tasks submission</figcaption></figure>

### 3.1.3 DAG for FHN

Figure 8 illustrates an explicative DAG for the StarPU version of FHN, considering a setup with 4 cells and 4 partitions. In order to focus only on the behaviour of tasks and their dependencies, the DAG is limited only to two timesteps. Tasks working on the same data partition, at different moments, are linked by data dependencies, while tasks on different data partitions can be scheduled at the same time.

### 3.1.4 ViTE traces FHN model

The StarPU version is able to replicate exactly the behaviour of the original DEMI version and the results are coherent. Once I had ensured the correctness of the StarPU version for the FHN model, I performed some simple performance tests using ViTE 2.4.2.

Figure 9, shows the traces generated by the same setup shown in Figure 8: 4 cells and 4 partitions (one cell per partition), executed on node sirocco17 1. As we can easily observe, the traces are completely red, which means that for the most part of the life of the program, tasks are not executed, and StarPU is not actually beneficial. In fact, as we mentioned in 2.4.2, the
<figure><img src="image_28.png" /><figcaption>Figure 8: DAG representing two iterations for the task implementation of the forward Euler method. Among tasks working on different partitions, there are no data dependencies, while there are between tasks working on the same data partition at different iterations. Image generated using ViTE 2.4.2</figcaption></figure>

amount of work executed by tasks must be relevant, and clearly this is not the case for the
current setup. The small dimension of the problem results in a light workload assigned to
tasks and their durations are shorter than the minimum threshold, 100 µs.

<figure><img src="image_29.png" /><figcaption>Figure 9: FHN computation for 4 cells and 4 partitions (one cell per partition), executed on node sirocco17 1. Traces are almost completely red: the current amount of work assigned to tasks is irrelevant and task durations are below the minimum threshold.</figcaption></figure>

To actually observe benefits deriving from StarPU, it is necessary to feed tasks with ap-
propriate amounts of work; this means increasing the size of the problem and the number of
considered cells. Figure 10, shows the traces generated from a setup much larger than the
previous one: FHN model for 360 milions of cells, subdivided into 36 partitions (10 millions
cells per tasks), and solved on node Bora1 5. Again, a large part of the screen appears red, but
a considerable part of green is now visible. For this case, StarPU is actually inactive during
the assembly of structures: this part was not addressed with tasks, and therefore we should
not expect anything different for it; but, once all structures are ready, and the program has en-
tered in the numerical solution part, engaged with tasks, the life time of the program is mostly
occupied by actual tasks execution.

Overall, this result is positive: with a high number of cells that we can consider realistic for a real-world case, it is possible to implement a correct and worth task-based version of the numerical methods used to compute the transmembrane potential.

**3.1.5 Time performance analysis FHN model**

To better underline how the granularity of workload is relevant, some successive experiments
were conducted, and the time executions were registered and compared. In these experiments,
two execution steps were considered:

* A setup phase performed by the solver, called "Evaluation of the source term"

* Forward Euler method on the whole set of cells.

These two metrics were collected for different problem sizes, using both the serial version of DEMI and the StarPU/DEMI version; the node Sirocco15 1 was used for the experiments and all the results collected are reported in Table 2. The StarPU/DEMI version implements an asynchronous and aggressive partitioning: if the number *n* of cells is lower than the available *m* cores on the node, *n* partitions are created, otherwise *m*. This policy is kind of trivial and does not consider the amount of work that needs to be addressed: as a consequence, for small size problems, StarPU is not beneficial and the overhead introduced for the management of tasks is much larger than the benefits derived by the introduced parallelism.

Figure 11 provide a visual representation of the results of the experiments. It easy to ob-
serve how for the smallest problem sizes, due to a bad partitioning, the time executions for the
Forward Euler method are clearly greater for the StarPU/MPI version. Using this aggressive
partitioning, the StarPU/DEMI version becomes interesting only for a problem size of 2²⁰ cells,
with 32768 cells for each of the 32 available cores. After that point, the StarPU/DEMI version
is able to perform better than the serial version of DEMI, and it provides a discrete speedup
(around 50x for 2³⁰ cells). The graph also reports the times needed for the evaluation of the
source term: currently, this part is not addressed with tasks, but since it seems to grow linearly
as the problem size increases, StarPU could be used to improve the time performance for this
step.

Finally, Figure 12 shows how the speedup change with respect to the problem size: be-
fore the "crossover point" is negative and the serial version is still to be preferred, after the
StarPU/DEMI version starts to become actually better. Probably, it would be possible, design-
ing a suitable partitioning policy, to achieve better performance even for smaller problem sizes,
but it was not the goal of this analysis.
<table><caption>Table 2: Execution times for the FHN model for different problem sizes, using the serial version of DEMI and the StarPU version.</caption><thead><tr><th>n cells</th><th>Evaluation source term (s)</th><th>Forward Serial (s)</th><th>Forward StarPU/DEMI (s)</th></tr></thead><tbody><tr><td>2<sup>0</sup></td><td>4.47e-06</td><td>3.31e-07</td><td>5.65e-03</td></tr><tr><td>2<sup>1</sup></td><td>5.53e-06</td><td>3.67e-07</td><td>4.26e-03</td></tr><tr><td>2<sup>2</sup></td><td>4.47e-06</td><td>3.85e-07</td><td>1.65e-03</td></tr><tr><td>2<sup>3</sup></td><td>4.37e-06</td><td>3.78e-07</td><td>1.23e-02</td></tr><tr><td>2<sup>4</sup></td><td>3.96e-06</td><td>5.64e-07</td><td>1.65e-02</td></tr><tr><td>2<sup>5</sup></td><td>5.49e-06</td><td>9.63e-07</td><td>1.66e-02</td></tr><tr><td>2<sup>6</sup></td><td>4.65e-06</td><td>1.34e-06</td><td>2.16e-02</td></tr><tr><td>2<sup>7</sup></td><td>8.71e-06</td><td>1.99e-06</td><td>1.44e-02</td></tr><tr><td>2<sup>8</sup></td><td>4.77e-06</td><td>3.68e-06</td><td>1.24e-02</td></tr><tr><td>2<sup>9</sup></td><td>1.19e-05</td><td>7.38e-06</td><td>1.44e-02</td></tr><tr><td>2<sup>10</sup></td><td>2.66e-05</td><td>1.43e-05</td><td>2.14e-02</td></tr><tr><td>2<sup>11</sup></td><td>1.90e-05</td><td>2.84e-05</td><td>2.44e-02</td></tr><tr><td>2<sup>12</sup></td><td>3.74e-05</td><td>5.53e-05</td><td>3.45e-02</td></tr><tr><td>2<sup>13</sup></td><td>6.61e-05</td><td>1.12e-04</td><td>1.61e-02</td></tr><tr><td>2<sup>14</sup></td><td>1.32e-04</td><td>2.23e-04</td><td>2.76e-02</td></tr><tr><td>2<sup>15</sup></td><td>2.51e-04</td><td>4.66e-04</td><td>1.32e-02</td></tr><tr><td>2<sup>16</sup></td><td>5.31e-04</td><td>8.79e-04</td><td>1.83e-02</td></tr><tr><td>2<sup>17</sup></td><td>9.64e-04</td><td>1.75e-03</td><td>3.09e-02</td></tr><tr><td>2<sup>18</sup></td><td>1.96e-03</td><td>3.52e-03</td><td>1.87e-02</td></tr><tr><td>2<sup>19</sup></td><td>2.95e-03</td><td>7.10e-03</td><td>1.79e-02</td></tr><tr><td>2<sup>20</sup></td><td>4.95e-03</td><td>2.71e-02</td><td>2.34e-02</td></tr><tr><td>2<sup>21</sup></td><td>9.91e-03</td><td>5.42e-02</td><td>4.90e-03</td></tr><tr><td>2<sup>22</sup></td><td>3.17e-02</td><td>1.13e-01</td><td>2.63e-03</td></tr><tr><td>2<sup>23</sup></td><td>6.15e-02</td><td>2.05e-01</td><td>2.00e-02</td></tr><tr><td>2<sup>24</sup></td><td>1.35e-01</td><td>4.38e-01</td><td>9.92e-03</td></tr><tr><td>2<sup>25</sup></td><td>2.87e-01</td><td>9.13e-01</td><td>3.18e-02</td></tr><tr><td>2<sup>27</sup></td><td>1.10e+00</td><td>3.83e+00</td><td>8.62e-02</td></tr><tr><td>2<sup>28</sup></td><td>2.23e+00</td><td>7.25e+00</td><td>1.65e-01</td></tr><tr><td>2<sup>29</sup></td><td>4.49e+00</td><td>1.89e+01</td><td>3.44e-01</td></tr><tr><td>2<sup>30</sup></td><td>9.01e+00</td><td>3.12e+01</td><td>6.67e-01</td></tr><tr><td>2<sup>31</sup></td><td>1.86e+01</td><td>6.53e+01</td><td>1.31e+00</td></tr></tbody></table>
<figure><img src="image_30.png" /><figcaption>Figure 10: FHN computation for 360 milions of cells, subdivided in 36 partitions (10 millions of cells per tasks), executed on node Bora1 1. The prevalence of green segments in correspondence of the numerical solving phase indicates that the tasks are fed with enough work to be worthy</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Figure 11: Time executions for the FHN model. The StarPU version works on the Forward Euler method using partitions. The evaluation of the source term is serial and it is the same for both versions. The crossover point is placed where the StarPU version starts to have better time performance than the serial version of DEMI for the Forward Euler Method. The graph has two different scales on the vertical axis: the green line for partitions is referred with respect to the scale on the right side, while the time executions with respect to the scale placed on the left side.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 12: Gained speedup using the StarPU version for the Forward Euler method. Before the crosspoint, at $2^{20}$ cells, the serial version performs better: tasks are not fed with enough amount of work and the scheduling of task is more relevant than their execution. Beyond the crosspoint, the StarPU version performs better than the serial version and the speedup is positive.</figcaption></figure>
## 3.2 DIFFUSION model

demiDIFFUSION model solves time-dependent diffusion equations in a two- or three-dimensional case. The time-dependent diffusion equation is a linear, elliptic partial differential equation, which describes the time evolution of some quantity or concentration, which is subject to diffusion, meaning that it flows from regions of high concentration to regions of lower concentration [cen25a]. The linear system arising for each time step is solved with the conjugate gradient (CG) method, which is an iterative numerical method used to solve linear systems of equations of type $Ax = b$, where $x$ is the vector of unknowns; a simplified pseudocode for it is reported in Algorithm 1.

My principal goal for the demiDIFFUSION model was to redesign the CG using tasks, implementing a single codebase on top of StarPU capable of handling both single-node and multi-node computations. The demiDIFFUSION model was more challenging than the previous ones; indeed, in addition to the single-node case, I also addressed the StarPU-MPI approach. Introducing distributed computation can expand the possibility of software, but this comes with further requirements: correct partitioning and distribution of data structures. A relevant part of the work on the demiDIFFUSION model was dedicated to the data handling of the sparse structures involved in the CG method, and new StarPU datatypes were defined. Actually, even if this would have been necessary only once distributed computation with MPI had been addressed, for the reasons that I will explain later, I decided to deal with that firstly; I took this decision to become more familiar with the structures used in DEMI from the beginning and to be able to proceed more smoothly later on.

In Section 3.2.1, I will talk about the process of defining a new type in StarPU and explain why it was necessary in more detail. After that, in Section 3.2.2, I will present the single node version for the demiDIFFUSION model; this intermediate step was useful to outline the workflow of tasks required by CG. Once done that, I moved to the multiple node case: in Section 3.2.3, I will explain how I partitioned and distributed data among MPI processes; in Section 3.2.4, I will present the final result of the StarPU-MPI approach, while in 3.2.5 I will discuss some related ViTE traces. Finally in 3.2.6, I will compare the time performance of all the different implementations of the model.
<caption>**Algorithm 1:** Conjugate Gradient Method, based on residual norm criteria</caption>

Input: Matrix $A$, vector $b$, initial guess $x_0$, max_iter, $\varepsilon$
Output: Approximate solution $x$
1
2 $r \leftarrow b - Ax_0$ // compute initial residual
3 $p \leftarrow r$ // compute initial direction
4
5 $r_{\text{norm}} \leftarrow \sqrt{\langle r, r \rangle}$ // compute initial residual norm
6
7 /* Initial convergence test */
8 $\lfloor$ return $x_0$ // $x_0$ is already a solution
9
10 for $k \leftarrow 0$ to $\text{max_iter}$ do
11 $t \leftarrow \delta t$
12 $t \leftarrow \gamma A p + t$
13 $\alpha \leftarrow (r, r)/(p, t)$
14 $r_{\text{old}} \leftarrow r$
15 $r \leftarrow r - \alpha t$ // update residual
16 $r_{\text{norm}} \leftarrow \sqrt{\langle r, r \rangle}$ // compute residual norm
17 $x \leftarrow x + \alpha p$ // update solution
18
19 /* convergence test */
20 $\lfloor$ return $x$ // $x$ is a solution
21
22 $\beta \leftarrow (r, r)/(r_{\text{old}}, r_{\text{old}})$
23 $p \leftarrow r + \beta p$
24
25 $\lfloor$ if $k == \text{max_iter}$ and $r_{\text{norm}} \ge \varepsilon$ then
26 $\lfloor$ return NOT CONVERGED
27 else
28 $\lfloor$ return $x$

3.2.1 Datatype definitions

The CG method works with two types of algebraic structures: vectors and matrices. In DEMI, there exist particular data structures for them: *demivector* and *demihostspblassymv*, the definition of which is reported in the Appendices 5. They are quite complex structures, capable of supporting sparse format and distributed computation. Since they are structures designed specifically for DEMI, they are not naturally supported in StarPU and this implies that to completely handle them with StarPU it is necessary to define a new datatype for them: data types are needed by StarPU to allocate, replicate and move pieces of data of that type. Defining a new datatype in StarPU is not mandatory, but it can be convenient, especially in a distributed context, and I did it for the structures *demivector* and *demihostspblassymv*.

The definition of a new datatype in StarPU is a quite complex process and requires around
a few hundred lines of code; in the following I limit myself to underlining the most interesting
parts, reporting some extracts of the definition for the structure demivector.

**Datatype operations** A new datatype has to specify some functions that will be used by StarPU for its work. The functions are listed in a special StarPU structure, called `starpu_data_interface_ops`; in Listing 10, the operations structure used for demivector is reported. For each field of the structure, I defined a specific function; in the following paragraph, I will discuss the ones that are the most interesting for the purpose of the discussion: `starpuDemivectorallocate_data_on_node` 3.2.1, `starpuDemivector_get_size` 3.2.1 and finally `starpuDemivector_pack_data` 3.2.1. Again, to avoid tediousness, I do not report completely their definition, but only the parts that I really consider useful for a better understanding of the successive sections, especially for those ones where the multinode framework is discussed.

<figure><img src="image_33.png" /><figcaption>Listing 10: Operations for demivector datatype</figcaption></figure>

**Allocate data on Node** Each piece of data, to be correctly registered, needs to be allocated on a physical memory region, a StarPU node. The function `starpuisu allocate data on node` specifies how to do that: with respect to the structure attributes and using C memory management, the strictly necessary amount of memory is allocated on the node. Listing 11, shows a simplified example of the function used for the demivector interface.
<figure><img src="image_34.png" /><figcaption>Listing 11: Allocation on node</figcaption></figure>

**Get size** To perform data MPI communications, StarPU needs to know the size of the structure. For this reason, it is convenient to have a support function for the computation of the size. In Listing 12, the example for `demivector`. Note that there is also the size for some attributes not present in 11.

<figure><img src="image_35.png" /><figcaption>Listing 12: Get size of the handle</figcaption></figure>

**Pack data** Before a data transfer from a StarPU node to another, it is necessary "to pack" data and prepare them: this means copy non-contiguous pieces of data into contiguous buffer. Therefore, the content of the subhandle owned by the node is a copy to a pointer. In Listing 13, the pack procedure for our current example.
<figure><img src="image_36.png" /><figcaption>Listing 13: Pack data</figcaption></figure>

### 3.2.2 CG on single node

Again, one fundamental step in dealing with StarPU is the design of code using tasks. Therefore, as I did for the transmembrane models 3.1, firstly I analysed the CG operations and designed a task workflow for them. Then, I implemented the different tasks and checked the correctness for single-node computations. Figure 13, generated with ViTE, shows the workflow of tasks for a single iteration of the internal loop of the CG method in the current implementation; the colours distinguish the type of operation, while the text inside the node specifies the structures involved. The reader could notice that not all the expected data dependencies are shown; this is true, and it is due to the fact that some `MPI_All_Reduce` routines are present in my code and are underlined by the black nodes. These reductions ensure correct results, and they are necessary for the computation of variables after a dot product. Aware of this, I decided to not explicitly specify some dependencies involving variables because it was not strictly necessary, but probably there will be inserted in future works. Of course, these MPI reductions limit the parallelism, but again, we can relies on handles partitioning, as I show in 9. Unfortunately, due to time constraints, I did not implement the partition for the `demiDIFFUSION` model, preferring to proceed with the multinode case.
<figure><img src="image_37.png" /><figcaption>Figure 13: Graph with the tasks implemented for a single iteration of the main loop of CG; the colour node distinguish the operation type, while the label specifies the structures involved. N.B. this is the first iteration of the internal loop, and therefore some pre-loop operations do not appear</figcaption></figure>

It is also interesting to look at two consecutive iterations of CG; I report them in Figure 14, where tasks for the first iteration are indicated by orange, while purple is used for the second iteration. With respect to the Forward Euler Tasks, discussed in 3.1.3, among consecutive iterations, data dependencies are less strict, and it is possible to schedule tasks of different iterations and working on the same partition, at the same time.
<figure><img src="image_38.png" /><figcaption>Figure 14: In orange tasks of the first iteration, in purple tasks of the second iteration and in black synchronization routines</figcaption></figure>

The correctness of the single-node StarPU version was checked and the results obtained were compared with those of the original serial version of DEMI. Meshes with different sizes and different dimensions were tested, and the StarPU version was able to perfectly replicate the desired results. Finally, it is important to remember that all the test cases addressed were "ideal", i.e. square or cubic meshes. Figure 15 shows an experiment in which the two versions were compared in a 2d case, while Figure 16 shows another experiment with the same setup, but in 3d; both show the final instant of the solution.
<figure><img src="image_39.png" /><figcaption>(a) Final instant of the solution computed with DEMI original version</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(b) Final instant of the solution computed with StarPU-DEMI single node version</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(a) Final instant of the solution computed with DEMI original version</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(b) Final instant of the solution computed with StarPU-DEMI single node version</figcaption></figure>

<caption>Figure 16: Comparison experiment results of DEMI and StarPU-DEMI. Diffusion problem solved for an ideal 3 dimensional case</caption>

### 3.2.3 Data partitioning and distribution

In this part, I will explain how I partitioned and distributed the data among MPI processes for the multinode approach.

**Partitioning** The term partitioning is used differently here: in 9 we discuss the StarPU partitioning of handles inside the same StarPU node, while here we refer to the partitioning of data structures among MPI processes, on different nodes. I chose a simple partitioning: each process manages a certain number of rows of the matrix A, all the related column entries and the corresponding entries in the vectors. Figure 17 shows the split considering a simple case where only two MPI processes are involved.
<figure><img src="image_43.png" /><figcaption>Figure 17: Partitioning of the structures involved in CG for a 4x4 linear system among two MPI processes</figcaption></figure>

**Distribution** At the time of my work, the MPI version of DEMI was still in development and for this reason I built the StarPU-MPI version starting from the serial version of the demiDIF-FUSION code. Therefore, while the current DEMI exploits METIS in combination with MPI to distribute the structures, my version adopts a more simple and manual approach: the rank 0 is in charge of the assembly and distribution of all the algebraic structures; currently this is the main bottleneck present in the StarPU-MPI version and it will be addressed in future works. I decided to opt for this approach for the following reason: the assembly of structures is based on the nodes of the mesh; therefore, to permit a process to assemble its part of structures, it would be necessary to distribute correctly the sub-meshes. This can be very complicated: meshes can have particular geometries and ghost elements have to be taken into account; for these reasons, dedicated tools such as METIS [ sof25] (a software package for graph partitioning) are used. Using METIS, it is possible to distribute the nodes of the meshes and have independent assembly procedures, reducing the synchronisation routines to only those necessary for the solution of CG. My current implementation does not rely on any tool for the distribution: Rank 0, at each time step, with respect to the current solution and the nodes of the mesh, assembles and distributes the structures serially. As I said, this surely represents a limitation, but since the MPI-METIS implementation was already in development, I preferred to focus my efforts on the numerical part of the application, while designing tasks that could be able to work even in combination with other distribution strategies: their current implementation logic is orthogonal at the distribution procedure, and this will allow their functionality also once the assemble part will be addressed more efficiently.

### 3.2.4 CG multiple nodes

To illustrate the behaviour of the CG procedure in the multi-node approach, we will follow an example with two MPI ranks working on the setup reported in 17. Figure 18, shows a single iteration of the method. Expect for the `MPI_AllReduce` routine and the shared task, the two MPI processes can proceed independently and work safely on their local data. The shared task deserves particular attention: it is not simply a task, but is a *mpi task*. This means that each
process submits the task locally, but part of needed data is located in another physical node; StarPU, looking at the specified data dependencies, handles the MPI data communications automatically and transparently.

Figure 19, shows the two different points of view of the two processes when they submit the task discussed and the operations they have to perform: they have to update their entries for the vector *t*, but to do so, they also need some ghost entries of *p*. For simplicity, only in the figure we assume that *t* is the null vector and *γ* is equal to zero.

<figure><img src="image_44.png" /><figcaption>Figure 18: First iteration of the internal loop with 2 MPI processes.</figcaption></figure>
<figure><img src="image_45.png" /><figcaption>Figure 19: Each process needs to update its own entries of vector *t* and ghost entries of vector *p* are needed. For simplicity, only in the figure, we assume that *t* is the null vector and *γ* is equal to zero</figcaption></figure>

In Listing 14, I also report the task submission code. It is interesting to observe that each
process submits the task for each sub-handle, but then the computation for a specific subpart
of *t* is performed by the MPI process actually working on the physical node where that subpart
is stored, and data MPI transfers of subparts of *p* are inferred consequently.

<figure><img src="image_46.png" /><figcaption>Listing 14: Codelet mul_matrix_vector</figcaption></figure>

Finally, the results were validated using the same approach as for the single-node case. Figures 20 show the comparison for the same test cases that have already been seen; again, all the results are coherent.

3.2.5 ViTE analysis for CG multiple node

As I said, since this work represents only the first attempt of the use of DEMI on top of StarPU,
the correctness of results was the main concern. Anyway, some ViTE analyses were conducted
and I think that the generated traces can very well underline the current limitations and pro-
vide more awareness for the following work. Figure 21, shows the traces generated by the
execution of the demiDIFFUSION models, using 4 MPI processes. The white region is the
result of the overlapping of many white arrows, which are used to indicate the MPI data trans-
fers that occur between processes. As is easy to understand, they represent a relevant part of
the whole life of the program, and this naturally has an impact on performance.
<figure><img src="image_47.png" /><figcaption>(a) Final instant of the solution computed with StarPU-DEMI multiple node version, for the 2 dimensional case</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>(b) Final instant of the solution computed with StarPU-DEMI multiple node version, for the 3 dimensional case</figcaption></figure>

Figure 20: Results obtained with StarPU-DEMI multiple node version; they are identical to 15 and 16

It is also interesting to see the directions of the white arrows in more detail; Figure 22, shows a zoom on a specific area of the previous image (i.e. Figure 21). Even if there are also some intra-node transfers, the majority are represented by those involving rank 0 and another process. This reflects what we expect from the current distribution design, and it underlines its limitations.

### 3.2.6 Time performance analysis DIFFUSION model

Also, for the demiDIFFUSION model, some time statistics were collected and the different versions evaluated. Again, using the meshes chosen by the DEMI developers for the testing of the model, different problem sizes were addressed and all the collected results are reported in Tables 3 4 5 6. The meshes are simply referred by numerical identification, and their features are reported in 2. Note that the StarPU tables does not have the statistic for the solver setup: this part is still serial in any implementation and even if it is performed by all the code versions, I reported it only for the serial one, in Table 3. For the experiments I used the node Sirocco15 1.

<table><caption>Table 3: Execution times for the serial version of DEMI.</caption><thead><tr><td>Mesh</td><td>Nodes</td><td>Elements</td><td>Solver Setup [s]</td><td>Assembly [s]</td><td>DOT Product [s]</td><td>CG [s]</td></tr></thead><tbody><tr><td>Mesh3</td><td>143</td><td>645</td><td>26.701716</td><td>6.7e-2</td><td>4.10e-6</td><td>1.37e-4</td></tr><tr><td>Mesh4</td><td>676</td><td>3492</td><td>26.0744693</td><td>5.62e-2</td><td>5.96e-6</td><td>2.91e-3</td></tr><tr><td>Mesh5</td><td>4002</td><td>22505</td><td>26.303911</td><td>6.87e-2</td><td>4.61e-2</td><td>4.61e-2</td></tr><tr><td>Mesh6</td><td>27383</td><td>162665</td><td>26.267617</td><td>1.90e-1</td><td>1.36e-4</td><td>6.61e-1</td></tr><tr><td>Mesh7</td><td>200992</td><td>1228034</td><td>35.275308</td><td>1.09</td><td>9.67e-4</td><td>1.48e1</td></tr><tr><td>Mesh8</td><td>1533356</td><td>9499616</td><td>101.471332</td><td>1.08e1</td><td>5.89e-3</td><td>2.42e2</td></tr></tbody></table>

Figures 23 24 provide graphical representations of the results.

Figure 23, shows a time comparison between the different configurations for the assembly of structures used by the CG method. The StarPU versions requires always more time than the serial version to perform this part; actually, this is expected and the results are not overall so negative: currently, even the StarPU versions rely on a serial assemble procedure and
<figure><img src="image_49.png" /><figcaption>Figure 21: Generated traces by the execution of the demiDIFFUSION model using 4 MPI processes. The white region is the result of the overlapping of many MPI data communication, indicated in ViTE by white edges.</figcaption></figure>

<table><caption>Table 4: Assembly times for StarPU versions of demiDIFFUSION model</caption><thead><tr><th>Configuration</th><th>Mesh3</th><th>Mesh4</th><th>Mesh5</th><th>Mesh6</th><th>Mesh7</th><th>Mesh8</th></tr></thead><tbody><tr><td>Single node</td><td>8.10e-2</td><td>5.82e-2</td><td>1.06e-1</td><td>1.71e-1</td><td>1.12</td><td>1.08e1</td></tr><tr><td>2 MPI</td><td>7.38e-2</td><td>1.15e-1</td><td>1.64e-1</td><td>2.43e-1</td><td>1.10</td><td>1.09e1</td></tr><tr><td>4 MPI</td><td>8.32e-2</td><td>1.54e-1</td><td>2.00e-1</td><td>4.56e-1</td><td>3.18</td><td>2.19e1</td></tr></tbody></table>

<table><caption>Table 5: Dot product times for StarPU versions of demiDIFFUSION model</caption><thead><tr><th>Configuration</th><th>Mesh3</th><th>Mesh4</th><th>Mesh5</th><th>Mesh6</th><th>Mesh7</th><th>Mesh8</th></tr></thead><tbody><tr><td>Single node</td><td>8.93e-3</td><td>6.06e-5</td><td>7.56e-5</td><td>1.62e-4</td><td>7.78e-4</td><td>9.57e-3</td></tr><tr><td>2 MPI</td><td>1.57e-5</td><td>1.43e-4</td><td>1.42e-4</td><td>1.86e-4</td><td>5.13e-4</td><td>4.74e-3</td></tr><tr><td>4 MPI</td><td>3.58e-2</td><td>1.28e-2</td><td>1.27e-2</td><td>1.31e-1</td><td>2.94e-2</td><td>1.95e-2</td></tr></tbody></table>

<table><caption>Table 6: CG times for StarPU versions of demiDIFFUSION model</caption><thead><tr><th>Configuration</th><th>Mesh3</th><th>Mesh4</th><th>Mesh5</th><th>Mesh6</th><th>Mesh7</th><th>Mesh8</th></tr></thead><tbody><tr><td>Single node</td><td>8.89e-3</td><td>3.42e-2</td><td>7.37e-2</td><td>9.13e-1</td><td>2.40e1</td><td>4.49e2</td></tr><tr><td>2 MPI</td><td>6.30e-3</td><td>2.41e-2</td><td>1.04e-1</td><td>9.49e-1</td><td>1.35e1</td><td>2.42e2</td></tr><tr><td>4 MPI</td><td>1.11</td><td>5.65</td><td>1.14e1</td><td>2.30e1</td><td>5.91e1</td><td>4.92e2</td></tr></tbody></table>

they have to deal in addition with the overhead introduced by the use of StarPU and the MPI communications used for the distribution. As we can see, the overhead tends to become more negligible as the mesh size is increased. However, as we mentioned before, the assemble of structures will need to be parallelised: as the size of the problem is increased, the procedure becomes more expensive, and a serial approach will surely represent a limitation.

Figure 24 shows a time comparison between the different configurations for the CG method. As in the previous graph, the differences tend to reduce as the mesh size increases. However, for the CG procedure, we would like to clearly observe some benefits derived by the use of StarPU, but this is not the case. Actually, the algebraic operations involved in the CG method are rather simple, and a parallelisation working on the different parts of the structures be-
<figure><img src="image_50.png" /><figcaption>Figure 22: Specific section of the image 21. The majority of MPI communications are point-to-point messages between the Rank 0 and another rank.</figcaption></figure>

comes effective and worthy only after a certain threshold size. Anyway, this can be achieved
only with greater test cases, but this also implies having more demanding requirements for the
solver setup and the assembly of structures. In conclusion, we are in a sort of stalemate: at the
same time we would like to test the capabilities of the StarPU versions for bigger test cases, but
the serial setup and the serial assembly prevent us from doing so, since the experiments be-
come almost impractical. For all these reasons, I would suspend this performance evaluations
for now and take up again once tasks have been used to address larger portions of code.
<figure><img src="image_51.png" /><figcaption>Figure 23: Assembly times comparison for different configurations. This procedure is always carried out serially, but the StarPU versions are penalised by additional overheads.</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>Figure 24: CG times comparison for different configurations.</figcaption></figure>
# 4 Conclusions

Cardiac electrophysiology is the branch of medicine dedicated to electric activity in the heart. The conduction of electricity in the human heart is crucial, and it is well established that abnormalities in cardiac conduction are associated with fatal consequences: disorders of the heart's electrical synchronisation system, named arrhythmias, are responsible for half of the deaths related to cardiovascular disease. Fortunately, thanks to many research efforts, accurate mathematical modelling and computer models, cardiac electrophysiology has made great progress in the last decades; still, many questions are without answer and for some specific investigation current computer models seem to not be sufficient. Several current cardiac computer models adopt the so-called "continuous" approach: cardiac cells are grouped and the same behaviour is assumed for all the cells belonging to the same group. In many cases, this approximation offers great results: it contains the computational demands without sacrificing reliability.

However, sometimes this simplistic assumption does not offer enough precision and consequently current computer models lack in accuracy, and they cannot provide an extensive replication of some events determined by reduced electrical coupling between cells and related to cardiac disease; therefore, other cardiac models, named cell-based cardiac models, were proposed to fill this gap. These models have shown their capabilities, and they have been appreciated in situations where continuous classical models are not sufficient. Naturally, their higher accuracy comes with an increase in terms of computational cost, which can become extremely high considering the large number of cells present in the human heart.

Special infrastructures, named supercomputers or High Performance Platforms (HPC), can provide the necessary resources needed by cell-based models. Anyway, HPC resources management can be challenging: competencies, portability and efficiency are required. In fact, HPC platforms are complex and variegate systems, composed of different type of hardware; there is not uniformity among different platforms and they evolve constantly; finally, efficiency is obtainable only with intelligent management of the resources, but again, this is not trivial. For all these motivations, it can be convenient to rely on runtime systems: special software in charge of the execution of scientific applications. They represent a valid help for the programmer, who can exploit a defined general interface avoiding to take direct decisions about the hardware management.

Runtime systems are largely used for scientific applications and cardiac cell-based models can benefit from them as well. This project aimed to investigate this possibility and this document describes the first attempt ever to combine a cardiac cell-based model, DEMI, with a runtime system, StarPU. DEMI is composed of many subparts; part of them were completely handle, while others were only partially, and in the next paragraphs I will propose how the work could continue in the short, medium, and long term.

**Short term** Surely the most immediate improvement would be the creation of handles for the variables that are the results of dot products tasks in the CG procedure. As I mentioned in 3.2.2, even if these data dependency exists, MPI routines ensure the correct behaviour even if the former are not detected. Anyway, even if managing this aspect within a task is not strictly necessary to ensure correctness, it would be better to handle it with StarPU for greater flexibility in scheduling and in overlapping communications and computations. Indeed, an implicit MPI reduction using StarPU data dependencies would be able to start as soon as some of the results are available from the reduction step, instead of waiting on a barrier.

Moreover, for the short term, I will also include the porting of the *demiDOMAIN model*. This model is just the combination of a transmembrane model with the demiDIFFUSION model. Once the distribution routine has been adapted, it will be able to run and execute.
Without significant efforts, another model DEMI model could be ported with StarPU.

**Medium term** In the medium term, I will suggest a redesign of the distribution mechanism using METIS. As I said in 3.2.3, at each iteration, CG needs the assembly of the involved structures, and this is done with respect to the nodes of the mesh; more processes could perform the assembly for their own local structures, but only if the mesh is subdivided into submeshes and only if then they are distributed correctly; this can be very challenging. Fortunately, there exist means like METIS for this purpose, and the most updated version of DEMI already uses it. At the time of the beginning of my work, the combination for DEMI of MPI and METIS was still in development and, for this reason, a more naive approach, described in 3.2.3, was adopted. The current procedure limits the parallelism since it needs the serial assembly by one process and then the serial distribution of structures. Naturally, as shown in Figures 21 22, this represents the main bottleneck of the current implementation. At some point, it will be necessary to combine the new version of DEMI with the current StarPU-DEMI version. Aware of this issue, I designed them to be orthogonal to any distribution procedure and they will be able to work also in combination with METIS. Additionally, for the medium term, I will also suggest the implementation of task also for the wrapping of the METIS functions: in this way, it would also be possible to deal with this onerous part with StarPU. Unfortunately, distributing structures with StarPU also requires their registration: a new datatype for the structure used in DEMI for meshes will have to be defined.

**Long term** Only in the long term, DEMI will be built completely on top of StarPU: before that, it will be necessary to tackle the remaining submodels, which are even more complicated than those addressed in this document. However, the remaining models are still built on top of those addressed, and I am confident that the work done will be useful. Finally, at some point, GPUs will be part of the implementation: real test cases need sufficient power, and GPUs can provide it. Moreover, a relevant part of execution in DEMI could exploit massive parallelism and therefore GPUs will be very suitable means.
References

[TP06] Kirsten HWJ Ten Tusscher and Alexander V Panfilov. "Alternans and spiral breakup in a human ventricular tissue model". In: *American Journal of Physiology-Heart and Circulatory Physiology* 291.3 (2006), H1088-H1100.

[Sun+07] Joakim Sundnes et al. *Computing the electrical activity in the heart*. Vol. 1. Springer Science & Business Media, 2007.

[Agu+17] Emmanuel Agullo et al. "Achieving High Performance on Supercomputers with a Sequential Task-based Programming Model". In: *IEEE Transactions on Parallel and Distributed Systems* (2017), pp. 1-1. DOI: 10.1109/TPDS.2017.2766064.

[Tve+17] Aslak Tveito et al. "A cell-based framework for numerical modeling of electrical conduction in cardiac tissue". In: *Frontiers in Physics* 5 (2017), p. 48.

[Jæg+19] Karoline Horgmo Jæger et al. "Properties of cardiac conduction in a cell-based computational model". In: *PLoS computational biology* 15.5 (2019), e1007042.

[Aum20] Olivier Aumage. "Instruments of Productivity for High Performance Computing". Habilitation à diriger des recherches. Université de Bordeaux (UB), France, Dec. 2020. URL: https://inria.hal.science/tel-03105625.

[Rei+23] Joyce Reimer et al. "Physiological accuracy in simulating refractory cardiac tissue: the volume-averaged bidomain model vs. the cell-based EMI model". In: *arXiv preprint arXiv:2304.04115* (2023).

[CT24] Xing Cai and James Trotter. *Enabling Cell-Based Simulations of Cardiac Electrophysiology with High-Performance Computing*. Baltimore, MD, Mar. 2024. URL: https://meetings.siam.org/sess/dsp_talk.cfm?p=132650.

[Jæg+24] Karoline Horgmo Jæger et al. "Evaluating computational efforts and physiological resolution of mathematical models of cardiac tissue". In: *Scientific Reports* 14.1 (2024), p. 16954.

[Boo25] JUPITER Booster. *JUPITER Booster, website*. 2025. URL: https://www.fz-juelich.de/en/news/archive/press-release/2025/jupiter-supercomputer-propels-european-computing-power.

[Bor25] INRIA Bordeaux. *ViTE website*. 2025. URL: https://solverstack.gitlabpages.inria.fr/vite/index.html.

[cen25a] Simula lab center. *DEMI software, Github repository*. 2025. URL: https://github.com/simulahpc/demi.

[cen25b] Simula lab center. *Simula lab center, website*. 2025. URL: https://www.simula.no/.

[Fug25] Supercomputer Fugaku. *Supercomputer Fugaku, website*. 2025. URL: https://www.f-ccs.riken.jp/en/fugaku/about/.

[HPC25] HPC6. *HPC6, website*. 2025. URL: https://www.eni.com/en-IT/actions/energy-transition-technologies/supercomputing-artificial-intelligence/supercomputer.html.

[INR25] Team STORM INRIA. *Team STORM, website*. 2025. URL: https://team.inria.fr/storm/.
[LaB25] laboratory LaBri. LaBri, website. 2025. URL: https://sin.u-bordeaux.fr/les-unites-de-recherche/labri.

[Pro25] TOP 500 Project. TOP 500 Project, website. 2025. URL: https://top500.org/site/51016/.

[sof25] METIS software. METIS software, Github repository. 2025. URL: https://github.com/KarypisLab/METIS.

[Bora] Centre Inria de Bordeaux. Centre Inria de Bordeaux website. URL: https://www.inria.fr/fr/centre-inria-universite-bordeaux.

[Borb] Centre Inria de Bordeaux. Software using StarPU. URL: https://starpu.gitlabpages.inria.fr/software.html.

[Borc] Centre Inria de Bordeaux. StarPU publications. URL: https://starpu.gitlabpages.inria.fr/publications.html#PublicationsOnApplications.

[Bord] INRIA Bordeaux. StarPU repository. URL: https://gitlab.inria.fr/starpu/starpu.

[Bore] INRIA Bordeaux. StarPU website. URL: https://starpu.gitlabpages.inria.fr/.

[Borf] Université de Bordeaux. Université de Bordeaux website. URL: https://www.u-bordeaux.fr/.

[pro] MICROCARD project. MICROCARD project. URL: https://microcard.eu/index-en.html.

[web] PlaFRIM website. PlaFRIM platform. URL: https://www.plafrim.fr/.
# 5 Appendices

## 5.1 DEMI structures

<caption>Listing 15: struct demivector</caption>

```c
/*
 * demivector represents a vector stored as a contiguous array of
 * elements in full or packed storage format. The vector may also be
 * partitioned and distributed across multiple processes.
 */
struct demivector
{
    //nparts is the total number of partitions of the vector
    int nparts;

    //parttag is a tag that may be used to identify the subvector.
    int parttag;

    //nprocs is the number of processes for distributed vector
    int nprocs;

    //npparts is the number of subvectors assigned to the current process.
    int npparts;

    //ownerrank is the rank of the process that owns the subvector
    int ownerrank;

    //ownerpart is the partition number of the subvector
    int ownerpart;

    //size is the number of vector elements.
    demi_idx_t size;

    //pointertype to determine the mean
    enum demipointertype pointertype;

    //x is an array containing the (nonzero) vector values
    double *x;

    //d.x is a pointer to a device-side array containing the
    // (nonzero) vector values, see x.
    double *d.x;

    //num.nonzeros is the number of explicitly stored vector entries
    demi_idx_t num_nonzeros;

    //idxbase
    int idxbase;

    //idx is an array of length num_nonzeros, containing the
    //offset of each nonzero vector entry.
    demi_idx_t *idx;

    //d_idx is a pointer to a device-side array of length num_nonzeros, //containing the offset of each nonzero vector entry
    demi_idx_t *d_idx;

    //num_ghost_nonzeros is the number of additional ghost entries that are //stored in the vector but are ignored during
    //certain operations
    demi_idx_t num_ghost_nonzeros;
};
```

<caption>Listing 16: struct demifhn</caption>

```c
/*
 * demifhnstate is a data structure representing the state (e.g.,
 * initial condition or solution) of an ionic cell model at a given
 * time point
 */
struct demifhnstate
{
    double t;
    demi_idx_t ncells;
    double *w;
};
```

<caption>Listing 17: struct demihostspblassymv</caption>

```c
/*
 * demihostspblassymv is a data structure used as a workspace for
 * sparse symmetric matrix-vector multiply operations.
 *
 * More specifically, it is used to expand and store the nonzeros of a
 * sparse matrix from upper triangular storage to a full storage
```
7 * format, which includes both upper and lower triangular parts.
8 *
9 * The symmetric sparse matrix is stored in a packed storage format,
10 * which consists only of the upper triangle of the symmetric sparse
11 * matrix in compressed row format.
12 *
13 * In some cases, it is desirable to store the matrix in a "full"
14 * storage format, where both upper and lower triangular elements are
15 * present in the same compressed row storage. In particular, sparse
16 * matrix-vector multiplication can then be performed rowwise in
17 * parallel without data conflicts between threads.
18 */
19 struct demihostspblassymv
20 {
21 /*
22 * full (lower and upper triangular) matrix storage
23 */
24
25 /*
26 * fnpnzs is the number of nonzeros in the submatrix when it is
27 * stored in full storage format.
28 */
29 demi_idx_t fnpnzs;
30
31 /*
32 * fnzidx is an array of length fnpnzs containing, for every
33 * nonzero in the full storage format, the location of the
34 * corresponding nonzero in the symmetric, upper triangular
35 * storage of the matrix.
36 */
37 demi_idx_t * fnzidx;
38
39 /*
40 * frowptr is an array of length nprows+1 containing the
41 * offsets to the first nonzero of each row of the matrix when it
42 * is stored in full storage format.
43 */
44 demi_idx_t * frowptr;
45
46 /*
47 * fcolidx is an array of length fnpnzs containing the column
48 * number of the column j in the pair (i,j) for every nonzero
49 * of the submatrix when it is stored in full storage format.
50 */
51 demi_idx_t * fcolidx;
52
53 /*
54 * fa is an array of length fnpnzs, which contains the nonzero
55 * matrix values of the matrix in full storage format.
56 */
57 double * fa;
58 };