Automated adjoints of coupled PDE-ODE systems
Patrick E. Farrell Johan E. Hake Simon W. Funke Marie E. Rognes



Straight quantum layer  with impurities inducing resonances  
     Sylwia Kondej
    
=============================================================




Mathematical models that couple partial differential equations (PDEs)
and spatially distributed ordinary differential equations (ODEs) arise
in biology, medicine, chemistry and many other fields. In this paper
we discuss an extension to the FEniCS finite element software for
expressing and efficiently solving such coupled systems. Given an ODE
described using an augmentation of the Unified Form Language (UFL) and
a discretisation described by an arbitrary Butcher tableau, efficient
code is automatically generated for the parallel solution of the
ODE. The high-level description of the solution algorithm also
facilitates the automatic derivation of the adjoint and tangent
linearization of coupled PDE-ODE solvers. We demonstrate the
capabilities of the approach on examples from cardiac
electrophysiology and mitochondrial swelling.



  finite element methods, adjoints, coupled PDE-ODE, FEniCS, dolfin-adjoint, code generation



  65L06, 65M60, 65M32, 35Q92





§ INTRODUCTION



In this paper we discuss solvers for systems involving spatially
dependent ordinary differential equations (ODEs) of the following
form: given an initial condition y_0(x), find y = y(x, t) such
that

    y(x, t) = f(y, x, t),    y(x, t_0) = y_0(x),

for all x in a point set X ⊂Ω⊆^d, where
d is the spatial dimension. The subscript t refers to
differentiation in time.  The right-hand side f cannot depend on
spatial derivatives of y; the ODE is decoupled at different points.
Problems of this form often arise when discretizing time-dependent
mathematical models that couple PDEs with spatially distributed
systems of ODEs via operator splitting. Examples of application areas
include cardiac electrophysiology in
general <cit.> and cardiac ion channel
modeling in particular <cit.>, mitochondrial
swelling <cit.>, groundwater flow and
contamination <cit.>, pulmonary gas
transport <cit.>, and
plasma-enhanced chemical vapor deposition <cit.>.
We also discuss the automated derivation of the adjoint and tangent
linearization of such models:
these
can be used to identify the sensitivity of the solution to model parameters,
solve inverse problems for unknown parameters, and characterize
the stability of trajectories.

For concreteness, we present two biological examples of coupled
PDE-ODE systems where problems of the form (<ref>)
arise. We return to numerical results for these examples in
Section <ref>.

[The bidomain equations]
As our first motivating example, we will consider the bidomain
equations for the propagation of an electrical signal in a
non-deforming domain Ω <cit.>: find the
transmembrane potential v = v(x, t), the extracellular potential
u_e(x, t) and additional state variables s = s(x, t) such that for
t ∈ (0, T]:

  

    v -  (M_i  v + M_i  u_e)    = - I_ ion(v, s) + I_s
         in Ω, 
    
      -  (M_i  v + (M_i + M_e)  u_e  )    = 0
         in Ω, 
    s   = F(v, s)
         ∀ x ∈Ω.


In (<ref>), I_ ion is a given nonlinear function
describing ionic currents and F defines a system of nonlinear
functions, while M_i and M_e are the intracellular and
extracellular conductivity tensors, respectively, and I_s is a given
stimulus current. The function F cannot depend on spatial
derivatives of v and s; it defines a pointwise system of ODEs. The
specific form of I_ ion and F are typically prescribed by a
given cardiac cell model and may vary greatly in complexity:
from involving a single state variable s such as the FitzHugh-Nagumo
model <cit.> to models with e.g. 41 state variables
such as the model of O'Hara et al. <cit.>. The
system (<ref>) is closed with appropriate initial and
boundary conditions. After the application of operator splitting, the
ODE step is decoupled from the PDE step and is a system of the form of
(<ref>).


[Mitochondrial swelling]
As a second example, we consider a model proposed
in <cit.> to describe the swelling of
mitochondria. Mitochondrial swelling plays a key role in the process
of programmed cell death (apoptosis) and thus for the life
cycle of cells. Mathematically, we consider the following
model <cit.>: find the calcium concentration u
= u(x, t) and the densities of mitochondria states N_i = N_i(x, t)
for i = 1, 2, 3 such that

  
  
    u   = d_1 Δ (|u|^q-2 u) + d_2 g (u) N_2
           in Ω, 
    N_1   = - f(u) N_1
           ∀ x ∈Ω, 
    N_2   = f(u) N_1 - g(u) N_2
           ∀ x ∈Ω, 
    N_3   = - g(u) N_2
           ∀ x ∈Ω.


where d_1 ≥ is a diffusion coefficient, d_2 ≥ 0 is a
feedback parameter, q ≥ 2 ∈ℕ determines the
nonlinearity of the diffusion-type operator, and g and f are
prescribed functions. The functions N_1, N_2 and N_3
describe the densities of unswollen, swelling and completely
swollen mitochondria respectively. Equation (<ref>) is a
spatially-coupled PDE, while equations
(<ref>)-(<ref>) define pointwise ODEs. The system
is closed with initial conditions and Dirichlet boundary
conditions for u. Again, after operator splitting, a subproblem of
the form (<ref>) results.


Over the last decade, there has been a growing interest in
computational frameworks for the rapid development of numerical solvers
for PDEs such
as the FEniCS Project <cit.>, the Firedrake
Project <cit.>, and
Feel++ <cit.>. The rapid development of
solvers is achieved by offering high-level abstractions
for the expression of such problems. Each of these projects
provides a domain-specific language for specifying finite element
variational formulations of PDEs and associated software for their
efficient solution. The methods presented in this paper are
implemented in the the FEniCS Project and
dolfin-adjoint <cit.>; dolfin-adjoint automatically
derives discrete adjoint and tangent linear models from a FEniCS
forward model. These allow for the efficient computation of functional
gradients and Hessian-vector products, and are essential ingredients in stability
analyses, parameter identification and inverse problems.

However, these systems do not currently efficiently extend to the kind
of coupled PDE-ODE systems arising in computational biology.  While
monolithic finite element discretizations of coupled PDE-ODE systems
such as (<ref>) or (<ref>) may easily be
specified and solved using the current software features in FEniCS,
this approach will involve an unreasonably large computational
expense, as the monolithic system is highly nonlinear and the solver
cannot exploit the fact that the ODE is spatially decoupled.

Operator splitting is the method of choice for such coupled PDE-ODE
systems <cit.>, and is implemented as standard in
many hand-written codes, see
e.g. <cit.> in the context of cardiac
electrophysiology. Operator splitting decouples a PDE-ODE system into
a spatially-coupled system of PDEs and a spatially-decoupled
collection of ODE systems. This collection of ODE systems then
typically takes the form (<ref>), which can be solved
using well-known temporal discretization methods. Heretofore it has
not been possible to specify spatially-decoupled ODE systems in
high-level PDE frameworks such as the FEniCS or Firedrake projects.

