Newton–Krylov-BDDC deluxe solvers for non-symmetric fully implicit time discretizations of the bidomain model

In the last decade, the increasing need for understanding the intrinsic mechanisms behind cardiac diseases has resulted in the growth of inter-disciplinary studies, where, for example, physiological phenomena are translated into mathematical models [3, 13, 36, 40, 47]. Modern medicine employs computational tools and models to simulate, predict and analyze risk situations with a non-invasive approach. An increasing number of studies, for example, have addressed the understanding of dysfunctions of the heart and the interaction between bio-electrical and mechanical phenomena [8, 9, 11, 12, 32].

However, simulation of these models represents a tough challenge for personal laptops, as a huge amount of computational resources is required for numerical calculations. Even more so if one wants to use realistic heart geometries, where there is the need to represent accurately all the many facets. For this reason, the employment of supercomputers and large scale system architectures in these studies has quickly spread (see for example works on cardiac mechanics [10, 24, 26]).

In this perspective, the present work seeks to design, theoretically analyze and validate numerically a Newton–Krylov solver for fully implicit time discretizations of the bidomain model, preconditioned by a Balancing Domain Decomposition by Constraints (BDDC) algorithm. This model consists of a degenerate parabolic system of two non-linear reaction diffusion Partial Differential Equations (PDEs), which describes the propagation of the electric signal in the cardiac tissue [7, 38, 40]. This system is coupled, by means of a non-linear reaction term, with a model of ionic current flows and associated gating variables, modeled as a system of Ordinary Differential Equations (ODEs).

The main contribution of this paper is a novel theoretical analysis for the convergence rate bound of the preconditioned non-symmetric operator coming from a coupled approach for the solution of the non-linear system arising from a fully implicit time discretization of the cardiac electrical model.

Common alternatives in the literature use semi-implicit time discretizations [10, 49] and/or operator splitting [4, 5, 44], as fully implicit schemes are more expensive from a computational point of view if complex and high-dimensional non-linear ionic models (e.g. [14, 29, 45]) are coupled with the bidomain system. Other choices consist in decoupling strategies, where the ionic model is solved prior to the bidomain, for example in the previous work of the author [23] or in [15, 33, 34, 43]. In Ref. [35], an attempt at developing a solver for fully implicit time discretizations of the bidomain model was done, in the framework of additive Schwarz preconditioners. Nevertheless, fully implicit time discretizations allows to avoid constraints on the time steps in order to ensure stability of the numerical scheme. In this perspective, they are more flexible in case time adaptivity is also considered, see e.g. Ref. [33].

In this work we extend this solution strategy to the class of dual–primal Domain Decomposition (DD) algorithms, with particular focus on BDDC preconditioners.

BDDC preconditioners were introduced by [16] as an alternative to FETI-DP (Dual–Primal Finite Elements Tearing and Interconnecting [20]) methods for scalar elliptic problems and then analyzed by [30, 31]. Within this field of applications, BDDC algorithms have been employed for the solution of the linearized semi-implicit bidomain system in [49, 50] and for cardiac mechanics in Refs. [10, 37].

Instead of using a non-linear BDDC algorithm (such as the non-linear FETI-DP and BDDC proposed in [25]) or, more generally, non-linear preconditioning strategies (e.g. [28]), we present hereby Newton–Krylov-BDDC approach for the solution of a coupled solution strategy for fully implicit time discretizations of the bidomain system, including the ionic model. At each time step we solve and update a non-linear problem, where the Jacobian system arising from the linearization of the non-linear problem is non-symmetric, thus forcing us to use a Generalized Minimal Residual (GMRES) [42] method for its solution, preconditioned by BDDC algorithms in order to accelerate the convergence.

We propose a theoretical estimate of the convergence rate, based on the work in [18] (which provided a theoretical bound for the residual of the GMRES iterations) and in [48], whose BDDC preconditioner addressed the solution of non-symmetric systems arising from the discretization of advection–diffusion PDEs. This analysis is enriched by the employment of the recently-introduced deluxe scaling [17]. The robustness and efficiency of the proposed solver is then confirmed by extensive parallel numerical tests on the bidomain model, using the Portable, Extensible Toolkit for Scientific Computation (PETSc) library [1], thereby encouraging further investigations with realistic heart geometries and the tailoring of these kinds of solvers for the solution of the electro-mechanical model.

We highlight that a practical advantage of these methods with respect to the more popular algebraic multigrid preconditioner (e.g. [39]) as well as other DD preconditioners such as Multilevel Schwarz (which has been considered in many previous works e.g. [33, 34, 43]) is that they can be easily extended to unstructured meshes, since they do not need the implementation of an inter-grid operator, thus allowing the applicability to realistic heart geometries.

The work is structured as follows. In Sect. 2, we introduce the bidomain system, describing the propagation of the electric signal in the cardiac tissue. In Sect. 3 we give an insight into the space discretization and we formulate the fully implicit time scheme. Moreover we provide some properties related to the system arising from the discretization. A brief overview of non-overlapping DD spaces and objects as well as an introduction to BDDC preconditioner is provided in Sect. 4. The novel convergence rate estimate is then proved in Sect. 5, followed by extensive parallel numerical tests in Sect. 6.

2.1 Bidomain model

We consider here the macroscopic bidomain representation of the cardiac tissue, which is represented as two interpenetrating domains [7, 38]. These two anisotropic continuous media, named intra- and extracellular domains, are assumed to coexist at every point of the cardiac tissue and to be connected by a distributed continuous cellular membrane which fills the complete volume.

The cardiac tissue consists of a setting of fibers that rotates counterclockwise and that is arranged in laminar sheets running radially from the epi- to the endocardium (the outer and inner surface of the heart respectively).

At each point \(\textbf{x}\) of the cardiac domain \(\varOmega \) it is possible to define an orthonormal triplet of vectors \(\textbf{a}_l(\textbf{x})\) parallel to the local fiber direction, \(\textbf{a}_t(\textbf{x})\) and \(\textbf{a}_n(\textbf{x})\) tangent and orthogonal to the laminar sheets respectively and transversal to the fiber axis ([27]). Moreover, if \(\sigma _{l, t, n}^{i,e}\) are conductivity coefficients in the intra- and extracellular domain along the corresponding direction, it is possible to define the conductivity tensors \(D_i\) and \(D_e\) of the two media as

which describe the anisotropy of the intra- and extracellular media. For our theoretical purpose, we assume that \(\sigma _{l,t,n}^{i,e}\) are constant in space.

With these premises, we can obtain the parabolic-parabolic formulation of the bidomain model as the following non-linear parabolic reaction-diffusion system,

being \(u_i\) and \(u_e\) the intra and extracellular electric potentials, \(v = u_i-u_e\) the transmembrane potential and w the gating variables. The system is defined in \(\varOmega \times [ 0,T ]\), where \(\varOmega \) is the computational domain (representing an approximation of the left ventricle tissue), and [0, T] is the time interval considered.

Moreover, \(C_m\) is the membrane capacitance for unit area of the membrane surface, \(i_\text {ion}\) the ionic membrane current for unit area of the membrane surface, \(\chi \) is the membrane surface to volume ratio and \(I^{i,e}_\text {app}\) are the intra- and extracellular applied currents. The equations describing the propagation of the electric signal through the cardiac tissue are coupled through the reaction term to a system of Ordinary Differential Equations (ODEs) which describes the ionic currents flowing inward and outward the cell membrane. Regarding the boundaries, we assume that the heart is electrically insulated by requiring zero-flux boundary conditions

and compatibility condition \( \int _{\varOmega } I_{\text {app}}^i dx = \int _{\varOmega } I_{\text {app}}^e dx, \) where \(I_{\text {app}}^{i,e}\) are the intra- and extracellular applied currents, with initial values

See [6] for results on existence, uniqueness and regularity of the solution of (1).

2.2 Ionic current model

In this work we consider a phenomenological ionic model, derived from a modification of the renowned FitzHugh–Nagumo model [21, 22]: indeed, the Roger–McCulloch ionic model [41] overcomes the hyperpolarization of the cell during the repolarization phase by adding a non-linear dependence between the transmembrane potential and the gating. In this case, \(i_\text {ion}(v,w)\) and R(v, w) are given by

