GPU-accelerated matrix-free solvers for the efficient solution of cardiac electrophysiology in lifex

Computational cardiology is becoming an increasingly popular field to unveil novel mechanisms of the heart and to support the clinical decision making process. In this framework, accurate and efficient numerical simulations exploiting biophysically detailed mathematical models play a major role. The entire cardiac function is driven by electrophysiology, which is characterized by a sharp and fast wavefront propagation. This multiscale problem can be modeled via the monodomain equation (a nonlinear reaction–diffusion equation) coupled with a system of nonlinear ordinary differential equations for both the ionic species and the gating variables. In this thesis, we choose a space discretization based on the high–order Spectral Element Method (SEM) as it allows high accuracy for this specific mathematical problem; the solution is therefore expressed as a linear combination of elementwise local polynomials weighted by their corresponding degrees of freedom (DoFs). We introduce two different time discretizations for the monodomain equation, namely the 2nd-order Backward Differentiation Formula (BDF2) and the 2nd-order Runge-Kutta scheme (RK2). The former choice leads to a semiimplicit time scheme yielding at each time-step a linear system solved by the conjugate gradient method (CG). On the contrary, the latter allows to compute the solution at each time-step with a much faster DoF-by-DoF update, provided that Legendre-Gauss-Lobatto (LGL) quadrature formulas are used to approximate integrals yielding SEM with Numerical Integration (SEM-NI). The application of SEM to ionic model allows, also in this case, to update the solution DoF-by-DoF as the equations are independent of each other, once they have been discretized in space. Additionally, we discretize the ionic model in time thanks to the BDF2 method. Several CPU-based solvers for cardiac electrophysiology are currently available. The exploitation of graphics processing units (GPUs) has proven to be an effective alternative, even though few matrix-based GPU solvers for the cardiac function have been proposed in the literature. We implement novel scalable and efficient, GPU-based and CPU-based, matrix-free solvers for cardiac electrophysiology, by considering semi-implicit and explicit time discretizations. We perform a thorough comparison among the different solvers under p-refinement, i.e. by increasing the polynomial degree of the basis functions. The BDF2-BDF2 GPU-based solver has been proven to be more efficient than the BDF2-BDF2 CPU-solver for serial execution, while it achieves better performances with high polynomial degrees for parallel execution. On the other hand, the RK2-BDF2 CPU-based solver, which does not require the assembly of a linear system, is slower than the BDF2-BDF2 CPU-based solver due to the repeated evaluation of different quantities for each time-step.

La cardiologia computazionale sta diventando un ambito di ricerca molto popolare perchè consente di svelare nuovi meccanismi del cuore e di aiutare i cardiologi nelle loro decisioni quotidiane. Per questa ragione, simulazioni numeriche accurate ed efficienti basate su modelli matematici che descrivono le dinamiche biofisiche rivestono un ruolo di primo piano. L’intera attività cardiaca è guidata dall’elettrofisiologia che è caratterizzata da un fronte d’onda estremamente veloce e con una forte pendenza. Questo fenomeno multiscala può essere modellizzato matematicamente tramite l’equazione del monodominio (un’equazione nonlineare di tipo reazione-diffusione) accoppiata ad un sistema di equazioni differenziali ordinarie nonlineari per calcolare le speci ioniche e le variabili di gating. In questa tesi scegliamo una discretizzazione spaziale basata sul metodo ad alto ordine degli elementi spettrali (SEM) perchè vantaggioso in accuratezza su questo problema matematico specifico; la soluzione può pertanto essere espressa come combinazione lineare di polinomi locali, elemento per elemento, pesati con i rispettivi gradi di libertà (DoFs). Introduciamo due discretizzazioni temporali differenti per l’equazione del monodominio, ovvero la Backward Differentiation Formula al secondo ordine (BDF2) e il metodo di Runge-Kutta al secondo ordine (RK2). La prima scelta porta ad uno schema temporale semi-implicito che produce, per ciascun istante temporale discreto, un sistema lineare che verrà risolto tramite il metodo del gradiente coniugato (CG). Al contrario, la seconda opzione ci permette di calcolare la soluzione a ciasun istante temporale discreto tramite un aggiornamento DoF per DoF, a patto che le formule di Legendre-Gauss-Lobatto (LGL) siano usate per risolvere gli integrali. Una discretizzazione spaziale basata su SEM con l’uso di formule LGL prende il nome di SEM con integrazione numerica (SEM-NI). Applicare SEM al modello ionico permette di operare un aggiornamento della soluzione DoF per DoF anche in questo caso, in quanto le equazioni discretizzate in spazio sono indipendenti tra loro. Infine discretizziamo in tempo il modello ionico usando BDF2. Diversi solver per l’elettrofisiologia cardiaca basati su CPU sono attualmente disponibili. L’impiego di processori grafici (GPU) si è dimostrato essere un’alternativa efficace, tuttavia in letteratura sono stati proposti pochi solver matrix-based basati su GPU per simulare la funzione cardiaca. In questa tesi vengono presentati due solver innovativi matrix-free, efficienti e scalabili, per l’elettrofisiologia cardiaca. Il primo è basato su GPU e nasce da una discretizzazione temporale semi-implicita, mentre l’altro è basato su CPU e presuppone una discretizzazione in tempo esplicita. Abbiamo condotto un’analisi accurata confrontando i diversi solver al variare del grado polinomiale delle funzioni base. Il solver GPU BDF2-BDF2 si è dimostrato più efficiente del solver CPU BDF2-BDF2 in un’esecuzione seriale, inoltre ottiene performance migliori per gradi polinomiali alti in un’esecuzione in parallelo. D’altra parte, il solver CPU RK2-BDF2, pur non richiedendo l’assemblaggio di un sistema lineare, è più lento del solver CPU BDF2-BDF2 siccome vengono valutate diverse quantità più volte per ciascun istante temporale discreto.