This work addresses the gap in available abstractions, algorithms
and software for arbitrary PDE-ODE systems.
We introduce high-level
domain-specific language constructs for specifying collections of ODE
systems and for specifying multistage ODE schemes via Butcher
tableaux. This has several major benefits. By automatically generating
the solver from a high-level description of the problem, practitioners
can flexibly explore a range of models; this is especially important
in biological problems where the model itself is uncertain.  Another
advantage is that a high-level description facilitates the
automated derivation of the associated tangent linear and adjoint
models.  This represents a significant saving in the time taken to
investigate questions of biological interest.

The main new contributions of this paper are (i) the extension of the
FEniCS finite element system to enable efficient the large-scale
forward solution of coupled, time-dependent PDE-ODE systems via
operator splitting, and (ii) the extension of dolfin-adjoint to
automatically derive and solve the associated tangent linear and
adjoint models, enabling efficient automated computation of functional
gradients for use in e.g. optimization or adjoint-based sensitivity
analysis.

This paper is organized as follows. In
Section <ref>, we describe the general operator
splitting setting, the resulting separate PDE and ODE systems, and
natural discretizations of these. We continue in
Section <ref> by briefly describing the
well-established ODE schemes that we consider in this work: the
multistage and Rush-Larsen families. We discuss the adjoint and
tangent linear discretizations of a general operator splitting scheme
and derive the adjoint and tangent linear models for the multistage
and Rush-Larsen schemes in Section <ref>. Key
features of our new implementation are described in
Section <ref>. We present numerical results for two
different application examples in Section <ref>, and
use these to evaluate the performance of our implementation. In
Section <ref>, we provide some concluding remarks and
discuss current limitations and possible future extensions.




§ OPERATOR SPLITTING FOR COUPLED PDE-ODE SYSTEMS





 §.§ Operator splitting



A classical approach to the solution of the bidomain
system (<ref>) <cit.> and similar systems
is to apply operator splitting. This approximately solves the full
system of equations by alternating between the solution of a system of
PDEs and a system of ODEs defined over each time interval. The
advantage of this approach is that it decouples the solution of the
typically highly nonlinear ODEs from the spatially coupled PDEs at
each time step. We will consider this general setting, but
use (<ref>) as a concrete example.

Applying a variable order operator split to (<ref>), the
resulting scheme reads: given initial conditions v^0, s^0, an order
parameter θ∈ [0, 1], and time points {t_0, …, t_N}
with an associated timestep _n = t_n+1 - t_n, then for
each time step n = 0, 1, …, N - 1:


  * Compute v^∗ and s^∗ by solving the ODE system
  
    
    
    v   = - I_ ion (v, s) 
    s   = F(v, s)

  
  over Ω× [t_n, t_n + θκ_n] with
  initial conditions v^n, s^n.

  * Compute v^† and u^n+1_e by solving the PDE system
  
    
    
    v -  (M_i  v + M_i  u_e)    = I_s 
    
          -  (M_i  v + (M_i + M_e)  u_e)    = 0

  
  over Ω× [t_n, t_n+1] with initial condition
  v^∗.

  * If θ < 1, compute solutions v^n+1 and s^n+1
  solving (<ref>) over Ω× [t_n + θκ_n, t_n+1] with initial conditions v^† and
  s^∗.

The split scheme relies on the repeated solution of a
nonlinear system of ODEs and (in this case) a linear system of
PDEs. The main advantage of the approach is that it allows for the
separate discretization and solution of the ODEs and PDEs, with the
respective solutions as feedback into the other system.
For θ = 1/2, the resulting Strang splitting scheme is
second-order accurate; for other values of θ the
resulting scheme is first-order accurate.



 §.§ Discretization of the separate PDE and ODE systems


Suppose now that the relevant PDE system
(e.g. (<ref>)) is discretized in space by a finite
element method defined over a mesh _h of the domain Ω,
and in time by some suitable temporal discretization. Efficient
solution algorithms for such discretizations are well-established
(such as <cit.>) and will not be detailed further
here.

At each iteration the solution of the PDE system relies on the
solution of the ODE system (e.g. (<ref>)), or at
least the ODE solution evaluated at some finite set of points X =
{x_i}_i=1^|X| in space. For instance, the set of points X may
be taken as the nodal locations of the finite element degrees of
freedom, or the quadrature points of the mesh. As a consequence of this and the
spatial locality of ODEs, a natural approach to discretizing the
system of ODEs is to step the ODE system forward in time at this set
of points X. Typically, |X| is very large and the efficient
repeated solution of these systems of ODEs is key. This is the setting
that we focus on next.




§ SOLUTION SCHEMES FOR THE ODE SYSTEMS



The initial value problem (<ref>) decouples in space, and
henceforth we consider its solution for a fixed x. With a minor
abuse of notation, let y = y(x) and f(y, t) = f(y, x, t) so
that (<ref>) reads as the classical ODE problem: find y
∈ C^1([T_0, T_1]; ℝ^m) for a certain m ∈ℕ
such that

    y(t) = f(y, t),    y(T_0) = y_0.

for t ∈ [T_0, T_1]. There exists a wide variety of solution
schemes for (<ref>) including but not limited to multistage,
multistep, and IMEX schemes, see e.g. <cit.>. In this
work, we focus on two classes of schemes that are widely used in
computational biology: multistage schemes and so-called Rush-Larsen
schemes. The Rush-Larsen schemes are commonly used in cardiac
electrophysiology in general and for discretizations of the bidomain
equations (<ref>) in particular. These two classes of
schemes are detailed below. Keeping the general iterative operator
splitting setting in mind, we present the schemes on a single
time-step [T_0, T_1] with κ = T_1 - T_0 for brevity of
notation.



 §.§ Multistage schemes



An s-stage multistage scheme for (<ref>) is defined by a set
of coefficients a_ij, b_i and c_i for i, j = 1, …, s,
commonly listed in a so-called Butcher tableau; for more
details, see e.g. <cit.>. Given y^0 at T_0, the
scheme finds the stage variables k_i for i = 1, …, s
satisfying

    k_i = f(y^0 + κ∑_j=1^s a_ij k_j, T_0 + c_i κ),

and subsequently sets the solution y^1 at T_1 via

    y^1 = y^0 + κ∑_i=1^s b_i k_i.

Note that (<ref>) defines a system of (non-linear)
equations to solve for the stage variables k_i (i = 1, …, s)
if a_ij≠ 0 for any j ≥ i. Our implementation demands
that a_ij = 0 for j > i, i.e. does not allow the computation of
earlier stages to depend on the values of later stages.



 §.§ The Rush-Larsen scheme and its generalization

Recall that y = {y_i}_i=1^m and f(y) = {f_i(y)}_i=1^m. The
original Rush-Larsen scheme employs an exponential integration scheme
for all linear terms for each component f_i (for i = 1, 2, …,
M), and a forward Euler step for all non-linear
terms <cit.>. Let J_i = ∂
  f_i/∂ y_i(y^0, T_0) be the i'th diagonal component of the