where G, \(v_{th}\), \(v_p\), \(\eta _1\) and \(\eta _2\) are constant coefficients. From now on, except ambiguity occurs, we consider the notation \(I_{\text {ion}}(v,w) = \chi i_\text {ion}(v,w)\).

3.1 Weak formulation and space discretization

Let the cardiac domain \(\varOmega \subset \mathbb {R}^3\) be a bounded open Lipschitz set with Lipschitz continuous boundary. Consider the functional space \( V = H^1(\varOmega ) \) and define the elliptic bilinear form associated with the intra- and extracellular conductivity tensors

Then, find \(u_{i,e} \in L^2(0,T; V) \) and \(w \in L^2(0,T; L^2(\varOmega ))\) such that \(\forall t \in (0,T)\)

\(\forall {\hat{u}}_{i,e} \in V\) and \(\forall {\hat{w}} \in L^2(\varOmega )\).

The model (2) is discretized in space by the finite element method, where the domain \(\varOmega \) is discretized by a structured quasi-uniform grid of hexaedral isoparametric \(Q_1\) elements. Denote by \(V_h \subset V\) be the associated finite element space, with the same basis functions \(\left\{ \varphi _p \right\} _{p=1}^{N_h}\) for all variables \(u_{i,e}\) and w and denote by \(A_{i,e}\) and M be the stiffness and mass matrices with entries

With these choices, we thus need to solve at each time step, the semi-discrete bidomain model

where

For simplicity, from now on we write \( \sum _l := \sum _{l = 1}^{N_h} \) assuming that we are adding contributions from all the \(N_h\) nodes of the discretization.

3.2 Fully implicit time scheme

In the literature, common alternatives takes into account implicit-explicit (IMEX) time discretization schemes [10, 49], where the diffusion term is treated implicitly while the remaining terms are treated explicitly, or, more generally, operator splitting [4, 5, 44]. Others effective strategies rely on a decoupling strategy (see e.g. Refs. [15, 23, 33, 34, 43]) where at each time step the microscopic and macroscopic models are solved successively. We propose here a fully implicit time discretization of the bidomain system, in the same fashion as in [35].

This choice allows larger time steps and, more generally, no constraints on the time step in order to guarantee stability. Nevertheless, the implicit treatment of the diffusion term and/or backward schemes for the parabolic part are widely popular and a well-established choice in the literature; moreover, since the ionic sodium currents are responsible for the stiffness of the ODEs, these terms can be included as well in the implicit discretization, allowing to overcome stability issues.

At the n-th time step,