Jacobian at time T_0 for i = 1, …, M, and assume that J_i
is non-zero. For brevity, denote f_i(y^0, T_0) = f_i^0. The
Rush-Larsen scheme then computes y_i^1 at T_1 by

    y_i^1 = y_i^0 + {[ J_i^-1 f_i^0 (e^κ J_i-1 )                         ,;                   κ f_i^0                         . ] .


For stiff systems of ODEs, the forward Euler step in (<ref>)
can become unstable for large κ. This motivates a generalized
version of the Rush-Larsen scheme <cit.>. In this
generalization, the exponential integration step is used for all
components of y reducing the scheme to

    y_i^1 = y_i^0 + J_i^-1 f_i^0 (e^κ J_i-1 )

Both (<ref>) and (<ref>) are first order accurate in
time <cit.>.

Both the original and the generalized Rush-Larsen schemes can be
developed into second order schemes by repeated use as follows. The
solution y^1/2 at T_1/2 = T_0 +
κ/2 is computed in the first step and is used in f and
its linearization J to compute y^1 in the second step. Let
f_i^1/2 = f_i(y^1/2, T_1/2), and let
J_i, 1/2 = ∂ f_i/∂
  y_i(y^1/2, T_1/2) be the diagonal component of
the Jacobian at time T_1/2 for i = 1, …, M. More
precisely, the second order version of (<ref>)
and (<ref>) is then given by

  
  
    y_i^1/2   = y_i^0 + {[ J_i^-1 f_i^0 (e^κ/2  J_i - 1 )                               ;                      κ/2 f_i^0                                ].
    
        y_i^1   = y_i^0 + {[ J_i, 1/2^-1 f_i^1/2 (e^κ J_i, 1/2 - 1 )                                        ;                               κ f_i^1/2                                         ].


for the original Rush-Larsen scheme and

  
  
    y_i^1/2   = y_i^0 + J_i^-1 f_i^0  (e^κ/2 J_i-1 ) 
    
        y_i^1   = y_i^0 + J_i, 1/2^-1 f_i^1/2( e^κ J_i, 1/2 - 1 )


for the generalized Rush-Larsen scheme. Here, (<ref>) is
simplified compared to the scheme in <cit.>: instead
of evaluating f and J at y̅^(i) (cf. eq. (7)
in <cit.>) we use y_i^1/2. 



§ ADJOINTS AND TANGENT LINEARIZATIONS OF COUPLED PDE-ODE SYSTEMS


Our aim is to automatically derive the adjoint and tangent linear
equations of operator splitting schemes in a manner that allows for
efficient solution of the resulting systems. We proceed as follows: we
first state the adjoint and tangent linear system for general
problems, then consider the special case of mixed PDE-ODE systems and
finally discuss the specific adjoint and tangent linear versions of
multistage and Rush-Larson schemes.

We begin by considering a general system of discretized equations in the form:
find y ∈ Y ⊂ℝ^M such that

    F(y) = 0,

The adjoint equation of (<ref>) associated with a real-valued
functional of interest J: Y → is: find the adjoint
solution y̅∈ℝ^M such that

    ∂ F/∂ y^∗(y) y̅ = ∂ J/∂ y(y)

where the superscript ∗ denotes the adjoint operator. The
associated tangent linear equation with respect to some auxiliary
control parameter m is given by: find the tangent linear solution
ẏ such that

    ∂ F/∂ y (y)ẏ
        = ∂ F/∂ m(y)


We now turn our attention to the case where F is a coupled PDE-ODE
system. For a guiding example, we consider again the bidomain
equations (<ref>)–(<ref>) and,
without loss of generality restrict ourselves to a single iteration
(N = 1).  For brevity, we denote the ODE operator
(<ref>) as O and the PDE operator
(<ref>) as P. We can then rewrite the scheme in the
form (<ref>):

    
    v_0 - v^0    = 0,       (Set initial condition)
    
            s_0 - s^0    = 0,       (Set initial condition)
    
            O(v^*, s^*; v^0, s^0)    = 0,       (Solve bidomain ODEs (<ref>))
    
            P(v^†, u_e^1; v^*)    = 0,       (Solve bidomain PDEs (<ref>))
    
            O(v^1, s^1; v^†, s^*)    = 0.      (Solve bidomain ODEs for final state)

    

The first pair of equations in (<ref>) set the
initial conditions, the following equations represents a collection of
nonlinear ODEs, the third equation a system of PDEs, and the last
equation again a collection of ODEs.

From (<ref>), the adjoint system for
(<ref>) derives as:

    (
        [                  2c|2*I  2|c|2*∂O^*/∂(v^0, s^0)                  2c|2*0                   2c2*0;                      xx                                                                                                                                                ;                  2c|2*0   2c|2*∂O^*/∂(v^*, s^*)   2c|2*∂P^*/∂(v^0, s^0)           2c2*∂O^*/∂s^*;                                                                                                                                                                                                ;                  2c|2*0                  2c|2*0 2c|2*∂P^*/∂(v^†, u_e^1)           2c2*∂O^*/∂v^†;                                                                                                                                                                                                ;                  2c|2*0                  2c|2*0                  2c|2*0    2c2*∂O^*/∂(v^1, s^1);                                                                                                                                                                                                 ])
      [   v̅^0;   s̅^0;   v̅^∗;   s̅^∗;   v̅^†; u̅_e^1;   v̅^1;   s̅^1;        ]
        =
      [         0;       2*⋮;          ;          ;          ;          ; ∂ J/∂ v^1; ∂ J/∂ s^1 ]

For brevity, the left hand side matrix combines 2×2 blocks, and
the functional of interest is assumed to only depend on the final
state variables v^1, s^1. The adjoint system
(<ref>) is a linear coupled PDE-ODE system with
upper-triangular block structure, which can efficiently by solved by
backwards substitution.  The last block row represents a linear
adjoint ODE system, preceded by a adjoint PDE system, and preceded by
another adjoint ODE system, and finalised by a variable assignment in
the first block row which represents the adjoint version of setting
the initial conditions.

The derivation of adjoint and tangent linear equations for finite
element discretizations of PDEs, such as the third block row in
(<ref>), is well-established, see
e.g. <cit.>. However, the automated derivation of
adjoint and tangent linear systems for the multistage and Rush-Larsen
discretizations of the generic ODE system (<ref>), such as the
second and fourth block row in (<ref>), is less
so, and these derivations are presented below.



 §.§ Adjoints and tangent linearizations of multistage schemes


Consider a multistage discretization of (<ref>) as described in
Section <ref> with a given Butcher tableau
a_ij, b_j, c_i, i, j = 1, …, s. Taking s = 3 for
illustrative purposes, we can
write (<ref>)–(<ref>) in form
(<ref>) as:

    [ y^0; k_1; k_2; k_3; y^1 ]
      -
      [                         y_0; f(y^0 + κ w_1, T_0 + c_1 κ); f(y^0 + κ w_2, T_0 + c_2 κ); f(y^0 + κ w_3, T_0 + c_3 κ);      y^0 + κ∑_i=1^3 b_i k_i ]
        =0,

with

    w_i = y^0 + κ∑_j=1^s a_ij k_j.

Assuming that the functional of interest J is independent of
the internal stage values k_i, we can derive the adjoint problem as:

    [ y̅^0; k̅_1; k̅_2; k̅_3; y̅^1 ]
      -
      [                   0       ∂ f_1/∂ w_1^*       ∂ f_2/∂ w_2^*       ∂ f_3/∂ w_3^*                   I;                   0 κ a_11∂ f_1/∂ w_1^* κ a_21∂ f_2/∂ w_2^* κ a_31∂ f_3/∂ w_3^*               κ b_1;                   0                   0 κ a_22∂ f_2/∂ w_2^* κ a_32∂ f_3/∂ w_3^*               κ b_2;                   0                   0                   0 κ a_33∂ f_3/∂ w_3^*               κ b_3;                   0                   0                   0                   0                   0;                     ][ y̅^0; k̅_1; k̅_2; k̅_3; y̅^1 ]
      =
      [ ∂ J/∂ y^0;         0;         0;         0;      y̅_0 ]
      ,

where y̅_0 is the terminal condition for the adjoint solution
at T_1.

Generalizing to s stages, we see that we first solve for the adjoint
stage values and then compute the adjoint solution at time T_0 via

  
  
    (I - κ a_ii∂ f_i/∂ w_i^*) k̅_i    = κ b_i y̅^1 + ∑_j=i + 1^s κ a_ji∂ f_j/∂ w_j^* k̅_j, 
    y̅^0    = y̅^1 + ∑_i=1^s ∂ f_i/∂ w_i^* k̅_i.



Following a similar calculation, the tangent linearization of the
multistage scheme with s stages is: given y_0 at T_0 compute
k̇_i for i = 1, …, s and then ẏ^1 at T_1 via

  
  
    (I - κ a_ii∂ f_i/∂ w_i) k̇_i     = ∂ f_i/∂ w_iẏ^n + ∑_j = 1^i-1κ a_ij∂ f_i/∂ w_ik̇_j + ∂ f_i/∂ m, 
    ẏ^n+1   = ẏ^n + ∑_i=1^s κ b_i k̇_i.





 §.§ Adjoints and tangent linearizations of Rush-Larsen schemes


We now turn to consider adjoints and tangent linearizations of the
(generalized) Rush-Larsen schemes. We here derive the equations for
the first order generalized Rush-Larsen scheme given
by (<ref>):

    y_i^1 = y_i^0 + J_i^-1 f_i^0 (e^κ J_i-1 )

for each component i = 1, …, M of the state variable.

We can write (<ref>) in form (<ref>) as

    [          y^0 - y_0; y^1 - L(y^0) - y^0;                    ]
    = 0

where y_0 is the given initial condition at T_0 and L(y^0) = { L_i(y^0) }_i=1^M with

    L_i(y^0)
      = J_i^-1 f_i^0 ( e^κ J_i-1  )
      ≡∂ f_i/∂ y_i(y^0, T_0)^-1 f_i(y_0, T_0) (e^κ∂ f_i/∂ y_i(y^0, T_0) - 1) .

The adjoint system of the first order generalised Rush-Larsen scheme
is then given by

    [ y̅^0; y̅^1;      ]
        -
      [                0 I +  ∂ L/∂ y^0^∗;                0                0;                  ][ y̅^0; y̅^1;      ]
    =
    [ ∂ J/∂ y^0;      y̅_0;           ]

where y̅_0 is a given terminal condition for the adjoint at
T_1. Similarly, the tangent linear system (<ref>)
with respect to an auxiliary parameter m is given by

    [ ẏ^0; ẏ^1 ]
        -
      [             0             0; I + ∂ L/∂ y^0             0;               ][ ẏ^0; ẏ^1 ]
      =
      [          ẏ_0; ∂ L(y^0)/∂ m ]

where ẏ_0 is a given initial condition at T_0.

The adjoint and tangent linear equations for the other Rush-Larsen
schemes are derived following the same steps, and we therefore omit
the details here.




§ FENICS AND DOLFIN-ADJOINT ABSTRACTIONS AND ALGORITHMS



The FEniCS Project defines a collection of software components
targeting the automated solution of differential equations via finite
element methods <cit.>. The components 
include the Unified Form Language (UFL) <cit.>, the
FEniCS Form Compiler (FFC) <cit.> and the finite
element library DOLFIN <cit.>. The separate
dolfin-adjoint project and software automatically derives the discrete
adjoint and tangent linear models from a forward model written in the
Python interface to DOLFIN <cit.>. This section
presents our extensions of the FEniCS form language and the FEniCS
form compiler, and other new FEniCS and dolfin-adjoint software
features targeting coupled PDE-ODE systems in general and collections
of systems of ODE in particular.



 §.§ Variational formulation of collections of ODE systems


Consider a spatial domain Ω⊂^d tessellated by a mesh
_h, and a collection of general initial problems as defined
by (<ref>) over a set of points X defined relative to
_h. For x_i ∈ X, denote by δ_x_i the Dirac delta
function centered at x_i ∈^d such that

    f(x_i)
      = ∫_^d f  δ_x_i

for all continuous, compactly supported functions f. We will also
write

    ∑_x_i ∈ X∫_^d f  δ_x_i≡∑_x_i ∈ X∫ f (x_i)
      ≡∑_x_i ∈ Xf1_x_i≡f1_X


Drawing inspiration from our context of finite element variational
formulations, we can then write (<ref>) as: find
y(·, t) such that

    y_t(t)ψ_i_X = f(y, t)ψ_i_X

for all (basis) functions (or distributions) ψ_i such that
ψ_i(x_i) = 1 and ψ_j(x_i) = 0 for j ≠ i.

The remainder of this section describes the extensions of the FEniCS
and Dolfin-adjoint systems to allow for in particular abstract
representation and efficient forward and reverse solution of systems
of the form (<ref>).



 §.§ Extending UFL with vertex integrals


UFL is an expressive domain-specific language for abstractly
representing (finite element) variational formulations of differential
equations. In particular, the language defines syntax for integration
over various domains. Consider a mesh _h of geometric
dimension d with cells {𝒯}, interior facets
{ℱ^i} and boundary facets {ℱ}^b. Interior
facets are defined as the d-1 dimensional intersections between two
cells and thus for any interior facet ℱ^i we can write
ℱ^i = 𝒯^+∩𝒯^-. UFL defines
the sum of integrals over cells, sum of integrals over boundary facets
and sum of integrals over interior facets by the , ,
and  measures, respectively. For example, the following linear
variational form defined in terms of a piecewise polynomial u
defined relative over _h:

    L(u) = ∑_{𝒯}∫_𝒯 u  
      + ∑_{ℱ^i}∫_ℱ^i u |_T^+
      + ∑_{ℱ^b}∫_ℱ^b u

is naturally expressed in UFL as:

L = u*dx + u('+')*dS + u*ds


To allow for point evaluation over the vertices of a mesh
(cf. (<ref>)), we have introduced a new
vertex integral (or vertex measure in UFL terms) type
with default instantiation :

from ufl import Measure
dP = Measure('vertex')
L = f*dP

In agreement with (<ref>), the vertex
integral is defined by:

    𝚏*𝚍𝙿≡f1_𝒱(_h) = ∑_x_i ∈𝒱(_h) f (x_i),

where V(_h) is the set of all vertices of the mesh
_h. Vertex measures restricted to a subset of vertices can be
defined as for all other UFL measure types.

This basic extension of UFL crucially allows for the abstract
specification of collections of ODEs such
as (<ref>) in a manner that is consistent and
compatible with abstract specification of finite element variational
formulations of PDEs. It also allows for the native specification of
e.g. point sources in finite element formulations of PDEs in FEniCS.



 §.§ Vertex integrals in the FEniCS Form Compiler


The FEniCS Form Compiler FFC generates specialized C++ code <cit.> from the symbolic UFL
representation of variational forms and finite element
spaces <cit.>. To accommodate the new vertex
integral type, we have extended the UFC interface with a class
 defining the interface for the tabulation of
the element matrix corresponding to the evaluation of an expression at
a given vertex, listed below:

/// This class defines the interface for the tabulation of
/// an expression evaluated at exactly one point.
class vertex_integral: public integral

public:

  /// Constructor
  vertex_integral() 

  /// Destructor
  virtual  vertex_integral() 

  /// Tabulate tensor for contribution from local vertex
  virtual void tabulate_tensor(double * A,
                               const double * const * w,
                               const double * coordinate_dofs,
                               std::size_t vertex,
                               int cell_orientation) const = 0;
;

The FFC code generation pipeline has been correspondingly extended to
allow for the generation of optimized code from the UFL
representation of variational forms involving vertex integrals. This
allows for subsequent automated assembly of variational forms
involving single vertex integrals (as illustrated above) and also in
combination with other (cell, interior facet, exterior facet)
integrals.



 §.§ DOLFIN features for solving collections of ODE systems




  §.§.§ Assembly of vertex integrals


We have extended DOLFIN with support for automated assembly of variational forms that include vertex
integrals. The support is currently limited to forms
defined over test and trial spaces with vertex-based degrees of
freedom only. The assembly algorithm follows the standard finite
element assembly pattern by iterating over the vertices of the mesh,
computing the map from local to global degrees of freedom, evaluating
the local element tensor based on generated code,
and adding the contributions to the global tensor. The vertex integral
assembler runs natively in parallel via MPI.



  §.§.§ Specification and generation of multistage schemes


Since version , DOLFIN supports the specification of multistage
schemes of the
form (<ref>)–(<ref>) for the
solution of collections of ODE systems of the
form (<ref>), via their Butcher tableaus. The DOLFIN
class  takes as input the right hand
side expression f, the solution y, the Butcher tableau specified
via a, b and c, and secondary variables such as a time variable,
the (integer) order of the scheme. From this specification, DOLFIN automatically generates a
variational formulation of each of the separate stages in the
multistage scheme, with each stage in accordance
with (<ref>)-(<ref>). DOLFIN
can also automatically generate the variational formulations
corresponding to the adjoint scheme (<ref>) and
to the tangent linear scheme (<ref>) on
demand. A set of common multistage schemes are predefined including
Crank-Nicolson, explicit Euler, implicit Euler, 4^th-order
explicit Runge-Kutta, ESDIRK3, ESDIRK4 <cit.>.

Similar features have also been implemented for easy specification of
Rush-Larsen schemes cf. (<ref>)–(<ref>) and the
automated generation of the corresponding adjoint and tangent linear
schemes cf. e.g. (<ref>) and
(<ref>) through the class
. Both
 and 
subclass the  class.



  §.§.§ PointIntegralSolvers


To allow for the efficient solution of the multistage schemes for
collections of systems of ODEs, we have introduced a targeted
 class in DOLFIN . This solver
class takes a  as input and its main
functionality is to compute the solutions over a single time step.

The solver iterates over all vertices of the mesh and solves for the
relevant (stage and/or final) variables at each vertex. The solution
algorithm for each stage depends on whether the stage is explicit or
implicit. For implicit stages, a custom Newton solver is invoked that
allows for Jacobian reuse across stages, across vertices, and/or
across time steps on demand. As the resulting linear systems are
typically small and dense, a direct (LU) algorithm is used for the
inner solves. For explicit stages, a simple vector update is
performed.

The point integral solver runs natively in parallel via MPI. For a
cell-partitioned mesh distributed between N processes, each cell is
owned by one process and that process performs the solve for each
vertex on the cell (once). As little communication is required
between processes, the total solve is expected to scale
linearly in N.



 §.§ Extensions to dolfin-adjoint


Dolfin-adjoint has been extended to support the new features including
automatically deriving and computing adjoint and tangent linear
solutions for . The symbolic derivation of
the variational formulation for the adjoint and tangent linear
equations for the  are based
on (<ref>), (<ref>)
for the Butcher tableau defined schemes
and (<ref>), (<ref>) and its
analogies for the Rush-Larsen schemes.

The overall result is that operator splitting algorithms as described
in Section <ref> can be fully specified and
efficiently solved within FEniCS. Moreover, the corresponding adjoint
and tangent linear models and first and second order functional
derivatives may be efficiently computed using dolfin-adjoint with only
a few lines of additional code.



§ APPLICATIONS



In this section, we demonstrate the applicability and performance of
our implementation by considering the two examples presented in the
introduction, originating from the computational modelling of
mitochondrial swelling and cardiac electrophysiology respectively. The
complete supplementary code is openly available, see <cit.>.



 §.§ Application: mitochondrial swelling


We consider the mathematical model defined by (<ref>).
Inspired by <cit.>, we
let x = (x_0, x_1) ∈Ω = [0, 1]^2 and t ∈ [0, T] with T
= 35, and consider the initial conditions N_1, 0 = 1, N_2,
  0 = 0, N_3, 0 = 0 and

    u_0(x) = 30/∫_Ω M(x)  M(x),    M(x_0, x_1) =
      1/2 π e^- 0.5 (X_0^2 + X_1^2),

with X_0 = (α - β) x_0 + β, X_1 = (α - β) x_1
+ β taking α = 3 and β = -1. We consider

    f(s) =
      
        0     s < C^-, 
    
        f^∗    s > C^+, 
    f^∗/2 (1 -  cos (s - C^-/C^+ - C^-π )  )    otherwise,

and

    g(s) =
      
        g^∗    s > C^+, 
    g^∗/2 ( 1 - cos ( s/C^+π )  )     otherwise.

with C^- = 20, C^+ = 200, f^∗ = 1, g^∗ = 0.1,
d_1 = 2 × 10^-6, and d_2 = 30.

We are interested in the total amount of completely swollen
mitochondria (whose density is given by N_3) at the final time T and its
sensitivity to the initial condition for u (u_0). Thus, our
functional of interest J is given by

    J = ∫_Ω N_3(T)

and our goal is to compute the derivative of J with respect to u_0.

We use a second order Strang splitting scheme
for (<ref>), a Crank-Nicolson discretization
in time and continuous piecewise linear finite elements in space.
The scheme then
reads as: given initial conditions u_0, N_1,0, N_2,0,
N_3,0, and time points {t_0, t_1, …, t_N} with time step
κ_n = t_n+1 - t_n, then for each n = 0, 1, …, N - 1:


  * Compute u^∗ and N_i, 0^∗ i = 1, 2, 3 solving
  
    u = d_2 g(u) N_2,  N_1 = - f(u) N_1,  N_2 = f(u) N_1 - g(u) N_2,  N_3 = g(u) N_2,

  over Ω× [t_n, t_n + 1/2κ_n] with initial
  conditions u^n, N_i^n for i = 1, 2, 3.

  * Compute a solution u^† of: find u_h ∈ V_h such that
  
    ∫_Ω ( u_h - u^∗ ) v
         + κ_n  𝒜 (u_h)  · v  = 0

  for all v ∈ V_h with u = 1/2 (u + u^∗ ) and where 𝒜(u) = u or 𝒜(u) =
  |u|^q-2 u.

  * Compute solutions u^n+1 and N_i^n+1, i = 1, 2, 3,
  solving (<ref>) over Ω× [t_n +
    1/2κ_n, t_n+1] with initial conditions
  u^† and N_i^∗, i = 1, 2, 3.

We choose to discretize (<ref>) via the ESDIRK4
method in time and set a tolerance of 10^-10 for the inner Newton
solves. Further, we take q = 3, take κ_n = 0.5, and let
_h be a uniform tessellation of Ω with N_x × N_x
× 2 triangles.

The initial and final solutions are presented in
Figure <ref> while the L^2-gradient of the
objective functional J given by (<ref>) with
respect to the initial calcium concentration u_0 is presented in
Figure <ref>. 



To verify that the computed gradient is correct, we performed a Taylor
test at a given spatial and temporal resolution with a given
perturbation seed. The results are listed in
Table <ref> and demonstrate the expected orders of
convergence, indicating a correctly computed gradient.


In a multi-stage scheme, a number of intermediate stage solutions are
computed during the solution process. To avoid excessive memory usage,
dolfin-adjoint does not store these stage solutions. Thus, in order to
compute the adjoint solution, the stage solutions are recomputed for
each time-step. As a consequence, the minimal ratio of adjoint runtime 
to forward runtime for multistage ODE solves using this strategy
is approximately 2. For a
linear PDE solve, the optimal ratio of adjoint run time to forward run
time, assuming all linear solves take the same amount of time, is 1. For a nonlinear PDE
solve via a Newton iteration, the optimal ratio of adjoint run
time to forward run time, assuming all linear solves take the same
amount of time, is
1/N where N is the number of average Newton iterations in the
forward solve. These optimal ratios for linear and nonlinear PDE
solves have observed for typical dolfin-adjoint
usage <cit.>. For this example, two Newton iterations
were required on average to solve the PDEs to within a tolerance of
10^-10 both for N_x = 40 and 80, and assembly dominated the
PDE solve runtime. Thus for this example with two multistage ODE
solves and a nonlinear PDE solve in each timestep, the optimal ratio
of adjoint runtime to forward runtime is expected to be in the range
0.5-2, and closer to the upper bound than the lower bound depending
on the distribution of computational cost between the PDEs and ODEs.

Experimentally observed timings for this example are listed in
Tables <ref>–<ref> for N_x
= 40, 80, 160. From Table <ref>, we observe that
the adjoint-to-forward runtime ratio for the total solve is in the
range 1.38-1.95 and decreasing with increasing problem size. The
corresponding gradient-to-forward runtime is in the range 1.44 -
2.28, again decreasing with increasing problem size as expected. From
Table <ref>, we observe that adjoint-to-forward
runtime ratio for the ODE solves is in the range 2.35-2.41 for this
set of problem sizes, also decreasing with increasing problem size. We
also note that the ODE solve runtime, both for the forward and adjoint
solves, appears to scale linearly with the problem size M = N_x^2 as
anticipated and as is optimal. The ODE solves account for 37-40% of
the total forward runtime for this example.



We also conducted numerical experiments varying the multistage scheme
used in the ODE solves, including the 1-stage Backward Euler (BDF1),
the 2-stage Crank-Nicolson (CN2), the 3-stage ESDIRK3 scheme in
addition to the previously reported 4-stage ESDIRK4 scheme. The
results (using N_x = 40) are reported in
Table <ref>. We observe that the
adjoint-to-forward runtime ratio ranges from 2.04 (BDF1) to 2.41
(ESDIRK3, ESDIRK4) and the general trend is that this ratio increases
with the number of stages as expected, but stabilizes around
2.3-2.4.




 §.§ Application: cardiac electrophysiology (2D)



In this example, we consider the bidomain
equations (<ref>) over a two-dimensional rectangular
domain Ω = [0, 50] × [0, 50] (mm) with coordinates (x_0,
x_1). We will consider a set of different cardiac cell models of
increasing complexity: a reparametrized FitzHugh-Nagumo (FHN) model
with 1 ODE state variable <cit.>, the Beeler-Reuter (BR)
model with 7 ODE state variables <cit.>, the ten
Tusscher and Panfilov (TTP) epicardial cell model with 18 ODE state
variables <cit.>, and the Grandi et al (GPB)
cell model with 38 ODE state variables <cit.>. The
description of each of these cell models is available via the CellML
repository, and our implementation of the ionic current I_ ion
and F arising in (<ref>) was automatically generated
from the CellML models for the BR, TTP and GPB models. We refer to the
supplementary code <cit.> for the precise description of the
models including our choice of FitzHugh-Nagumo parameters.

Our parameter setup is otherwise as follows. We let χ = 140
(mm^-1), C_m = 0.01 (μF/mm^2), and let M_i =
diag (g_ if, g_ if) and M_e = diag
(g_ ef, g_ es) where g_ ef = 0.625/(χ C_m),
g_ es = 0.236/(χ C_m) and g_ if = 0.174/(χ
C_m). We let I_s = 0 and set the initial condition

    v_0 = 10  ( x_0/50 ) ^2 + 10

for v. The other state variables are initialized to the default
initial conditions given by the respective CellML models.

We use a second-order Strang splitting scheme (i.e. θ = 0.5
for (<ref>)–(<ref>)), a
Crank-Nicolson discretization in time and continuous piecewise linear
finite elements in space for the PDE system, and different Rush-Larsen
type discretizations (RL1, GRL1, RL2, GRL2) for the time
discretization of the ODE systems[We observed spurious wave
  propagation results when using a first-order splitting scheme and
  implicit Euler for with the BR model.]. The coupled linear systems
that arise at each PDE step (<ref>) were solved with
a block-preconditioned GMRES scheme with a relative tolerance of
10^-10. The complete details of the solver are detailed in the
supplementary code <cit.>.

We take κ_n = 0.1 (ms), and consider a uniform mesh of the
computational domain with N_x × N_x × 2 triangles. The
coarsest mesh considered used N_x = 160. For this case,
simulations converged using any of the cell models and schemes, and
gave qualitatively correct results, except the Grandi cell model
discretized by the RL1 scheme for which the numerical solution scheme
failed to converge.



We are interested in computing the gradient with respect to the
extracellular fiber conductivity g_ ef of the following
objective functional

    J(v) = ∑_i=1^5 ∫_Ω ( v(t_i)^2 + s(t_i)^2  )

for a equidistributed set of time points t_i. The model initial condition and solution for the
transmembrane potential at T = 100, together with the L^2-gradient
of the objective functional with respect to the fiber conductivity are
illustrated in Figure <ref>.

The following numerical experiments were performed on the Abel
supercomputer at the University of Oslo. All runtime experiments were
repeated three times, of which the minimum timing values are reported
here.



  §.§.§ Verifying the correctness of the discrete gradient



To verify that the computed functional gradient is correct, we
performed a Taylor test in a random perturbation direction at given
spatial and temporal resolutions. The results with the Fitz-Hugh
Nagumo cell model and the GRL1 scheme are listed in
Table <ref> and demonstrate the expected
orders of convergence without a computed gradient (order 1) and when
using the computed discrete gradient (order 2) in the Taylor
expansion.

We also performed similar experiments for the other cell models (BR,
TTP, and GPB) and other Rush-Larsen solution schemes (RL1, GRL2,
RL2). We used the same settings as in
Table <ref>, except for the Grandi cell model
with any other Rush-Larsen scheme than GRL1, where the time step had to
be reduced to 0.01 in order for the forward solver to converge. We
obtained the expected convergence order for all combination of cell
models and schemes tested. These results indicate that the
automatically derived and computed adjoints and gradient are correct
for all the Rush-Larsen schemes implemented.



  §.§.§ Adjoint runtime performance


As the bidomain PDEs (<ref>) are linear, a rough
theoretical estimate of the optimal PDE adjoint to PDE forward runtime
ratio is 1, assuming that the solution times in the adjoint and
forward PDE solves are equal and dominated by the finite element
assembly and linear solvers. Similarly, the ODEs are solved with
explicit Rush-Larsen schemes with a theoretically optimal
adjoint-to-forward ratio of 1 as well.

In Table <ref>, we list the total runtime, and
the runtime components corresponding to only the ODE solves, only the
PDE solves and an intermediate variable update (merge) step, for both
the forward and the adjoint solution computation for a test case with
T = 10.0, κ_n = 0.1, N_x = 160, the GRL1 scheme and the
four different cardiac cell models (FHN, BR, TTP, GPB).



From Table <ref>, for the forward solves, we
observe that the total runtime increases with cell model
complexity. The forward PDE step runtime is essentially independent of
the cell model, while the forward ODE step runtime increases with the
cell model complexity, as expected: the ODE runtime corresponds to
approximately 10 % of the PDE runtime for the simplest
FitzHugh-Nagumo model, and approximately 114 % of the PDE runtime
for the most complex Grandi model. The runtime of the forward merge
step is insignificant in comparison to the ODE and PDE steps, but
increases with cell model complexity.

We observe that also the adjoint PDE runtimes are comparable for all
cell models. The resulting adjoint-to-forward PDE step ratios are in
the range 0.79-0.82. The decrease in compute time for the adjoint
PDE solves compared to the forward PDE solves are attributable to the Krylov
solvers: fewer Krylov solver iterations were necessary for convergence
for the adjoint solves compared to the forward solves.

For the ODE solves, we observe that the adjoint-to-forward ratio for
the ODE solves is in the range 1.88-2.58 for the different cell
models. Additionally, we observe that the runtime of the adjoint merge step is
significantly larger than the corresponding time for the forward merge
step and with an increasing adjoint-to-forward ratio with increasing
cell model complexity 3.76-5.40. We suggest that the reason for the
increased and increasing adjoint runtime is that the adjoint merge
step involves additional assembly of variational forms of complexity
comparable to that of the cell model.

Overall, we observe that the adjoint-to-forward total runtime ratios
are in the range 1.15-1.66 for the different cell models considered,
with increasing adjoint-to-forward ratio with increasing cell model
complexity. Since the PDE runtime dominates the total runtime for the
simpler cell models, the total adjoint-to-forward ratio is closer to
1 for these models. On the other hand, for the more complicated cell
models, the ODE runtimes dominate, and so we observe
adjoint-to-forward ratios closer to 2 for these.

In the last column of Table <ref>, we observe
that the combined ODE and PDE solve time contributes to between 97%
(for FitzHugh-Nagumo) and 86% (for Grandi) of the total forward
runtime for the cell models considered. The remaining forward runtime
cost is primarily caused by the initialization routines such as
loading the computational mesh, and creating the required function
spaces. For the corresponding adjoint runtimes, the combined ODE and
PDE solve times contribute to between 86% (for FitzHugh-Nagumo) and
85% (for Grandi) of the total run time, and the percentage
contribution decreases with the complexity of the cell model.  The
remaining runtime cost can primarily be attributed to recording the
forward states during the adjoint solve (which is included in the
reported total run time but not in the ODE or PDE runtimes here). We
also timed the computation of the discrete gradient and noted that the
cost of additionally computing the gradient was negligible (less than
0.1% of the adjoint solve runtime).



  §.§.§ Forward and adjoint parallel scalability



In this section, we discuss the weak and strong parallel scalability
of the forward and adjoint ODE solvers and the scalability of the
point integral solvers in particular applied to this 2D
electrophysiology test case. We focus on the scalability of the point
integral solver particular as this is the primary new feature
described in this paper. Previous studies have discussed the parallel
scalability of FEniCS in general <cit.>.

We ran a series of weak scaling experiments on Abel, hosted by the
Norwegian national computing infrastructure (NOTUR), for the
previously considered series of cell models (FHN, BR, TTP, GPB) on an
increasing number of CPUs (1, 4, 16, 64) and with an increasing mesh
resolution.  We timed the total runtime of the point integral solves,
both for the forward solves and the adjoint solves. For each cell
model and resolution, we ran three experiments and extracted the
minimal run times. The results are listed in
Table <ref>. For both the forward and adjoint
run times, we computed the parallel efficiency (PE) as the ratio of
the 1-CPU runtime to the n-CPU run time for n = 4, 16, 64, also
listed in Table <ref>. The optimal PE would be
1. The last column of the Table lists the computed the
adjoint-to-forward ratio for these solves.

For the forward runtimes, we observe that the parallel efficiencies
range from 0.66 to 0.95 for all cell models and resolutions, with
efficiencies higher than 0.8 for all but the least complex cell
model. In general, we observe that the parallel efficiency increases
with increasing cell model complexity and that it decreases moderately
with the number of CPUs. Some of the loss of parallel efficiency (from
1) is attributable to only a near-perfect distribution of mesh
vertices across CPUs; we observed approximately 10% variation between
CPUs in the number of local vertices.

For the adjoint runtimes, we observe parallel efficiencies ranging
from 0.75 (FHN on 64 cores) to 0.95, again with efficiencies higher
than 0.8 for all but the least complex cell model. In general, we note
the same trends as for the forward run times: increasing parallel
efficiency with increasing cell model complexity and moderately
decreasing parallel efficiency with increasing number of CPUs. We
emphasize that we thus achieve the same parallel efficiency for
the adjoint point integral solves as for the forward point integral solves.

Comparing the forward and adjoint runtimes, we see that the
adjoint-to-forward ratio for the point integral solvers range from
1.19 to 1.56. Comparing these numbers with the adjoint-to-forward
ratios of the overall ODE solves as discussed in the previous, we
observe that the adjoint-to-forward point integral solver performance
is better. The adjoint-to-forward ratio is moderately increasing with
cell model complexity, but is stable with respect to the number of
CPUs. This latter point again illustrates that the adjoint point
integral solvers demonstrate the same parallel efficiency as the
forward point integral solves.


We also ran a series of strong scaling experiments on Abel for the
previously considered series of cell models (FHN, BR, TTP, GPB) on an
increasing number of CPUs (1, 4, 16, 64) and with a fixed mesh
resolution. We timed the total runtime of the point integral solver
steps for both the forward and adjoint solves. Again, for each cell
model and resolution, we ran three experiments and extracted the
minimal run times. The results are listed in
Table <ref>. For both the forward and adjoint run
times, we computed the (strong scaling) parallel efficiency (PE) as
the ratio of the 1-CPU runtime to n-times the n-CPU runtime for n
= 4, 16, 64, also listed in Table <ref>. The
optimal (strong scaling) PE would be 1. The last column of the Table
lists the computed the adjoint-to-forward ratio for these solves.

We observe that the parallel efficiencies are comparable to the
results from the weak scaling test case with efficiencies in the range
from 0.60 to 0.94. We note that the adjoint strong parallel efficiency
is comparable to the forward strong parallel efficiency. The
adjoint-to-forward ratios for this strong scaling test are also in the
same range as for the weak scaling test.





 §.§ Application: biventricular cardiac electrophysiology (3D)



In this example, we aim to compute the sensitivity of the squared
L^2-norm of the transmembrane potential with respect to its initial
condition, in order to illustrate the features discussed in this paper
for a more complex test case. This example also aims to illustrate how
different representations may seamlessly be used for the control
variables; this is useful for instance when computing sensitivities
with respect to spatially constant, regionally defined or highly
resolved control variables.

Inspired by <cit.>, we consider the monodomain
variation of (<ref>) over a mesh of a three-dimensional
bi-ventricular domain Ω with the ten Tusscher and Panfilov
epicardial cell model <cit.> as detailed in
Section <ref>. For the monodomain equations, we replace
the bidomain equations (<ref>)–(<ref>)
by: find v such that

    v_t -  M  v = - I_ ion(v, s)     in Ω.

with homogeneous Neumann boundary conditions and with v(t = 0) =
v_0. As before, we let χ = 140 (mm^-1), C_m = 0.01
(μF/mm^2), and for simplicity (ignoring realistic,
spatially-varying fiber directions), let M = diag (g_ f,
g_ s, g_ n) where:

    g_f = 0.255,    g_n = 0.0775,    g_s = 0.0775.


To solve the coupled equations (<ref>)
and (<ref>), we use a second-order Strang splitting scheme
(θ = 0.5),
a Crank-Nicolson discretization in time and continuous piecewise
linear finite elements in space for the resulting PDEs, and the
first-order generalized Rush-Larsen (GRL1) scheme for the resulting
ODEs. We take κ_n = 0.05 (ms). The moderately coarse mesh has
mesh cell diameters in the range [0.61, 2.70] (mm) and bounding box
[0.103, 40.2] × [-15.8, 36.8] × [-30, 13.9] (mm^3), with
45 625 vertices and 236 816 cells, and is illustrated in
Figure <ref>.

To quantify the transmembrane potential exceeding a zero threshold ,
we consider the following objective functional:

    J(v) = ∫_Ω_v v(T)^2 ,

at a fixed time T = 6.0 (ms) and where Ω_v = { x ∈Ω  |   v(x) ≥ 0 }. We are
interested in computing the gradient of J with respect to the initial
condition u_0.

We consider two different representations of the initial
condition to illustrate computing the gradient with respect to low and
high resolution fields. For both cases, we let

    v_0(x) = 
        -61.1     x_2 > z_ mid, 
    
        -61.3     x_2 ≤ z_ mid,

with z_ mid = -8.07, but consider (i) v_0 as represented by
two constants (c_0, c_1) i.e. the two-dimensional space ^2 and
(ii) v_0 represented by a continuous piecewise linear defined
relative to the mesh i.e. as a n_ verts-dimensional space V
where n_ verts denotes the number of mesh vertices. 


The initial condition and computed solution at t = 6.0 (ms) are
illustrated in Figure <ref>. (With reference
to the chosen perspective in this Figure, we denote the subdomain
where x_2 > z_ mid as the top part of the domain and the
subdomain x_2 ≤ z_ mid as the bottom part.)

For the case where the initial condition is represented at a low
resolution by v_0 = (c_0, c_1) ∈^2, we considered a series of
experiments starting with different initial conditions v_0^α =
(c_0 + α, c_1 + α) ∈^2 for α∈{ -0.1, -0.05,
0.0, 0.05, 0.1 }. The resulting gradients (in ^2) with respect
to v_0 ∈^2 are illustrated in
Figure <ref> (left). We observe that as we
increase α, the gradient with respect to the top initial
condition c_0 goes from large and positive to small and negative,
while the gradient with respect to the bottom initial condition c_1
varies from moderately negative to close-to-zero and rapidly to large
and positive. The rapid changes in the gradients illustrate the
nonlinear nature of the problem at hand.

The sensitivity of the objective functional J defined
by (<ref>) with respect to the initial condition v_0
∈ V is illustrated in Figure <ref>
(right). We observe that the gradient in the upper part of the domain
is large and negative, that the gradient in a boundary zone is large
and positive and the gradient in the lower part of the domain is close
to zero. This indicates that infinitesimal increases/decreases in the
upper part of the domain will lead to a large decrease/increase in the
functional value, infinitesimal increases/decreases in the boundary
zone domain will lead to a large increase/decrease in the functional
value, and that infinitesimal increases/decreases in the lower part of
the domain will have relatively small effect on the functional value.





§ CONCLUDING REMARKS



We have presented a high-level framework that allows for the
specification of coupled PDE-ODE systems and their efficient forward
and adjoint solution with operator splitting schemes. We have
illustrated the features of the framework with a series of examples
originating from cell modelling and computational cardiac
electrophysiology. Our numerical results indicate adjoint runtime
performance near optimal (measured in terms of adjoint-to-forward
runtimes), and parallel efficiency indices of 84-94% for the more
realistic cardiac cell models.

The main advantages of the framework are the speed of developing
solvers for new systems, and the ability to automatically derive their
adjoints from the high-level specification. This allows for flexible
exploration of models when the relevant governing equations are
uncertain, and for the identification of unknown parameters via the
solution of inverse problems as demonstrated
e.g. in <cit.>. Both of these capabilities are of
significant importance in numerous areas of scientific computing, and
in particular in computational medicine, physiology and
biology. Future work will focus on further improving the parallel
scalability of the solvers and its application to problems in
personalized medicine.


siamplain