1. 1.

   compute the intra- and extracellular potentials as well as the gating by solving the non-linear system \(F_{bido}(u_i^{n+1}, u_e^{n+1}, w^{n+1}) = 0\) derived from the Backward Euler scheme applied to the bidomain system,

   $$\begin{aligned}&F (u_i^{n+1}, \ u_e^{n+1}, \ w^{n+1}) := \begin{pmatrix} F_1 (u_i^{n+1}, u_e^{n+1}, w^{n+1}) \\ F_2 (u_i^{n+1}, u_e^{n+1}, w^{n+1}) \\ F_3 (u_i^{n+1}, u_e^{n+1}, w^{n+1}) \end{pmatrix} \nonumber \\&\quad =\begin{pmatrix} \chi C_m \left( v^{n+1}, \varphi _i \right) + \tau a_i \left( u_i^{n+1}, \varphi _i \right) + \tau \left( I_{\text {ion}}(v^{n+1}, w) , \varphi _i \right) \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad - \left[ \chi C_m \left( v^n, \varphi _i \right) + \tau \left( I_{\text {app}}, \varphi _i \right) \right] \\ -\chi C_m \left( v^{n+1}, \varphi _e \right) + \tau a_e \left( u_e^{n+1}, \varphi _e \right) - \tau \left( I_{\text {ion}}(v^{n+1}, w) , \varphi _e \right) \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad - \left[ \chi C_m \left( v^n, \varphi _e \right) + \tau \left( I_{\text {app}}, \varphi _e \right) \right] \\ \left( w^{n+1}, \varphi _w \right) - \tau \left( R(v^{n+1}, w^{n+1}), \varphi _w \right) - \left( w^n, \varphi _w \right) \end{pmatrix}, \end{aligned}$$

   (4)

   where \(\varphi _i, \varphi _e\) and \(\varphi _w\) are the test functions related to \(u_i\), \(u_e\) and w respectively.

   1. 2.1

      Apply an exact Newton method for the solution of the nonlinear system (4); given the initial guess \((u_i^0, u_e^0, w^0)\), at the \(k^{th}\) iteration of the Newton loop, solve the linear system of equations

      $$\begin{aligned} {\left\{ \begin{array}{ll} \sum _l \dfrac{\partial F_1}{\partial u_{i,l}} (u_i^k, u_e^k, w^k) s_{i,l}^{k+1}+ \sum _l \dfrac{\partial F_1}{\partial u_{e, l}} (u_i^k, u_e^k, w^k) s_{i,l}^{k+1}\\ \qquad \qquad \qquad \qquad + \sum _l \dfrac{\partial F_1}{\partial w_l} (u_i^k, u_e^k, w^k) s_{i,l}^{k+1}= -F_1(u_i^k, u_e^k, w^k) \\ \sum _l \dfrac{\partial F_2}{\partial u_{i,l}} (u_i^k, u_e^k, w^k) s_{e,l}^{k+1}+ \sum _l \dfrac{\partial F_2}{\partial u_{e, l}} (u_i^k, u_e^k, w^k) s_{e,l}^{k+1}\\ \qquad \qquad \qquad \qquad + \sum _l \dfrac{\partial F_2}{\partial w_l} (u_i^k, u_e^k, w^k) s_{e,l}^{k+1}= -F_2(u_i^k, u_e^k, w^k) \\ \sum _l \dfrac{\partial F_3}{\partial u_{i,l}} (u_i^k, u_e^k, w^k) s_{w,l}^{k+1}+ \sum _l \dfrac{\partial F_3}{\partial u_{e, l}} (u_i^k, u_e^k, w^k) s_{w,l}^{k+1}\\ \qquad \qquad \qquad \qquad + \sum _l \dfrac{\partial F_3}{\partial w_l} (u_i^k, u_e^k, w^k) s_{w,l}^{k+1}= -F_3(u_i^k, u_e^k, w^k) \end{array}\right. }, \end{aligned}$$

      (5)

      explicitly written as

      $$\begin{aligned} {\left\{ \begin{array}{ll} \chi C_m \left( s_i^{k+1}- s_e^{k+1}, \varphi _i \right) + \tau \left( \sum _l \frac{\partial I_{\text {ion}}}{\partial v_l} (v^k,w^k)\left( s_{i,l}^{k+1}- s_{e,l}^{k+1}\right) \psi _l, \varphi _i \right) \\ \qquad \qquad + \tau a_i \left( s_i^{k+1}, \varphi _i \right) + \tau \left( \sum _l \frac{\partial I_{\text {ion}}}{\partial w_l} (v^k,w^k)s_{w,l}^{k+1}\ \psi _l, \varphi _i \right) \\ \qquad \qquad = - F_1(u_i^k, u_e^k, w^k) \\ -\chi C_m \left( s_i^{k+1}- s_e^{k+1}, \varphi _e \right) - \tau \left( \sum _l \frac{\partial I_{\text {ion}}}{\partial v_l} (v^k,w^k)\left( s_{i,l}^{k+1}- s_{e,l}^{k+1}\right) \psi _l, \varphi _e \right) \\ \qquad \qquad + \tau a_e \left( s_e^{k+1}, \varphi _e \right) - \tau \left( \sum _l \frac{\partial I_{\text {ion}}}{\partial w_l} (v^k,w^k)s_{w,l}^{k+1}\ \psi _l, \varphi _e \right) \\ \qquad \qquad = - F_2(u_i^k, u_e^k, w^k) \\ \left( s_w^{k+1}, \varphi _k \right) - \tau \left( \sum _l \frac{\partial R}{\partial v_l} (v^k,w^k)\left( s_{i,l}^{k+1}- s_{e,l}^{k+1}\right) \psi _l , \varphi _w \right) \\ \qquad \qquad - \tau \left( \sum _l \frac{\partial R}{\partial w_l} (v^k,w^k)s_w^{k+1}\ \psi _l, \varphi _w \right) = - F_3(u_i^k, u_e^k, w^k) \end{array}\right. } \end{aligned}$$

      where

      $$\begin{aligned} s_i^{k+1}= \sum _l s_{i,l}^{k+1} = \delta u_i^{k+1}, \quad{} & {} s_e^{k+1}= \sum _l s_{e,l}^{k+1} = \delta u_e^{k+1},\\ {}{} & {} s_w^{k+1}= \sum _l s_{w,l}^{k+1} = \delta w^{k+1}\end{aligned}$$

      are the increments at time step k and \(\psi _l\) the l-th nodal basis function. In matricial form, this means to solve the linear system

      $$\begin{aligned} \textbf{JF}^k \textbf{s}^{k+1} = - \textbf{F} (\textbf{u}^k) \end{aligned}$$

      (6)

      where, by implying \(\dfrac{\partial g^k}{\partial t} = \dfrac{\partial g}{\partial t} (\textbf{v}^k, \textbf{w}^k) \), with \(g = \left\{ I_{\text {ion}}, R \right\} \) and \(t = \left\{ v, w \right\} \),

      $$\begin{aligned} \textbf{JF}^k= & {} \begin{bmatrix} \chi C_m M + \tau A_i + \tau M \dfrac{\partial I_{\text {ion}}^k}{\partial v} &{}- \chi C_m M - \tau M \dfrac{\partial I_{\text {ion}}^k}{\partial v} &{} \tau M \dfrac{\partial I_{\text {ion}}^k}{\partial w} \\ - \chi C_m M - \tau M \dfrac{\partial I_{\text {ion}}^k}{\partial v} &{}\chi C_m M + \tau A_e + \tau M \dfrac{\partial I_{\text {ion}}^k}{\partial v} &{}-\tau M \dfrac{\partial I_{\text {ion}}^k}{\partial w} \\ -\tau M \dfrac{\partial R^k}{\partial v} &{} \tau M \dfrac{\partial R^k}{\partial v} &{} \left( 1 - \tau \dfrac{\partial R^k}{\partial w} \right) M \end{bmatrix}, \\ \textbf{s}^{k+1}= & {} \begin{pmatrix} \mathbf {s_i}^{k+1} \\ \mathbf {s_e}^{k+1} \\ \textbf{w}^{k+1} \end{pmatrix}, \qquad \textbf{F} (\textbf{u}^k) = \begin{pmatrix} - M F_1(\mathbf {u_i}^k, \mathbf {u_e}^k, \textbf{w}^k) \\ -M F_2( \mathbf {u_i}^k, \mathbf {u_e}^k , \textbf{w}^k) \\ -M F_3( \mathbf {u_i}^k, \mathbf {u_e}^k , \textbf{w}^k) \end{pmatrix}, \end{aligned}$$

      with the same stiffness and mass matrices defined in (3).

   2. 2.2

      Update

      $$\begin{aligned} u_i^{k+1}= u_i^k+ s_i^{k+1}, \qquad u_e^{k+1}= u_e^k+ s_e^{k+1}, \qquad w^{k+1}= w^k+ s_w^{k+1}. \end{aligned}$$

We drop the index k from now on, unless an explicit ambiguity occurs.

Remark 1

In some cases, the construction of the Jacobian system (6) can be difficult, since several ionic models are not differentiable and complex. Instead of computing analytically the system, possible alternatives take in consideration the calculation of the derivatives via automatic differentiation or via finite differences.

3.3 Properties of the symmetric part of the bilinear form associated with the bidomain Jacobian system

The Jacobian system \(\textbf{JF}\) in (6) is non-symmetric, due to the inclusion of the ionic model. For this reason, the iterative solver must be addressed to the solution of such type of systems, such as the Generalized Minimal Residual (GMRES) method [42].

Following the work of [48] (where BDDC preconditioners are applied to the solution of non-symmetric problems arising from the discretization of advection–diffusion PDEs) in order to properly prove the convergence rate estimate of the solver, we need to associate to the linear system (6) a bilinear form and to analyze its symmetric and skew-symmetric parts.

We reformulate problem (5) in variational form: find \(s = \left( s_i, s_e, s_w\right) \in \mathbf {V_h}\), being \(\mathbf {V_h} = V_h \times V_h \times V_h\), such that

where

being \(\psi _l\) the l-th nodal basis function. The symmetric and skew-symmetric parts of \(a(s, \phi )\) respectively are denoted by

In the same way the system of linear equations (5) correspond to the finite element problem (7), we can denote by B and Z the symmetric and skew-symmetric parts of \(\textbf{JF}\), which correspond to the bilinear forms \(b(\cdot , \cdot )\) and \(z(\cdot , \cdot )\) respectively:

where we simplify the notation by writing \(\frac{\partial g}{\partial t} := \frac{\partial g}{\partial t} (\textbf{v}, \textbf{w})\), with \(g = \left\{ I_{\text {ion}}, R \right\} \) and \(t = \left\{ v, w \right\} \).

In the same spirit as in [23, 33], it is possible to show that \(b(\cdot , \cdot )\) is continuous and coercive with respect to an appropriate norm.

Lemma 1

Assume that

\(c_{1,2} \in \mathbb {R}^+\) and \(\forall l = 1, \dots , N\). Then the bilinear form \(b(\cdot , \cdot )\) is continuous and coercive with respect to the norm \(||| \cdot |||_{\tau , coup} \), defined as

\(\forall u = (u_1, u_2, u_3) \in \mathbf {V_h}\).

Remark 2

The norm \(||| \cdot |||_{\tau , coup} \) is well defined, as the quantity \(1-\tau \) is always positive (typical computational values for \(\tau \) are less than \(10^{-2}\)).

Remark 3

As in the case of the decoupled strategy (see [23, 33]), the hypothesis of non-negativity of the above Lemma is always satisfied for any time step \(\tau \le 0.37\) ms if we consider the Roger-McCulloch ionic model. Indeed, numerical computations of \(\chi C_m + \tau \frac{\partial I_\text {ion}}{\partial v}\) validate this assumption (see Fig. 1, above). Regarding the other two hypothesis, it is easy to compute analytically that \( 1 - \tau \frac{\partial R}{\partial w_l} = 1 + \eta _2 \tau \ge 0\) for any value of \(\tau \), being \(\eta _2\) a physiological parameter, while the last inequality is always satisfied for any \(v \geqslant 2 \cdot 10^{-4}\). This request is not restrictive, as for those values of the transmembrane potential the tissue is almost at rest.

Surface plots of \(\chi C_m + \tau \frac{\partial I_\text {ion}}{\partial v}\) (top) and \(\frac{\partial I_{\text {ion}}}{\partial w} - \frac{\partial R}{\partial v}\) (bottom), with \(C_m = 1 \frac{mF}{cm^3}\), \(\chi = 1\) and \(\tau =0.05\) ms, which are values usually employed in numerical experiments. Others parameters related to the Rogers–McCulloch ionic models are \(G=1.2\), \(v_{th}=13\), \(v_{p}=100\), \(\eta _1=4.4\) and \(\eta _2=0.012\)

Full size image

As an immediate consequence of the continuity and coercivity of the symmetric bilinear form \(b(\cdot , \cdot )\), it is possible to prove the following bounds.

Lemma 2

Assuming that the conductivity coefficients are constant in space, the bilinear form \(b(\cdot , \cdot )\) satisfies the bounds

where

and \(K_{M, *}\), \(K_{m, *}\), independent from the subdomain diameter H and the mesh size h.

Remark 4

This result is extensible to the case of conductivity coefficients almost constant over each subdomain.

4.1 Non-overlapping dual–primal algorithms

Let us decompose the cardiac domain \(\varOmega \) into N non-overlapping subdomains \(\varOmega _j\), with \(j=1,\dots ,N\), such that \( \overline{\varOmega } = \cup _{j=1}^{N} \overline{\varOmega }_j, \varOmega _j \cap \varOmega _k = \emptyset , \text {if } j \ne k, \) and the intersection between boundaries from different subdomains is either empty, a vertex, an edge or a face. We define the interface \(\varGamma \) as the set of points that belong to at least two subdomains,

where \(\partial \varOmega _j\), \(\partial \varOmega _k\) are the boundaries of \(\varOmega _j\) and \(\varOmega _k\) respectively. To our purposes, we assume the subdomains to be shape-regular with a typical diameter of size H; moreover let them be the union of shape-regular finite elements of diameter h.

Denote the associated local finite element spaces by \(W_j\) and partition it into its interior part \(W_I^{(j)}\) and the finite element trace space \(W_\varGamma ^{(j)}\), such that

In this dissertation, we consider variables on the Neumann boundaries \(\partial \varOmega _N\) as interior to a subdomain. By introducing the product spaces by

we define \({\widehat{W}} \subset W\) as the subspace of functions of W, which are continuous in all interface variables between subdomains. Similarly we denote by \({\widehat{W}}_\varGamma \subset W_\varGamma \), the subspace formed by the continuous elements of \(W_\varGamma \).

The idea of dual–primal algorithms is to solve iteratively in the space W, while imposing continuity constraints (primal constraints). Let \({\widetilde{W}}\) be the space of finite element functions in W, continuous in all primal variables, such that \({\widehat{W}} \subset {\widetilde{W}} \subset W\) and likewise \({\widehat{W}}_\varGamma \subset {\widetilde{W}}_\varGamma \subset W_\varGamma \).

Let \(W_\varPi ^{(j)} \subset W_\varGamma ^{(j)}\) be the primal subspace of functions, continuous across the interface, that will be subassembled between the subdomains that share \(\varGamma ^{(j)}\). We denote by dual, the subspace \(W_\varDelta ^{(j)} \subset W_\varGamma ^{(j)}\) which contains the finite element functions that can be discontinuous across the interface and which vanish at the primal degrees of freedom. Let \(W_\varPi \) and \(W_\varDelta \) be two subspaces such that

and \( W_\varGamma = W_\varPi \oplus W_\varDelta . \)

Schematic diagram of dual–primal subdivision between subdomains of interior (blue squares), dual (red circles) and primal (green triangle) degrees of freedom

Full size image

With this notation, it is possible to decompose \({\widetilde{W}}_\varGamma \) into a primal subspace \({\widehat{W}}_\varPi \) which has continuous elements only and a dual subspace \(W_\varDelta \) which contains finite element functions which are not continuous,

and in the same fashion

In this paper we will denote with subscripts I, \(\varDelta \) and \(\varPi \) the interior, the dual and the primal variables respectively. A schematic diagram can be found in Fig. 2.

In the non-overlapping framework, the global system matrix Eq. (6) is never formed explicitly, but a local matrix with the same structure is assembled on each subdomain, by restricting the integration set and by defining the local bilinear forms

and the symmetric and skew-symmetric counterparts

where \(\left( \cdot \ , \cdot \right) _{\varOmega _j}\) denotes the restriction of the \(L^2\)-inner product to the j-th subdomain.

As the proposed theory allows constant non-negative distribution of the diffusion coefficients among all subdomains, with large jumps aligned to the interfaces, these definitions are valid here.

Following the workflow in [48], after we introduce the space \({\widetilde{W}}\) of partially subassembled finite element space, we can define the corresponding bilinear forms by

We denote the partially subassembled matrices corresponding to the bilinear forms above with \({\widetilde{A}}\), \({\widetilde{B}}\) and \({\widetilde{Z}}\) respectively, and let

being \({\widetilde{R}}\) the injection operator from \({\widehat{W}}\) to \({\widetilde{W}}\).

In the same fashion as in [48], we define the truncated norms on the space \({\widetilde{W}}\)

In this work \(|| w ||_{L^2(\varOmega )} \) and \(| w |_{H^1(\varOmega )} \) for \(w \in {\widetilde{W}}\) always represent these truncated norms. As for construction the bilinear forms \(b^{(j)} (\cdot , \cdot )\) for \(j = 1, \dots , N\) are symmetric and positive definite on \(W^{(j)}\), it is possible to define

and

In dual–primal methods, the reordering of the degrees of freedom, lead to consider a reordered system matrix: assuming that the system (6) can be written as \(\mathscr {K} u = f\), then it is equivalent to write

where \( K_{II} = diag \left[ K_{II}^{(1)}, \dots , K_{II}^{(N)} \right] \) is a block-diagonal matrix. As in many iterative substructuring algorithms, we eliminate all the interior variables (step known as static condensation), obtaining the local Schur complement \(S_{\varGamma }^{(j)}\) on the j-th subdomain \(\varOmega _j\)

By defining the block-diagonal matrix \( S_\varGamma = diag \left[ S_\varGamma ^{(1)}, \dots , S_\varGamma ^{(N)} \right] \), \( {\widehat{S}}_\varGamma = R_\varGamma ^T S_\varGamma R_\varGamma , \) and the quantities \( g = \left[ g_I, \ {\widehat{g}}_\varGamma \right] ^T, \) \( {\widehat{g}}_\varGamma = R_\varGamma ^T g_\varGamma \), where \(R_\varGamma \) is the direct sum of local restriction operators \(R_\varGamma ^{(j)}\) (which returns the local interface components), \(S_\varGamma \) the unassembled global Schur complement system and g the right-hand side of the linear system, the resulting system which we need to solve is

Once this problem is solved, it is possible to retrieve the solution on the internal degrees of freedom (dofs) by using \(u_\varGamma \)

In the case of this application, the Schur complement system Eq. 8 is non-symmetric, thus it is necessary to apply a solver for non-symmetric problems, such as the GMRES iterative algorithm. For any \(u_\varGamma \in {\widetilde{W}}_\varGamma \), we define the harmonic extension to the interior of subdomains \(u_{\mathscr {A}, \varGamma }\) as

and analogously, it is possible to define its counterpart \(u_{\mathscr {A}, \varGamma }\in {\widehat{W}}\) for \(u_\varGamma \in {\widehat{W}}_\varGamma \).

We define the following bilinear forms for vectors in \({\widehat{W}}_\varGamma \) and \({\widetilde{W}}_\varGamma \)

We observe, as needed for further calculations, that [48, Lemma 7.2] follows from these definitions.

Following [48], it is useful to define \(B_\varGamma \) and \({\widetilde{B}}_\varGamma \) norms:

4.2 Preconditioning

As already mentioned, with the solution strategy described previously, the Schur complement matrix \({\widehat{S}}_\varGamma \) of the Jacobian bidomain system (6) is non-symmetric but positive semidefinite: therefore it is necessary to apply a solver for non-symmetric systems, such as the Generalized Minimal Residual method (GMRES) [42].

Additionally, in order to enable fast convergence, preconditioning occurs. We hereby present theoretical results related to the Balancing Domain Decomposition with Constraints (BDDC) preconditioning algorithm [17].

When working with these methods, an interface averaging is needed: the standard scaling (\(\rho \)-scaling) has weights built from the values of the elliptic coefficients in each substructure. On the contrary, the stiffness-scaling takes its weights from the diagonal elements of both local and global stiffness matrix, while the more recent deluxe-scaling (see [2, 17]) is based on the solution of local problems built from local Schur complements associated with the dual unknowns.

As in our previous work [23], we provide here a convergence rate estimate that holds both with the classic \(\rho \)-scaling and with the deluxe-scaling.

4.2.1 Restriction operators and scaling.

Before going into details of the proposed preconditioners, we define the restriction operators

and the direct sums \(R_\varDelta = \oplus R_\varDelta ^{(j)}\), \(R_\varPi = \oplus R_\varPi ^{(j)} \) and \({\widetilde{R}}_\varGamma = R_{\varGamma \varPi } \oplus R_{\varGamma \varDelta }\), which maps \(W_\varGamma \) into \({\widetilde{W}}_\varGamma \). We also need a proper scaling of the dual variables.

The \(\rho \)-scaling, originally proposed for Neumann-Neumann methods, can be defined for the coupled bidomain model at each node \(x \in \varGamma ^{(j)}\) as

where \(\mathscr {N}_x\) is the set of indices of all subdomains with x in the closure of the subdomain. We note that \(\mathscr {N}_x\) induces the definition of an equivalence relation that classifies the interface degrees of freedom into faces, edges and vertices equivalence classes.

The deluxe scaling [2, 17] computes the average \({\bar{w}} = E_Dw\) for each face \(\mathscr {F}\) or edge \(\mathscr {E}\) equivalence class as follows: suppose that \(\mathscr {F}\) is shared by subdomains \(\varOmega _j\) and \(\varOmega _k\). Denote by \(S_{\mathscr {F}}^{(j)}\) and \(S_{\mathscr {F}}^{(k)}\) the principal minors obtained from \(S^{(j)}_\varGamma \) and \(S^{(k)}_\varGamma \) by extracting all rows and columns related to the degrees of freedom of the face \(\mathscr {F}\). Let \(u_{j, \mathscr {F}}= R_{\mathscr {F}}u_j\) be the restriction of \(u_j\) to the face \(\mathscr {F}\) through the restriction operator \(R_{\mathscr {F}}\). Then, the deluxe average across \(\mathscr {F}\) can be defined as

The action of \( ( S_{\mathscr {F}}^{(j)}+ S_{\mathscr {F}}^{(k)})^{-1} \) can be computed by solving a Dirichlet problem over the two subdomains involved, with zero value on the right-hand side entries that correspond with the interior degrees of freedom.

It is possible to extend this definition when considering the deluxe average across an edge \(\mathscr {E}\). Suppose for simplicity that \(\mathscr {E}\) is shared by only three subdomains with indices \(j_1\), \(j_2\) and \(j_3\); the extension to more than three subdomains is straightforward. Denote by \(u_{j, \mathscr {E}}= R_{\mathscr {E}}u_j\) the restriction of \(u_j\) to the edge \(\mathscr {E}\) through the restriction operator \(R_{\mathscr {E}}\) and define \(S_{\mathscr {E}}^{(j_{123})}= S_{\mathscr {E}}^{(j_1)}+ S_{\mathscr {E}}^{(j_2)}+ S_{\mathscr {E}}^{(j_3)}\); the deluxe average across an edge \(\mathscr {E}\) is given by

The relevant equivalence classes, involving the substructure \(\varOmega _j\), will contribute to the values of \({\bar{u}}\). These contributions will belong to \({\widehat{W}}_\varGamma \), after being extended by zero to \(\varGamma \backslash \mathscr {F}\) or \(\varGamma \backslash \mathscr {E}\); the sum of all contributions will result in \(R^T_*{\bar{u}}_*\). We then add the contributions from the different equivalence classes to obtain

where \(E_D\) is a projection. We define its complementary projection by

We define the scaling matrix for each subdomain \(\varOmega _j\)

being \(k_1, \dots , k_j \in \varXi _j^*\), a set containing the indices of the subdomains that share the face \(\mathscr {F}\) or the edge \(\mathscr {E}\) and where the diagonal blocks are given by \(D^{(j)}_\mathscr {F} = ( S_{\mathscr {F}}^{(j)}+ S_{\mathscr {F}}^{(k)})^{-1} S_{\mathscr {F}}^{(j)}\) or \(D^{(j)}_\mathscr {E}= (S_{\mathscr {E}}^{(j_1)}+ S_{\mathscr {E}}^{(j_2)}+ S_{\mathscr {E}}^{(j_3)})^{-1} S_{\mathscr {E}}^{(j_1)}\).

Lastly, we define the scaled local restriction operators

\(R_{D, \varDelta }\) as direct sum of \(R_{D, \varDelta }^{(j)}\) and the global scaled operator \({\widetilde{R}}_{D, \varGamma } = R_{\varGamma \varPi } \oplus R_{D, \varDelta } R_{\varGamma \varDelta }\).

4.2.2 BDDC preconditioner

We recall that the Jacobian linear system (6) for the original non-linear reaction-diffusion problem has been reduced to the non-symmetric Schur complement system \({\widehat{S}}_\varGamma u_\varGamma = {\widehat{f}}_\varGamma \) for the subdomain interface variables.

BDDC algorithms were initially proposed by [17] for the solution of symmetric, positive definite problems, but its formulation can be equally applied to non-symmetric problems, as done in [48] for advection–diffusion equations.

We define the BDDC preconditioner using the restriction operators as

Thanks to this definition, BDDC algorithm for the solution of the Schur complement problem (8) can be defined as: find \(u_\varGamma \in {\widehat{W}}_\varGamma \) such that

The above system is non-symmetric, thus a GMRES iterative method is employed for its solution. Once the interface \(u_\varGamma \) is computed, we can retrieve the internal solution \(u_I\) by solving local Dirichlet problems.

The key-point of the proof for the convergence rate estimate relies in this Lemma, which can be also proved for the \(\rho \)-scaling.

Lemma 3

Assume that the primal space is spanned by the vertex nodal finite element functions and the edge cutoff functions. Let the projection operator be scaled by either the standard \(\rho \)-scaling or the deluxe-scaling. Then

holds \(\forall u \in {\widetilde{W}}_\varGamma \), with \(n = 2\) in case the \(\rho \)-scaling is applied, \(n=3\) in case of the deluxe-scaling.

Proof

We report here a sketch of the proof for the deluxe scaling.

We first define the deluxe scaling induced by the symmetric bilinear form \(b(\cdot , \cdot )\). Let the face \(\mathscr {F}\) be shared by \(\varOmega _j\) and \(\varOmega _k\); denote with \(B_{\mathscr {F}}^{(j)}\) and \(B_{\mathscr {F}}^{(k)}\) be the principal minors obtained from \(B^{(j)}_\varGamma \) and \(B^{(k)}_\varGamma \) by extracting all the contributions that are related to the degrees of freedom of the face \(\mathscr {F}\).

Let \(u_{j, \mathscr {F}}= R_{\mathscr {F}}u_j\) be the restriction of \(u_j\) to the face \(\mathscr {F}\) through the restriction operator \(R_{\mathscr {F}}\); we define the deluxe average induced by the bilinear form \(b(\cdot , \cdot )\) across \(\mathscr {F}\) as

We have denoted with \({\bar{u}}_{j,\mathscr {F}} \) the mean value of \(u_j\) over the face \(\mathscr {F}\), while, if we omit the index j, we denote the deluxe average.

In a similar manner, we can define the deluxe average induced by the symmetric bilinear form \(b(\cdot , \cdot )\) across an edge \(\mathscr {E}\).

Suppose for simplicity that \(\mathscr {E}\) is shared by only three subdomains with indices \(j_1\), \(j_2\) and \(j_3\); the extension to more than three subdomains is immediate. Let \(u_{j, \mathscr {E}}= R_{\mathscr {E}}u_j\) be the restriction of \(u_j\) to the edge \(\mathscr {E}\) through the restriction operator \(R_{\mathscr {E}}\) and define

then, the deluxe average induced by the bilinear form \(b(\cdot , \cdot )\) across an edge \(\mathscr {E}\) is given by

An initial crucial point, as usual in the substructuring framework, is to estimate the local contributions for the complementary projection \(P_D\),

where \(\varXi _j^*\) is the index set containing the indices of the subdomains that share the face \(\mathscr {F}\) or the edge \(\mathscr {E}\). We will denote by

the vector containing the mean value of the intra- and extra-cellular potentials over \(\mathscr {G} = \{ \mathscr {F}, \mathscr {E} \}\) on the subdomain j and null value corresponding to the gating component. Regarding the face contributions, by simple algebra, it is easy to bound the quantity \(u_{j, \mathscr {F}} - {\bar{u}}_\mathscr {F}\) with

\(\forall u_{j, \mathscr {F}} \in {\widetilde{W}}_\varGamma \), by using the inequalities that arise from the generalized eigenvalue problem \(B_{\mathscr {F}}^{(j)}\phi = \lambda B_{\mathscr {F}}^{(k)}\phi \) and by observing that all eigenvalues are strictly positive [2]

It is sufficient to estimate \(| u_{j, \mathscr {F}} - {\bar{u}}^{i,e,w}_{j,\mathscr {F}} |_{B_\mathscr {F}^{(j)}} ^2\) and \(| ( B_{\mathscr {F}}^{(j)}+ B_{\mathscr {F}}^{(k)})^{-1} B_{\mathscr {F}}^{(k)}( {\bar{u}}^{i,e,w}_{j,\mathscr {F}} - {\bar{u}}^{i,e,w}_{k,\mathscr {F}} ) |_{B_\mathscr {F}^{(j)}} ^2\); we highlight that, in case also the subdomain faces averages are included in the primal space, the latter is zero. The first term can be bounded by applying the ellipticity Lemma 2, the Poincarè-Friedrichs inequality and the Trace theorem, by

where \(\mathscr {H}^{i,e}_{j} u_j\) is the discrete Laplacian extension operator. Regarding the second term \(| ( B_{\mathscr {F}}^{(j)}+ B_{\mathscr {F}}^{(k)})^{-1} B_{\mathscr {F}}^{(k)}( {\bar{u}}^{i,e,w}_{j,\mathscr {F}} - {\bar{u}}^{i,e,w}_{k,\mathscr {F}} ) |_{B_\mathscr {F}^{(j)}} ^2\), let \(\mathscr {E} \subset \partial \mathscr {F}\) be a primal edge, such that \({\bar{u}}^{i,e,w}_{j,\mathscr {E} } = {\bar{u}}^{i,e,w}_{k,\mathscr {E} } \). Then, it is straightforward to see that

Combining the result of ellipticity Lemma (2), the Poincarè-Friedrichs inequality and the Trace theorem, we obtain

as the mean value of the gating component vanishes by construction Eq. 15 and the norm related to the potentials can be obtained with results from Ref. [46, Chapter 4, Lemmas 4.26 and 4.30]. To conclude, the face contribution for the bound of \(| P_D u |_{B_\varGamma ^{(j)}} ^2\) is given by

In a similar manner, the edge contribution can be obtained with

where for simplicity, by supposing that an edge \(\mathscr {E}\) is shared only by three substructures, each with indexes \(j_1\), \(j_2\) and \(j_3\) (the extension to the case of more subdomains is then similar), we define \( B_{\mathscr {E}}^{(j_{123})}:= B_{\mathscr {E}}^{(j_1)}+ B_{\mathscr {E}}^{(j_2)}+ B_{\mathscr {E}}^{(j_3)}\) and the average operator as

Proceeding in the same fashion as for the face contribution, it follows

By adding and subtracting \({\bar{u}}^{i,e,w}_{j_1,\mathscr {E} } \) (which assume the same value over the three subdomain, as we have included the edge averages into the primal space) we can get the counterpart estimate for the egdes:

In conclusion, the edge contribution for the estimate of \(| P_D u |_{B_\varGamma ^{(j)}} ^2\) is given by

where the index j collects all contributions from the subdomains that share edge \(\mathscr {E}\). \(\square \)

The convergence rate of the preconditioned GMRES iteration can be obtained using the result in [18] and following the proof techniques proposed in [48].

Theorem 1

Let H be the subdomain size and let the mesh size h be small enough.

Assume, for \(u \in {\widehat{W}}_\varGamma \), that there exists two positive constants c and C such that

hold, with \(c = \dfrac{c_0}{K^2} \) and \(C = \varPhi ^{\star , k} (H,h) \ K^2 \), where

where \(n=2,3\) depends on a standard \(\rho \)-scaling or deluxe scaling procedure and T is the preconditioned operator \( T = M^{-1}_{BDDC} {\widetilde{S}}_\varGamma = {\widetilde{R}}_{D, \varGamma }^T {\widetilde{S}}_\varGamma ^{-1} {\widetilde{R}}_{D, \varGamma } {\widetilde{R}}_\varGamma S_\varGamma {\widetilde{R}}_\varGamma ^T . \) Then

where \(r_m\) is the residual at the m-th iteration.

5.1 Proof of the upper bound

For the proof of the upper bound \(\langle Tu, Tu \rangle _{B_\varGamma } \le C \langle u, u \rangle _{B_\varGamma }\), we need the following results.

Lemma 4

There exists a constant \(C_1 >0\) such that \(\forall u_j, v_j \in W^{(j)}\) with \(j=1,\dots ,N\),

where

Proof

Thanks to the definition of the \(B^{(j)}\)-norm and by using the ellipticity Lemma 2, we can bound from below the norm

Therefore, \(\forall u_j, v_j \in W^{(j)}\)

from which

We can estimate the bound for the skew-symmetric bilinear form

The bound for the bilinear form \(a (\cdot , \cdot )\) follows easily from its decomposition

and from the continuity of both symmetric and skew-symmetric forms. \(\square \)

Lemma 5

There exists a constant \(C_2 > 0 \) such that \(\forall u, v \in {\widehat{W}}\)

with K defined in Lemma 4.

Lemma 6

Let \(C_3, C_4 > 0 \) be two positive constants, independent from H and h, such that \(\forall u_\varGamma , v_\varGamma \in {\widetilde{W}}_\varGamma \)

where K is defined in Lemma 4.

Proof

1. (i)

   Regarding the first bound, it is straightforward to see that \(\forall u_\varGamma , v_\varGamma \in {\widetilde{W}}_\varGamma \)

   $$\begin{aligned} | \langle u_\varGamma , v_\varGamma \rangle _{{\widetilde{Z}}_\varGamma } | \le K | u_{\mathscr {A}, \varGamma } |_{{\widetilde{B}}} | v_{\mathscr {A}, \varGamma } |_{{\widetilde{B}}} = K | u_\varGamma |_{{\widetilde{B}}_\varGamma } | v_\varGamma |_{{\widetilde{B}}_\varGamma } , \end{aligned}$$

   by applying Definition 10 and Lemma 4.

2. (ii)

   By using [48, Lemma 7.2] and 4, we obtain

   $$\begin{aligned} | \langle u_\varGamma , v_\varGamma \rangle _{{\widetilde{S}}_\varGamma } | = | \langle u_{\mathscr {A}, \varGamma }, v \rangle _{{\widetilde{A}}} | = | {\widetilde{a}} \left( u_{\mathscr {A}, \varGamma }, v \right) | \le K | u_\varGamma |_{{\widetilde{B}}_\varGamma } | v_\varGamma |_{{\widetilde{B}}_\varGamma } . \end{aligned}$$

\(\square \)

We are ready to prove the upper bound of Theorem 1.

Given \(u_\varGamma \in {\widehat{W}}_\varGamma \), by defining \(w_\varGamma = {\widetilde{S}}_\varGamma ^{-1} {\widetilde{R}}_{D, \varGamma } u_\varGamma \) and using [48, Lemmas 7.2 and 7.9] and Lemmas 3 and 6

from which we conclude

where

and \(n=2,3\) depends on which scaling we are considering.

5.2 Proof of the lower bound

Conversely, for the proof of the lower bound, we need Lemma 7.9 from [48] and the following results.

Lemma 7

There exists a constant \(C_5 > 0\) such that \(\forall u_\varGamma , v_\varGamma \in {\widehat{W}}_\varGamma \),

with K defined in Lemma 4.

Lemma 8

Given \(w_\varGamma = {\widetilde{S}}_\varGamma ^{-1} {\widetilde{R}}_{D, \varGamma } S_\varGamma u_\varGamma \) for \(u_\varGamma \in {\widehat{W}}_\varGamma \), there exists a constant \(C_6 > 0\) such that

with \(\varPhi \) and K defined in Eq. (16).

Proof

We make use of [48, Lemma 7.2] and Lemma 3 to obtain

where in the last inequality we combined Lemma 7.9 from [48] and Lemma 6. It follows

It is straightforward to conclude, by applying again [48, Lemma 7.9] and Lemma 6

\(\square \)

Lemma 9

Let \(w_\varGamma = {\widetilde{S}}^{-1}_\varGamma {\widetilde{R}}_{D, \varGamma } S_\varGamma u_\varGamma \), for \(u_\varGamma \in {\widehat{W}}_\varGamma \). Then, the following property holds

for all \(v \in {\widetilde{R}} ( {\widehat{W}} )\).

Proof

Given \(v \in {\widetilde{R}} ( {\widehat{W}} )\), we denote \(v_\varGamma \in {\widetilde{R}}_\varGamma ( {\widehat{W}}_\varGamma )\) its continuous interface part. Then, given \(u_\varGamma \in {\widehat{W}}_\varGamma \), by using Lemma 7.2 from [48] and the identity \({\widetilde{R}}_{\varGamma } {\widetilde{R}}_{D, \varGamma }^T v_\varGamma = v_\varGamma \), we get

\(\square \)

Lemma 10

For h sufficiently small, given \(w_\varGamma = {\widetilde{S}}_\varGamma ^{-1} {\widetilde{R}}_{D, \varGamma } S_\varGamma u_\varGamma \) for \(u_\varGamma \in {\widehat{W}}_\varGamma \), there exists a positive constant \(C_7 >0\) such that \(\forall u_\varGamma \in {\widehat{W}}_\varGamma \)

with \(\varPhi \) and K are defined in Eq. (16).

Proof

It is useful to observe that, being \(u_\varGamma \in {\widehat{W}}_\varGamma \), by defining \({\widetilde{R}}: {\widehat{W}} \longrightarrow {\widetilde{W}}\), then \({\widetilde{R}} u_\varGamma \in {\widetilde{R}} \left( {\widehat{W}} \right) \). By using the Poincaré-Friedrichs inequality, the ellipticity Lemma 2 and [48, Lemma 7.2], we can obtain

We focus on the right-hand side term:

We treat separately the two terms.

1. (i)

   Using Lemma 8 and [48, Lemma 7.2],

   $$\begin{aligned} \begin{aligned} \langle w_{\mathscr {A}, \varGamma }- {\widetilde{R}} u_{\mathscr {A}, \varGamma }, w_{\mathscr {A}, \varGamma }\rangle _{{\widetilde{B}}}&= | w_\varGamma |_{{\widetilde{B}}_\varGamma } ^2 - | {\widetilde{R}} u_{\mathscr {A}, \varGamma } |_{{\widetilde{B}}} ^2 = | w_\varGamma |_{{\widetilde{B}}_\varGamma } ^2 - | u_\varGamma |_{B_\varGamma } ^2 \\&\le C \left[ K^2 \varPhi ^{\star , k} (H,h) - 1 \right] | u_\varGamma |_{B_\varGamma } ^2 , \end{aligned} \end{aligned}$$

   where we observe that, for \(u_{\mathscr {A}, \varGamma }\in {\widehat{W}}\), \({\widetilde{R}} u_{\mathscr {A}, \varGamma }\in {\widetilde{R}} ({\widehat{W}})\) and, since \({\widetilde{B}}\) is symmetric, Lemma 9 holds:

   $$\begin{aligned} \langle {\widetilde{R}} u_{\mathscr {A}, \varGamma }, w_{\mathscr {A}, \varGamma }\rangle _{{\widetilde{B}}} = \langle w_{\mathscr {A}, \varGamma }, {\widetilde{R}} u_{\mathscr {A}, \varGamma }\rangle _{{\widetilde{B}}} = \langle {\widetilde{R}} u_{\mathscr {A}, \varGamma }, {\widetilde{R}} u_{\mathscr {A}, \varGamma }\rangle _{{\widetilde{B}}}. \end{aligned}$$
2. (ii)

   From Lemma 9, it holds

   $$\begin{aligned} \langle w_{\mathscr {A}, \varGamma }- {\widetilde{R}} u_{\mathscr {A}, \varGamma }, \ {\widetilde{R}} u_{\mathscr {A}, \varGamma }\rangle _{{\widetilde{B}}} = 0, \end{aligned}$$

   since \({\widetilde{R}} u_{\mathscr {A}, \varGamma }\in {\widetilde{R}} ( {\widehat{W}} )\).

Then

which leads to

\(\square \)

Lemma 11

Given \(v_\varGamma = {\widetilde{R}}_{D, \varGamma }^T w_\varGamma \) for \(w_\varGamma \in {\widetilde{W}}_\varGamma \), there exists a positive constant \(C >0\) such that

with \(\varPhi \) defined in Eq. (16).

Proof

As the quantity \(v_{\mathscr {A}, \varGamma }- w_{\mathscr {A}, \varGamma }\) has zero average on the interface, \(\forall w_\varGamma \in {\widetilde{W}}_\varGamma \), then thanks to Poincaré-Friedrichs inequality, the ellipticity Lemmas 2, 3 and [48, Lemma 7.2] we have

where we used the definition of projection \(E_D\),

and an inverse inequality. \(\square \)

Lemma 12

Set \(v_\varGamma = T u_\varGamma - u_\varGamma \), for \(u_\varGamma \in {\widehat{W}}_\varGamma \). For h sufficiently small, there exists a positive constant \(C>0\) such that for \(u_\varGamma \in {\widehat{W}}_\varGamma \),

with K, \(\varPhi \) defined in Equations (16).

Proof

Since \(T u_\varGamma = {\widetilde{R}}^T_{D, \varGamma } w_\varGamma \) and \({\widetilde{R}}^T_{D, \varGamma } {\widetilde{R}}_\varGamma = I\), then \( v_\varGamma = {\widetilde{R}}^T_{D, \varGamma } \left( w_\varGamma - {\widetilde{R}}_\varGamma u_\varGamma \right) . \) Therefore, thanks to Lemmas 11 and 10,

\(\square \)

Finally, we are able to conclude the proof of Theorem 1 by proving the lower bound of Theorem 1. By using Lemma 6 and [48, lemmas 7.2 and 7.9],

which means

Since [48, Lemma 7.2] and \(\langle u_\varGamma , u_\varGamma \rangle _{Z_\varGamma } = 0 \), we obtain

Therefore, (17) becomes

However, thanks to Lemmas 7 and 12,

from which it follows

In conclusion, we have the lower bound

where

We now validate our theoretical results with several parallel numerical experiments, performed on the supercomputer Galileo from the Cineca centre (www.hpc.cineca.it), a Linux Infiniband cluster equipped with 1084 nodes, each with 36 2.30 GHz Intel Xeon E5-2697 v4 cores and 128 GB/node, for a total of 28,804 cores. The last set of tests have been perfomed on Indaco cluster from University of Milan, a Linux Infiniband cluster with 16 nodes, each carrying 2 processors Intel Xeon E5-2683 V4 2.1 GHz with 16 cores each, for a total amount of 512 cores (https://www.indaco.unimi.it/). Our C code is based on the PETSc library [1] from the Argonne National Laboratory.

The numerical tests are run both on a thin Cartesian slab and on a idealized left ventricular geometry, modeled as a truncated ellipsoid (see Fig. 3).

Portion of left ventricle, idealized as a truncated ellipsoid (left) and slab geometry (right)

Full size image

The latter is described in ellipsoidal coordinates by the parametric equations

where \(a (r) = a_1 + r (a_2 - a_1)\), \(b (r) = b_1 + r (b_2 - b_1)\) and \(c (r) = c_1 +r (c_2 - c_1)\) with \(a_{1,2}\), \(b_{1,2}\) and \(c_{1,2}\) given coefficients defining the main axes of the ellipsoid.

We assume that fibers rotate intramurally linearly with the depth for a total amount of \(120^o\) proceeding counterclockwise from epicardium (\(r = 1\), outer surface of the truncated ellipsoid) to endocardium (\(r=0\), inner surface).

We choose the Roger-McCulloch ionic model for the description of the ionic current flows, whose physiological parameters on Table 1 are taken from [41].

Regarding the bidomain physiological coefficients in Table 2, we refer to [7].

In order to produce the propagation of the electric potential, we apply for 1 ms to the surface of the domain representing the endocardium an external stimulus of \(I_{\text {app}}= 100\) mA/cm\(^3\). Instead, on the slab geometry, the current is applied in one corner of the domain, over a spheric volume of radius 0.1 cm.

We report in Figs. 4 and 5 the time evolution of the transmembrane and extra-cellular potentials respectively from the endocardial side of the portion of the truncated ellipsoid. The external stimulus \(I_{\text {app}}\) is applied on five different sites on the endocardium layer, positioned at the apex of the idealized left ventricle.

The boundary conditions are for an insulated tissue, while the initial conditions represent a resting potential. The time step is fixed \(\tau = 0.05\) ms.

The values reported in the Tables (average number of Newton iterations, average number of linear iterations) are mean values over the time interval and are rounded at the nearest integer. We adopt a classical implementation of the Newton method, with residual stopping criterion with tolerance \(10^{-4}\). The non-symmetric linear system arising from the discretization of the Jacobian problem at each Newton step is solved with the Generalized Minimal Residual (GMRES) method, preconditioned by the BDDC preconditioner (included in the PETSC library^{Footnote 1}) and the Boomer Algebraic MultiGrid (bAMG, from the Hypre library [19]).

Snapshots (every 5 ms) of transmembrane potential v time evolution. For each time frame, we report the epicardial view of a portion of the left ventricle, modeled as a truncated ellipsoid

Full size image

Snapshots (every 5 ms) of extra-cellular potential \(u_e\) time evolution. For each time frame, we report the epicardial view of a portion of the left ventricle, modeled as a truncated ellipsoid

Full size image

Test 1: weak scaling. The first set of tests we report here is a weak scaling test on both slab and ellipsoidal domains. For both cases, we fix the local mesh size to \(12\cdot 12\cdot 12\) and we increase the number of subdomains from 32 to 256, thus resulting in an increasing slab geometry and in an increasing portion of ellipsoid. In this way, the dofs are increasing from 180k up to 1 million and a half. From Tables 3 and 4, it is evident how the dual–primal algorithm has a better performance respect to the bAMG: the average number of linear iteration per Newton iteration (lit) does not increase with the number of subdomains and is clearly lower. As a matter of fact, for the slab geometry we can observe an increasing reduction rate from \(85\%\) up to \(93\%\) for the average linear iterations, while for the ellipsoid it varies between \(65\%\) and \(90\%\). In contrast, BDDC’s average computational time is higher due to the need of communication between the implemented structures (Fig. 6).

Weak scaling test. Slab (left) and ellipsoidal (right) geometries, local mesh of \(12 \cdot 12 \cdot 12\) elements. Simulations of 2 ms of cardiac activation with \(dt = 0.05\) ms (40 time steps). Comparison of average CPU time per time step, in seconds

Full size image

Test 2: strong scaling. We perform a strong scaling test for the two geometries: we fix the global mesh to \(128 \cdot 128 \cdot 24\) elements (resulting in more than 1 millions of dofs) and we increase the number of subdomains from 32 to 256.

We can observe from Table 5 that, as the local number of dofs decreases, the preconditioner with the better balance in term of average linear iterations and CPU time per time step is BDDC: although slightly higher CPU times, the performance is balanced by the average number of linear iterations, which is certainly lower. We test the efficiency of the proposed solver on the parallel architecture by computing the parallel speedup \(= \frac{T_1}{T_N}\), which is the ratio between the runtime needed by 1 processor (\(T_1\)) and the average runtime needed by N processors (\(T_N\)) to solve the problem. In both cases, bAMG preconditioner is sub-optimal (see Fig. 7), while BDDC outperforms the ideal speedup. This super-ideal strong scaling might be explained by the strong reduction of the computational costs required by the PETSc LU factorization of the local problems, when they become smaller due to the increase of the number of processors.

Strong scaling test with global mesh of \(128\cdot 128\cdot 24\). Slab (left) and ellipsoidal (right) geometries. Simulations of 2 ms of cardiac activation with \(dt = 0.05\) ms, for a total amount fo 40 time steps. Comparison of actual parallel speedup (ideal speedup dotted)

Full size image

Test 3: optimality tests. Tables 6 and 7 report the results of optimality tests, for both slab and ellipsoid geometries, carried on Galileo cluster.

We fix the number of processors (subdomains) to \(4 \cdot 4 \cdot 4\) and we increase the local size H/h from 4 to 24, thus reducing the finite element size h. We consider both scalings (\(\rho \)-scaling on top, deluxe scaling at the bottom of each table) and we test the solver for increasing primal spaces: V includes only vertex constraints, V+E includes vertex and edge constraints, and V+E+F includes vertex, edge and face constraints. The deluxe scaling tests are up to a local size of \(20 \cdot 20 \cdot 20\) elements, due to limited computational resources. We consider a time interval of 2 ms during the cardiac activation phase. The time step is \(dt = 0.05\) ms, for a total amount fo 40 time steps.

Similar results hold for both geometries. The deluxe solver seems to be more robust while increasing the local mesh size, both in terms of average linear iterations (see also Fig. 8 ) and average CPU time per time step.

Optimality test with different scalings (dash-dotted \(\rho \)-scaling, continuous deluxe scaling) and primal sets (V = vertices, E = edge averages, F = face averages). Slab (left) and ellipsoidal (right) domains, \(4 \cdot 4 \cdot 4\) subdomains, increasing local size from \(4 \cdot 4 \cdot 4\) to \(24 \cdot 24 \cdot 24\). Average number of linear iterations per time step

Full size image

Full activation-recovery simulation on slab domain, time interval [0, 150]ms, 3000 time steps (on the left) and on ellipsoidal domain, time interval [0, 200]ms, 4000 time steps (on the right). Fixed number of subdomains \(8 \cdot 8 \cdot 2\) and fixed global mesh \( 128 \cdot 96 \cdot 24\). Comparison between bAMG and BDDC average number of linear iterations per time step. Performed on Indaco cluster

Full size image

Test 4: full activation-recovery interval simulations. In these last tests, we compare the performance of our dual–primal and the multigrid preconditioners during a full activation-recovery interval. These tests have been performed on Indaco cluster from the University of Milan.

We fix the number of subdomains to \(128=8 \cdot 8 \cdot 2\) and the global mesh size to \(128 \cdot 96 \cdot 24\), thus considering local problems of 8,619 dofs. We consider a time interval of [0, 200] ms for a total of 4000 time steps for a portion of ellipsoid defined by \(\varphi _{\min } = -\pi /2\), \(\varphi _{\max } = \pi /2\), \(\theta _{\min } = -3/8 \pi \) and \(\theta _{\max } = \pi /8\), while on the slab of dimensions \(1.92 \times 0.96 \times 0.96\) cm\(^3\) we perform the tests for 3000 time steps, on the time interval [0, 150] ms.

In Fig. 9 we report the trend of the average number of linear iteration per time step during the simulation. Firstly, we notice a huge difference between the bAMG and the BDDC preconditioners, with a reduction of more than \(85\%\) for the latter. The number of iterations remains bounded during the test.

Both preconditioned solvers seem to affected by the different stages of the action potential: an initial peak during the activation phase is followed by a constant elevated number of linear iterations as the electric signal propagates in the cardiac tissue, ending with a lower - but always constant - number of linear iterations during the resting phase.

Despite similar qualitative trends between the two geometries, it is undeniable that there are differences from the quantitative point of view, due to the complexity of the domain taken in consideration. Even if we observe comparable performances in terms of CPU time per time step (see Table 8) of both preconditioners, we notice a huge difference in the assembly time for the Jacobian matrix and for the right-hand side of the linearized system. In case of bAMG, the solver spends between \(9.67-14.91\%\) of the total CPU time assembling the Jacobian matrix, while \(0.8-1.25\%\) for the right-hand side, whether we are considering the ellipsoidal or the slab geometry. On the other hand, BDDC is much faster for these processes, since it takes only \(0.86\%\) and \(0.11-0.15\%\) of the total CPU time in order to compute the Jacobian matrix and the right-hand side vectors respectively, for both geometries.

We have designed a dual–primal Newton–Krylov solver for the monolithic solution strategy for fully implicit time discretizations of the cardiac bidomain model. Theoretical analysis for the convergence rate of the non-symmetric preconditioned operator has been provided. We have validated this result through extensive parallel numerical tests, showing the efficiency and robustness of the solver, thus encouraging further investigation using more complex ionic models and realistic heart geometries.

Further research would be needed to verify whether the current theoretical analysis and numerical validation could be extended to more advanced ionic models and realistic geometries, that represent the major limitations of the present study.

Additionally, possible future discussion should be devoted to the testing of alternative nonlinear solvers, such as NGMRES, as well as inexact Newton methods, in order to further improve convergence.