A p-adaptive, implicit-explicit mixed finite element method for reaction-diffusion problems
Mebratu Wakeni, Ankush Aggarwal, Lukasz Kaczmarczyk, Andrew McBride, Ignatios Athanasiadis, Chris Pearce, Paul Steinmann
AA p-adaptive, implicit-explicit mixed finite element method for reaction-di ff usionproblems Mebratu Wakeni a , Ankush Aggarwal a , Lukasz Kaczmarczyk a , Andrew McBride a , Ignatios Athanasiadis a , ChrisPearce a , Paul Steinmann a,b a Glasgow Computational Engineering Centre, University of Glasgow, Glasgow, G12 8QQ, United Kingdom b Institute of Applied Mechanics, Friedrich-Alexander University of Erlangen-Nuremberg, Paul-Gordan-Str. 3, D-91052, Erlangen, Germany
Abstract
A new class of implicit-explicit (IMEX) methods combined with a p-adaptive mixed finite element formulation isproposed to simulate the di ff usion of reacting species. Hierarchical polynomial functions are used to construct an H (div )-conforming base for the flux vectors, and a non-conforming L base for the mass concentration of the species.The mixed formulation captures the distinct nonlinearities associated with the constitutive flux equations and thereaction terms. The IMEX method conveniently treats these two sources of nonlinearity implicitly and explicitly,respectively, within a single time-stepping framework. The combination of the p-adaptive mixed formulation and theIMEX method delivers a robust and e ffi cient algorithm. The proposed methods eliminate the coupled e ff ect of meshsize and time step on the algorithmic stability. A residual based a posteriori error estimate that provides an upperbound of the natural error norm is derived. The availability of such estimate which can be obtained with minimalcomputational e ff ort and the hierarchical construction of the finite element spaces allow for the formulation of ane ffi cient p-adaptive algorithm. A series of numerical examples demonstrate the performance of the approach. It isshown that the method with the p-adaptive strategy accurately solves problems involving travelling waves, and thosewith discontinuities and singularities. The flexibility of the formulation is also illustrated via selected applications inpattern formation and electrophysiology. Keywords:
Implicit-explicit method; Mixed formulation; Hierarchical basis functions; H (div )-conforming;Reaction-di ff usion equation; p-adaptivity
1. Introduction
The spatio-temporal dynamics of multiple species interacting through a combination of two distinct mechanisms,namely reaction and di ff usion, can be described by reaction-di ff usion equations. Reaction refers to the inter / intraspecies interactions, resulting in the production and extinction of species. It is embodied in a term that is referred toas reaction kinetics f , a function of the mass concentration(s) m of the involved species. Di ff usion refers to the flowof substance (concentration) in space, and it is mathematically described by a flux h related to m (and / or its spatialgradient) through a constitutive equation. Reaction-di ff usion models are relevant in various important applications,including tissue morphogenesis and pattern formation [1–4], tissue remodelling [5–10], electrophysiology [11, 12],and epidemiology [13–15].The aforementioned applications motivate the need for robust and e ffi cient numerical methods for solving reaction-di ff usion problems. Various numerical methods have been proposed for approximating the solutions of reaction-di ff usion problems. Meshless methods in conjunction with operator-splitting techniques were used in [16] in one-and two-dimensions. In [17], Fitzhugh-Nagumo type models are solved using a multidomain algorithm based on apseudospectral approach. A review of some finite di ff erence based methods in one-dimension can be found in [18].Furthermore, a finite di ff erence scheme was constructed for the simulation of waves in excitable media using a two-variable reaction-di ff usion equation in [19]. However, finite di ff erence and spectral algorithms are suitable only forapproximations over relatively simple domains. Preprint submitted to arXiv January 5, 2021 a r X i v : . [ m a t h . NA ] J a n .2. Spatial discretisation The finite element method, due to its capabilities in handling arbitrary geometry and nonlinearities and its strongtheoretical foundation, is a natural choice for solving reaction-di ff usion problems. The majority of the finite elementnumerical approaches used in the literature are based on the standard, single-field formulation. The standard formu-lation for reaction-di ff usion was employed in a computational framework for the coupling of reaction-di ff usion andelasticity in [20]. In [21], a multilevel finite element approach with spatial and temporal adaptivity was constructedfor reaction-di ff usion problems. In [22], a projected finite element approach on stationary closed surface geometries,together with a backward Euler time integration, was used to simulate pattern-formation in biological applications. In[23], a moving mesh finite element method was constructed for simulating chemotaxis in two-dimensions. A multi-grid finite element method on stationary and evolving surfaces was proposed in [24]. A semi-linear multistep finiteelement was constructed in [25] for the two-dimensional simulation of pattern formation in ecological application. Insuch standard formulations, only the mass concentration m is solved for. The other physically important quantity, theflux h , is obtained as a post-processing step, decreasing the accuracy of its approximation.Solutions of reaction-di ff usion problems exhibit a variety of phenomena from the formation of travelling waves tocomplex structures like dissipative solitons. Some solutions may even involve low regularity features such as evolvingjump discontinuities and singularities. The H -conforming basis functions used in the standard formulation impose anunnecessarily high regularity requirement; for solutions displaying low regularity features the approximate solutionmay never converges to the true solution. Mixed finite element methods [see e.g. 26–28] o ff er an elegant solution forsuch problems. Mixed finite element methods are two-field formulations, which employ an H (div )-conforming basisfor the flux and a L -conforming basis for the mass concentration. This combination of basis functions relaxes theconformity requirements, allowing a wider class of solutions to be approximated accurately.Numerical studies of reaction-di ff usion type problems using mixed methods are, by comparison with standardformulations, relatively few. In [29], a stabilised mixed formulation in combination with a first-order implicit timeintegration was proposed for solving steady and unsteady state reaction-di ff usion problems. In this approach, H -conforming finite element spaces are used for m , and L -conforming spaces for h . With regard to the finite elementspace used for m , such a method has no particular advantage over the standard finite element method in terms ofaccuracy. A two-grid approach based on a variation of a mixed method with an implicit temporal integration wasproposed and analysed in [30]. Most numerical procedures for reaction-di ff usion equations that utilise finite elements, usually approach the tem-poral integration using either fully-implicit or fully-explicit methods. It is well-established that explicit methods canbe very e ffi cient and are easy to implement, however, they usually su ff er in terms of algorithmic stability, and im-pose severe time step restriction arising from the di ff usion term [31]. Implicit methods are known for their greaterstability, but can be challenging in terms of implementation, and are usually less e ffi cient as they lead to the solutionof a large system of algebraic equation. In addition, for nonlinear problems, it is necessary to derive and computetangent matrices that includes implicit nonlinearities at each time step, adding further ine ffi ciencies. Implicit-explicit(IMEX) methods mitigate such problems by combining the advantages of explicit and implicit methods [32]. Bytreating the non-local di ff usion term (involving a spatial derivative) implicitly, and the local reaction term (withouta spatial derivative) explicitly, one can eliminate the coupling e ff ect of spatial mesh size h and the time step size ∆ t on the stability condition. This allows the spatial mesh to be refined adaptively without the need for reducing thetime-step size. The application of IMEX methods appears conducive for reaction-di ff usion problems, however, mostof the literature on IMEX methods for such problems are limited to classical spatial discretisation techniques, such asfinite-di ff erence [31, 33, 34] and standard finite element [25, 35, 36].The mixed method allows for the nonlinearities that may appear in the flux constitutive equation and the reactionterm to be considered separately. For stability reasons, the flux constitutive equation must be treated implicitly, whilethe nonlinearities in the reaction term can be handled explicitly. In this presentation we present a robust and e ffi cientnumerical algorithm based on mixed formulation with IMEX temporal integration methods for problems of reaction-di ff usion type. 2 igure 1: Schematics of the domain Ω with its boundary partitions Γ E and Γ N The contribution is organised as follows. In Section 2, a general mathematical model of multi-species reaction-di ff usion systems is presented briefly. The weak formulation of the model using a mixed approach, is describedin Section 3. In Section 4, relevant aspects of the numerical procedure for the temporal discretisation using theIMEX method and the spatial approximation using mixed Galerkin approaches are presented. Finally, in Section 7the performance and capabilities of the proposed formulations are demonstrated using various numerical examples.Here, the performance of the mixed and the standard formulations are compared, and finally some selected examplesrelevant to pattern formation, ecology, and electrophysiology are simulated using the mixed method.
2. Model overview
Consider n species, each with mass concentration m i , where i = , , , . . . , n , interacting in an open, boundedregion Ω ⊂ IR d ( d = , , or 3). The local form of the mass balance, for each of the species, is given by˙ m i + div h i = f i ( m , . . . , m n ) , i = , . . . , n , (1)where h i denotes the concentration flux of the i th species, and f i is the chemical kinetics term that represents the rateof production or degradation of species concentration of the i th species as a result of its interaction with other species.In addition to the mass balance equation (1), a constitutive relation relating the flux h i to the mass concentration m i isrequired. A commonly used constitutive relation is given by h i = − D i ∇ m i , i = , . . . , n . (2)Here D i is a symmetric and positive-definite second-order tensor representing a potentially spatially varying di ff u-sivity / mobility of the i th species on the domain Ω . Let Γ N and Γ E be nonoverlapping portions of the boundary of Ω ,denoted by Γ (see Fig. 1), such that Γ N ∪ Γ E = Γ . The prescribed boundary conditions imposed on these partitions are m i = ¯ m i on Γ N , and (3) − h i · n = ¯ h i on Γ E , (4)3here n represents the unit outward normal vector to the boundary Γ . A complete description of the problem alsorequires the prescription of initial conditions for each m i , which read as m i ( x ) = m i ( x ) , at t = , ∀ x ∈ Ω .
3. Weak formulations
The focus here is on the mixed formulation. However, for the sake of completeness, the standard, single-fieldformulation is first briefly stated. Thereafter, a detailed presentation of the mixed formulation and its spatial andtemporal discretisation is given.In the context of initial-boundary value problems, such as analysed here, it is helpful to view functions of spaceand time as mappings from the time interval of interest I = [0 , T ] to the corresponding functional space. For example,a function u ∈ L ( Ω ; I ) is understood as the map u : I → L ( Ω ) (the space L ( Ω ) denotes the space of measurablefunctions which are square integrable over the domain Ω ). In addition to the functional space L ( Ω ; I ) we also makeuse of the space H ( Ω ; I ) and H (div , Ω ; I ), where H ( Ω ) = { m ∈ L ( Ω ) : ∇ m ∈ [ L ( Ω )] d } , and H (div , Ω ) = { h ∈ [ L ( Ω )] d : div h ∈ L ( Ω ) } . The natural norms endowed by H ( Ω ) and H (div , Ω ) are, respectively, given by (cid:107) m (cid:107) , Ω (cid:66) (cid:107) m (cid:107) , Ω + (cid:107)∇ m (cid:107) , Ω , and (cid:107) h (cid:107) , Ω (cid:66) (cid:107) h (cid:107) , Ω + (cid:107) div h (cid:107) , Ω , where (cid:107) · (cid:107) , Ω denotes the standard L -norm for scalar or vector-valued functions.The standard weak problem is defined as: Standard formulation find m i ∈ H ( Ω ), satisfying the boundary conditions (3), such thatdd t ( v , m i ) Ω + ( ∇ v , D i ∇ m i ) Ω = ( v , ¯ h i ) Γ E + (cid:96) i ( v ) , ∀ v ∈ H N ( Ω ) , (5)where H N ( Ω ) is a subspace of H ( Ω ) that contains functions whose trace on Γ N vanish. Here, it should be notedthat the test function v is time-independent. The functional (cid:96) i : H ( Ω ) → IR is defined by (cid:96) i ( v ) (cid:66) ( f i , v ) Ω . (6)The pairings ( · , · ) Ω and ( · , · ) Γ N / E represent the standard L inner product over the domain Ω and the boundary Γ N or Γ E , respectively.As can be seen from equation (5), the boundary condition (4) is incorporated into the weak form, while equation(3) is enforced as a constraint on the trial solutions. Thus, for the standard formulation, (4) is a natural boundarycondition and (3) is an essential boundary condition.For a mixed finite element formulation, in addition to m i the flux h i is an unknown variable. Thus, the constitutiverelation (2) is re-written as D − i h i + ∇ m i = . (7)The mixed weak form associated with equations (7) and (1) reads:4 ixed formulation find ( h i , m i ) ∈ H E (div , Ω ; I ) × L ( Ω ; I ) such that a i ( τ , h i ) − b ( τ , m i ) = ( τ · n , ¯ m i ) Γ N , ∀ τ ∈ H E (div , Ω ) (8)dd t c ( v , m i ) + b ( h , v ) = (cid:96) i ( v ) , ∀ v ∈ L ( Ω ) (9)where the spaces H E (div , Ω ) and H E (div , Ω ) are subspaces of H (div , Ω ). Also, the test functions v and t au are time-independent. The vector-valued functions in the former space satisfy the boundary condition (4), whereasfunctions in the latter space satisfy a vanishing normal component at the boundary Γ E .The bilinear forms a i : H (div , Ω ) × H (div , Ω ) → IR, b : H (div , Ω ) × L ( Ω ) → IR, and c : L ( Ω ) × L ( Ω ) → IR aredefined by a i ( τ , h i ) (cid:66) ( τ , D − i h i ) Ω , b ( τ , m i ) (cid:66) (div τ , m i ) Ω , c ( m i , v ) (cid:66) ( m i , v ) Ω , for any τ , h i ∈ H (div , Ω ) and m i , v ∈ L ( Ω ). Remarks.
1. When the L ( Ω ) space of test functions is replaced by its discrete counterpart, the test function v is chosen such thatit vanishes everywhere in the domain Ω except on a given element Ω e . This in turn implies that for v in (9) set to unityon the element Ω e , one obtains dd t (cid:90) Ω e m i d Ω + (cid:90) Ω e div h i d Ω = (cid:90) Ω e f i d Ω . As a result, such an approximation method is said to have a locally conservative property. That is, the conservationof mass (1) is satisfied on each element.2. The classification of boundary conditions in the mixed formulation is opposite to the standard single-field case. Inthe standard formulation, the boundary condition (4) is a natural one as it does not require a priori prescription on thespace of trial or test spaces. By contrast, it becomes an essential boundary condition in the mixed formulation sincethe trial and test functions require the normal flux at the boundary Γ E to be prescribed a priori. The role of equation (3) is also reversed, that is, it becomes essential in the standard formulation but natural in the mixed formulation.
4. Discretisation
The temporal discretisation of the weak formulations, (8) and (9), using a combination of implicit and explicitmethods is now presented. Then, the discrete counterparts of the spaces H (div , Ω ) and L ( Ω ) are detailed in thecontext of the hierarchical construction of shape functions over a triangular / tetrahedral mesh. Consider first the temporal discretisation of the mixed formulation using a class of IMEX methods. The timeinterval of interest is partitioned into subintervals [ t n − , t n ] with step size ∆ t n = t n − t n − . Note that the partition neednot be uniform, that is, the step-sizes need not be equal. As equation (8) is without a time derivative, we treat it fullyimplicitly at the current time t n . For the IMEX method only the second equation (9) involving time derivatives isrelevant. For clarity of notation, we drop the subscript i from the weak forms (8) and (9), and, in a general multistepcontext, replace each term by interpolation or extrapolation formulas as linear combinations of previous discrete5alues, as defined by dd t c ( • , m ) ≈ ◦ c ( • ) (cid:66) α r ∆ t n c ( • , m n ) + r − (cid:88) j = α j ∆ t n c ( • , m n + j − r ) , (10) b ( h , • ) ≈ (cid:98) b ( • ) (cid:66) β r b ( • , h n ) + r − (cid:88) j = β j b ( h n + j − r , • ) , (11) (cid:96) ( • ) ≈ ˜ (cid:96) ( • ) (cid:66) r − (cid:88) j = γ j (cid:96) n + j − r ( • ) , (12)where (cid:96) k ( • ), k = n − r , . . . , n −
1, represents the family of functionals (cid:96) ( • ) defined using the time discrete values of thereaction kinetics f k = f ( m k ), i.e., with reference to (6), (cid:96) k ( v ) = ( f k , v ) Ω . The coe ffi cients β , β , . . . , β r and α , α , . . . , α r correspond to the implicit interpolation formula (corresponding toequations (10) and (11)) for the value and its time derivative of a field at t n on the time interval [ t n − r , t n ]. γ , γ , . . . , γ r − are coe ffi cients of the explicit extrapolation (corresponding to the equation (12)) of a field at t n on the interval. Theinteger r represents the extent to which previous step solutions, starting from the current step, are included in thescheme. Remark 1.
Some of the commonly used IMEX schemes in the literature are: • IMEX schemes based on the Backward Di ff erentiation Formula (BDF)Second-order α = / , α = − , α = / ,β = , β = , β = ,γ = − , γ = , Third-order α = / , α = − / , α = − / , α = / ,β = / , β = − / , β = / , β = / ,γ = / , γ = − / , γ = / . • The second-order Crank-Nicholson – Adams-Bashforth scheme α = , α = − , α = ,β = , β = / , β = / ,γ = − / , γ = / . • The second-order additive Runge-Kutta scheme [37, 38] α = − , α = , α = β = , β = , β = γ = , γ = , Substituting the discrete approximations (10)-(12) into the weak formulation (9), together with the discrete equa-tion corresponding to equation (8) at the current time-step renders a n ( τ ) − b n ( τ ) = ( τ · n , ¯ m ) Γ N , ∀ τ ∈ H (div , Ω ) , and (13) ◦ c ( v ) + (cid:98) b ( v ) = ˜ (cid:96) ( v ) , ∀ v ∈ L ( Ω ) , (14)where a n ( τ ) (cid:66) a ( τ , h n ) , b n ( τ ) (cid:66) b ( τ , m n ) . Note, equation (13) and (14) constitute a boundary value problem at the time-step t n . IMEX methods can be viewedas multistep schemes involving r − r -order convergent in time.6 .2. Spatial discretisation Assume a regular decomposition T h of Ω into simplexes (triangles in 2D and tetrahedral in 3D). For a given k ∈ Z + (a non-negative integer), denote the set of all polynomials, on a given T ∈ T h , whose order is less than or equal to k by P k ( T ). For finite element methods which typically involve the use of non-uniform higher-order approximationson unstructured meshes, increasing the order of the polynomial space locally via p - and hp - adaptivity can lead tocomplications in enforcing global conformity of shape functions [39]. Hierarchical shape functions address suchproblems, as well as naturally supporting the use of p - and hp -adaptivity [40, 41]. The construction of hierarchicalshape functions of arbitrary order with various conformity conditions to obtain finite element subspaces for L ( Ω ), H ( Ω ), H (div , Ω ) on a general unstructured meshes is detailed in Ainsworth and Coyle [39]. An alternative con-struction of H-div conforming exact sequence element with arbitrary order has been proposed by Fuentes et al. [42].Note, however, that space in [39] consist of divergence-free zero normal functions was used in following numericalexamples in Section 7.Here, hierarchic shape functions are used to define the finite element spaces corresponding to the triangulation T h such that the test spaces, for concentration m and flux h , are defined as S h = { v h ∈ L ( Ω ) : v h | T ∈ P k ( T ) , where T ∈ T h } , and (15) V h = { τ h ∈ H (div , Ω ) : τ h | T ∈ [ P k + ( T )] dim and τ h · n = Γ E } , (16)where dim = , , or 3 refers to the spatial dimension. While the trial space for concentration m is also S h , the trialspace V h for the flux h is given by V h = { τ h ∈ H (div , Ω ) : τ h | T ∈ [ P k + ( T )] dim and − τ h · n = ¯ h on Γ E } . (17)It should be noted that to obtain a stable pair ( m h , h h ) the order of approximation for V h is required to be at least oneorder higher than that of S h , see, for example [26].For the approximation of ( h , m ), we employ the finite dimensional trial space ( V h , S h ) and test space ( V h , S h ),as defined in equations (15)-(17), following a Galerkin approach. Having specified the corresponding finite elementspaces, the spatio-temporal discrete form of the equations (13) and (14) assumes a block matrix system given by (cid:34) K BB T − σ M (cid:35) (cid:34) H n m n (cid:35) = (cid:34) FG (cid:35) , (18)where K IJ = a ( τ hI , τ hJ ) , B IL = − b ( τ hI , v hL ) , and M KL = c ( v hK , v hL ) . Here, I , J denote the global indices corresponding to the numbering of the basis elements of V h , while K , L corre-spond to that of S h . The right hand side F and G are given by F I = ( τ I · n , ¯ m h ) Γ E , and G K = β r r − (cid:88) j = [ γ j (cid:96) n + j − r ( v K ) − β j b ( v K , h n + j − rh )] , respectively . The vector H n and m n are the solution vectors containing the degrees of freedom (dof) associated with the currentvalues of h n and m n , respectively. The coe ffi cient σ = α r / [ ∆ t n β r ] is a shift coe ffi cient of the mass matrix M . Thematrices K and M are positive definite and symmetric. Solvability of the block system (18) also requires that B as a linear map is surjective (see, for example, [26] Section 3.3). The requirement that the order of the flux shapefunction should be at least one order higher than the mass concentration shape function is a su ffi cient condition forthe surjectivity of B .Formally, for su ffi ciently smooth solutions, the expected rate of convergence for the spatial approximation em-ploying the finite element spaces S h and V h will be of order k + L and the natural norms. Remarks. . One of the most important implications of the mixed formulation, from a computational perspective, is that thematrix M can be inverted locally on an element-by-element basis, and the inverse is sparse. This is due to the factthat there is no conformity requirement on m h ∈ S h ⊂ L ( Ω ) over element boundaries.2. The consequence of the above observation is that one can e ffi ciently solve the block system using a solver thatutilises a Schur complement preconditioner. More precisely, one can exactly compute the sparse Schur complement S = K + BM − B T in an e ffi cient manner.
5. A posteriori error estimators and p-adaptivity
Adaptive finite element methods are a fundamental numerical approach in science and engineering applications.The success of an adaptive algorithm relies on the availability of a good error indicator (or a posteriori error estimator)that provides an upper bound to the true approximation error, and the complexity of its implementation and computa-tion. In this section, we present a residual based error estimate that can be computed cheaply. This estimate togetherwith the hierarchical construction of the shape functions makes the method well suited for local p-adaptivity. Hierar-chical shape functions (and dofs) are associated with mesh entities such as vertices, edges, faces, and volumes, ratherthan nodes. For example, in 2D, if the local order of two adjacent faces F and F sharing an edge E are di ff erent,say order- k and order- k , respectively, then to satisfy the global H (div , Ω )-conformity one only needs to add localshape functions of order-max( k , k ) on the edge E . By contrast, since there is no continuity requirement for L ( Ω )along element interfaces, the shape functions are only associated with faces in the 2D case. Thus, the polynomialorder of L shape functions can be set independently in each element, thereby greatly simplifying the implementationof p-adaptivity.To underpin an e ff ective local p-adaptivity scheme, one requires a reliable a posteriori error estimate that providesan upper bound to the true error and, at best, decays with the same rate as the true error as the polynomial orderincreases uniformly. The energy norm on H (div , Ω ) × L ( Ω ), defined by (cid:107) h (cid:107) div , Ω + (cid:107) m (cid:107) , Ω (19)is the appropriate norm for measuring the magnitude of approximation errors in the mixed formulation (8) and (9).To develop the error estimator, and for the sake of simplicity, we shall consider a one-species mixed transientproblem with homogeneous boundary condition on m h , q , where the additional subscript used here is to denote theorder of the polynomial space, that is a ( τ h , k , h nh , k ) − b ( τ h , k , m nh , k − ) = , ∀ τ h , k ∈ V h , k (20) σ c ( v h , k − , m nh , k − ) + b ( h nh , k , v h , k − ) = ( g , v h , k − ) Ω , ∀ v h , k − ∈ S h , k − (21)where g = ˜ f − σ m n − h , k − , where ˜ f denotes an extrapolation of f at the previous time-step values of m h that is determinedby the specific type of the IMEX scheme, the remaining term in g is a contribution from temporal discretisation of˙ m h , k − = σ [ m nh , k − − m n − h , k ]. Recall that for stability the orders used in equations (20) and (21) for the flux and mass are k and k −
1, respectively ( k > h nh , k and m nh , k − are denoted as h k and m k − . The a posteriori error estimate corresponding to the energy norm (19) of the error is derived from the followingresidual and jump based errors. Consider first a su ffi ciently refined regular mesh T h of Ω , where h is the meshparameter. • Element residual errors corresponding to the constitutive and conservation of mass equations: Let K be anelement in T h , define η K , R , : = (cid:107) D ∇ m k − − h k (cid:107) , K , (22) η K , R , : = (cid:107) σ m k − + div h k − g (cid:107) , K . (23)8 Inter-element interface jump error: Let e be an edge shared by two adjacent elements K + and K − , then the jumperror is defined by η e , J : = h − / (cid:107) [[ m k − ]] (cid:107) , e , (24)where the jump operator is [[ m k − ]] = m + k − − m − k − , m + k − and m − k − are values of m k − at the edge e from the sideof K + and K − , respectively.The local error estimate over a given element K ∈ T h is thus defined by summing the error contributions from (22),(23) and (24), η K : = (cid:20) η K , R , + η K , R , + (cid:88) e ∈ ∂ K η e , J (cid:21) / , (25)and the global error estimate is given by η T h : = (cid:20) (cid:88) K ∈T h (cid:2) η K , R , + η K , R , (cid:3) + (cid:88) e ∈ Γ h η e , J (cid:21) / (26) The global residual error (26) provides an upper bound of the energy norm of the error. To show this, we makeuse of the following standard estimates1. Optimality estimate: Let ( h , m ) ∈ H (div , Ω ) × L ( Ω ) be the exact solution at the current time-step, i.e., t n , thenthere exists C >
0, such that (cid:107) h − h k (cid:107) div , Ω + (cid:107) m − m k − (cid:107) , Ω ≤ C (cid:26) inf τ k ∈V k (cid:107) h − τ k (cid:107) div , Ω + inf v k − ∈S k − (cid:107) m − v k − (cid:107) , Ω (cid:27) (27)2. Energy norm error estimate of the finite element solution: (cid:107) h k (cid:107) div , Ω + (cid:107) m k − (cid:107) , Ω ≤ C (cid:20) sup τ k ∈V k (cid:107) τ k (cid:107) div , Ω = (cid:8) a ( τ k , h k ) − b ( τ k , m k − ) (cid:9) + sup v k − ∈S k − (cid:107) v k − (cid:107) , Ω = (cid:8) σ c ( v k − , m k − ) + b ( h k , v k − ) } , (cid:21) (28)for some C > Saturation assumption : One of the most crucial ingredients towards the proof of an upper bound is the saturationassumption. Roughly, it states that the error norm decreases uniformly as we increase the order of approximationby one. More precisely, let ( h k + , m k ) and ( h k , m k − ) be approximate solutions of (20) and (21), then there is0 < β <
1, such that (cid:107) h − h k + (cid:107) div , Ω + (cid:107) m − m k (cid:107) , Ω < β (cid:2) (cid:107) h − h k (cid:107) div , Ω + (cid:107) m − m k − (cid:107) , Ω (cid:3) . (29)One can construct a mesh T h on which such a saturation estimate does not hold, however, for su ffi ciently refinedregular mesh it always hold true.Having the above results for the upper bound, it is su ffi cient to show that the error between successive approximationsis bounded from above. That is, there is a constant C > (cid:107) h k + − h k (cid:107) div , Ω + (cid:107) m k − m k − (cid:107) , Ω ≤ C η T h . (30)To show this, we first note that a ( τ k + , h k + − h k ) − b ( τ k + , m k − m k − ) = , ∀ τ k + ∈ V k + (31) σ c ( v k , m k − m k − ) + b ( h k + − h k , v k ) = . ∀ v k ∈ S k (32)9ence, by the estimate (28), we have, for some C > (cid:107) h k + − h k (cid:107) div , Ω + (cid:107) m k − m k − (cid:107) , Ω ≤ C (cid:20) sup τ k + ∈V k + (cid:107) τ k + (cid:107) div , Ω = (cid:8) a ( τ k + , h k + − h k ) − b ( τ k + , m k − m k − ) (cid:9) + sup v k ∈S k (cid:107) v k (cid:107) , Ω = (cid:8) σ c ( v k , m k − m k − ) + b ( h k + − h k , v k ) (cid:9)(cid:21) . (33)Now, since a ( τ k , h k + − h k ) − b ( τ k , m k − m k − ) = , ∀ τ k ∈ V k , (34)it follows for every τ k that a ( τ k + , h k + − h k ) − b ( τ k + , m k − m k − ) = a ( τ k + − τ k , h k + − h k ) − b ( τ k + − τ k , m k − m k − ) = − a ( τ k + − τ k , h k ) + b ( τ k + − τ k , m k − ) = (cid:88) K (cid:26) ( τ k + − τ k , h k ) K + (div ( τ k + − τ k ) , m k − ) K (cid:27) = (cid:88) K (cid:26) − ( τ k + − τ k , h k ) K + ( τ k + − τ k , ∇ m k − ) K (cid:27) + (cid:88) e ∈ Γ h ([ τ k + − τ k ] · n , [[ m k − ]]) e = (cid:88) K ( τ k + − τ k , h k − ∇ m k − ) K + (cid:88) e ∈ Γ h ([ τ k + − τ k ] · n , [[ m k − ]]) e ≤ (cid:88) K (cid:107) τ k + − τ k (cid:107) , K (cid:107) h k − ∇ m k − (cid:107) , K + (cid:88) e ∈ Γ h (cid:107) [ τ k + − τ k ] · n (cid:107) , e (cid:107) [[ m k − ]] (cid:107) , e = (cid:88) K η K , R , (cid:107) τ k + − τ k (cid:107) , K + (cid:88) e ∈ Γ h η e , J h / (cid:107) [ τ k + − τ k ] · n (cid:107) , e . Hence, we obtain thatsup τ k + ∈V k + (cid:107) τ k + (cid:107) div , Ω = (cid:8) a ( τ k + , h k + − h k ) − b ( τ k + , m k − m k − ) (cid:9) ≤ C (cid:26) (cid:88) K η K , R , + (cid:88) e η e , J (cid:27) / (35)for some constant C >
0. A similar argument also leads tosup v k ∈S k (cid:107)(cid:107) , Ω = (cid:8) σ c ( v k , m k − m k − ) + b ( h k + − h k , v k ) (cid:9) ≤ C (cid:26) (cid:88) K η K , R , (cid:27) / . (36)Therefore, for C = max( C , C ) we obtain the estimate (30). Employing the saturation estimate (29) and (30), it thenfollows that (cid:107) h − h k (cid:107) div , Ω + (cid:107) m − m k − (cid:107) , Ω ≤ C − β η T h . (37)Here the constant C depends only on the approximation order k and h .
6. Adaptive p-refinement strategy
Once the problem is solved with a given distribution of polynomial orders over the mesh entities, and the local aposteriori error η K over each element K ∈ T h is calculated, the next step is to apply a p-refinement strategy inspiredby the well-known bulk-chasing D¨orfler’s criterion [43]. The refinement algorithm is characterised by two parameters θ min and θ max , where 0 ≤ θ min < θ min ≤
1, and is performed in three stages:10 tage 1.
Given a posteriori error estimate η K on each element K ∈ T h and η MAX = max K ∈T h { η K } , the polynomialorder over element K is raised by one if η K ≥ θ max η MAX , or reduced by one if η K ≤ θ min η MAX . After applying this first stage of the adaptive process it may happen that the polynomial order distributionover adjacent elements be greater than one order. Numerical experiments (not presented here) reviled thathetrogeneity of polynomial order distribution results in undesirable oscillatory feature of the approximatedsolution. Hence, following this step, certain smoothing of polynomial order over the mesh is required, whichleads us to the next stage.
Stage 2.
To smooth the polynomial order distribution, we force the di ff erence in polynomial order between twoadjacent elements K and K (cid:48) to not exceed one, by resetting the order on the element with smaller degree to thatof with the higher degree minus one. That is, suppose order( K ) + < order( K (cid:48) ), then we reset order( K ) (cid:66) order( K ) − Stage 3.
This stage is responsible to maintain the H (div )-conformity of the space of flux functions after weexecute the above two stages. For each interface entity E shared by two elements K and K (cid:48) , we set the order asthe maximum of the polynomial orders over K and K (cid:48) .The adaptive p-refinement algorithm consisting of the above three stages is summarised in Algorithm 1. Followingthe above p-adaptive stages, one also needs to adjust the quadrature rules over the mesh entities appropriately in orderto match the polynomial order distributions optimally. Algorithm 1:
Adaptive p-refinements algorithm
Input: η K on each K ∈ T h , , θ max and θ min Stage 1.
Setting order on each K ; for K ∈ T h doif η K ≥ θ max η MAX then raise polynomial order on K by one; endif η K ≤ θ min η MAX then decrease polynomial order on K by one but not less than a minimum order, say 1; endendStage 2. Order smoothing; for E ∈ Γ h shared by two elements K and K (cid:48) in T h such that the order in K is greater than that of K (cid:48) by morethan 1 do set: order in K (cid:48) equals order in K minus 1; endStage 3. Setting order on the interfaces; for E ∈ Γ h doif E is on the boundary with only one adjecent element K then set order on E to be equal to that of on K ; else find adjacent elements K and K (cid:48) sharing E ;Set order on E to be maximum of orders on K and K (cid:48) ; endend . Numerical examples Two groups of numerical examples are presented: the first compares the convergence of results of the mixedscheme and the standard single-field formulation in approximating important aspects of the solution. These includesolutions involving singularities, and computation of the speed of travelling wave solutions. Suitability of the p-adaptive mixed formulation in terms of the features of the solution is also investigated. The second group showcasesthe capabilities of the p-adaptive, IMEX mixed formulation in simulating problems of practical importance: patternformation and electrophysiology. We investigated several IMEX schemes, for the examples presented in this sectionwe opted for the second-order additive Runge-Kutta scheme.The computer implementation of the proposed numerical scheme is carried out using the open-source libraryMOFEM [44]. The library integrates and utilises other open-source libraries such as MOAB, a mesh-oriented database[45, 46], and PETSc [47, 48]. The MOAB library is used to store and manage mesh related data, while PETSc is usedfor parallel operations involving linear algebra.The IMEX methods presented in Section 4 are implemented using the PETSc (Portable, Extensible Toolkit forScientific computations [47, 48]) time solvers [49].
Two cases are considered. In the first, a spatially smooth solution is considered with a piece-wise temporal profilethat stabilises after some specified time. Thus, the approximation error after a su ffi ciently long simulation time isassociated entirely with the spatial discretisation. The second case considers the approximation of a one-speciesFisher’s type problem on a square domain Ω with heterogeneous di ff usivity. a) Smooth manufactured solution It is well-known that both standard and the mixed finite element formulations are optimal in terms of convergencein the L -norm, i.e., O ( h p + ), for su ffi ciently smooth solutions, where p is the order of the finite element space. Notingthat the mesh size parameter h is inversely proportional to the number of degrees-of-freedom to the power dim , where dim = , g ( x , y ) = + sin(2 π x ) sin(2 π y ) , for ( x , y ) ∈ Ω . (38)Consider first a one-species reaction-di ff usion system over the domain Ω = [ − , (so that g vanishes on theboundary) with isotropic di ff usivity D = d I , d =
1, and then assume the exact (manufactured) solution for the massconcentration m m ( x , y , t ) = t g ( x , y ) , t < t ∗ , g ( x , y ) , t ≥ t ∗ . (39)for some given t ∗ . The right-hand-side source term f is given by the residual of the exact solution, i.e., f (cid:66) ˙ m + div ( D ∇ m ) . Note that ˙ m and div ( D ∇ m ) are also piecewise in time, i.e.,˙ m = g ( x , y ) , t < t ∗ , , t > t ∗ , and div ( D ∇ m ) = t div ( D ∇ g ) , t < t ∗ , div ( D ∇ g ) , t ≥ t ∗ . Consequently, the source term is also temporally piece-wise which stabilises to a time-independent profile after t ∗ . Note also that for t < t ∗ , ¯ m = t on the boundary Γ , and for t > t ∗ , ¯ m =
1. With this set up, the temporaldiscretisation error after a time t su ffi ciently greater than t ∗ will be negligible, and the total error is dominated bythe spatial approximation. In other words, it amounts to the approximation of the steady state case ( ˙ m =
0) with themanufactured solution m = g .The convergence results presented in Fig. 3 and 4 are obtained by a successive refinement of an initial uniformmesh with h = /
5, and the inflection time t ∗ and a uniform time step length ∆ t are chosen to be 1 and 0 .
1, respectively,12 igure 2: Exact (manufactured) solution for m at t = corresponding to each mesh. The simulations are run up to t =
10; a su ffi ciently long time to ensure that the temporaldiscretisation error is negligible.With the same construction of the manufactured solution (39), the convergence rate of the mixed formulation withrespect to the H -norm, given by (cid:2) (cid:107) m − m h (cid:107) , Ω + (cid:107) h − h h (cid:107) , Ω (cid:3) / , is expected to be one order higher than that of the standard formulation, as demonstrated in Fig. 4. This is due to thefact that the flux h h for the standard formulation is obtained by postprocessing from m h , unlike the mixed formulation,wherein the flux is directly approximated as a primary field variable.The next set of examples in this group aims at demonstrating the e ff ectiveness of the p-adaptive mixed formulationin resolving fine features of solutions e ffi ciently. For smooth solutions such as (38), the variability of the solution isalmost uniform on the larger scale. In this case, the application of p-adaptivity is less e ff ective since the error isdistributed almost uniformly. This is demonstrated in the convergence result displayed in Fig. 5. It shows that thep-adaptive strategy with parameters θ min = .
02 and θ max = .
8, representing a quite conservative adaptive strategy,produces a convergence trend which is not generally better than that of the uniform p-refinement. As expected, ateach adaptive step, as shown in Fig. 6, the error is distributed almost uniformly, which leads to the marking of mostof the elements for refinement. This corresponds to the convergence result shown in Fig. 6, which is not better thanthe uniform p-adaptive strategy.In contrast, when the solution is characterised by the presence of sudden spatial changes over the domain, suchas travelling waves, the p-adaptive algorithm becomes most e ff ective. To demonstrate this, we consider a smoothanalytical solution (manufactured) replacing the g in equation (38) by the bump function over the square domain Ω = [ − , , g ( x , y ) = exp (cid:0) − x − y r − x − y (cid:1) if x + y < r , . (40)where r = .
75. As in the previous examples, the di ff usion parameter is chosen such that d =
1, and t ∗ = ∆ t = .
1. The error is computed at t =
10, as in the previous simulation, as it is far enough from the inflection time t ∗ = ff erentiable hence smooth. However, ascan be seen from Fig. 7, along the circle x + y = r the solution changes drastically from zero to some finite non-zerovalue within a relatively small distance in the radial direction.Thus it is expected that most of the error of the finite element approximation concentrates along the circle. Thisbecomes evident in the distributions of the error as well as polynomial order over the mesh as displayed in the p-13 -5 -4 -3 -2 -1 -8 -7 -6 -5 -4 -3 -2 (a) (b) Figure 3: Comparison of convergence rates of the mixed (MFE) and standard (SFE) formulations with respect to the L error norm. Figure (a) isfor orders of approximation p = ,
2, while (b) is for p = ,
4. The legend ‘order’ stands for the absolute value of the slope of the convergencecurve once a consistent slope is established between consecutive refinements. -3 -2 -1 -6 -4 -2 (a) (b) Figure 4: Comparison of convergence rate the mixed (MFE) and standard (SFE) formulations with respect to the H -norm. Figure (a) is for orderof approximation p = ,
2, for m while (b) is for p = ,
4. The legend ‘order’ stands for the absolute value of the slope of the convergence curveonce a consistent slope is established between consecutive refinements. Degrees-of-freedom -6 -4 -2 E rr o r True error (uniform ref)A posteriori error (uniform ref)True error (adaptive ref)A posteriori error (adaptive ref)
Figure 5: Comparison of the p-adaptive and uniform refinement strategies for smooth and slowly varying solution.Figure 6: Distribution of error and polynomial order over the mesh at each p-adaptive iteration igure 7: Bump function (left) and norm of the gradient (right) Degrees-of-freedom -1 E rr o r True error (uniform ref)A posteriori error (uniform ref)True error (adaptive ref)A posteriori error (adaptive ref)
Figure 8: Convergence comparison between uniform and adaptive p-refinements for the bump function solution. adaptive sequence as shown in Fig. 9. The corresponding convergence result, shown in Fig. 8, exhibits the betterperformance compared to the uniform p-adaptive strategy.
Rough solution with singularities
Consider a one-species reaction di ff usion problem on the square domain Ω = [ − , , with reaction kinetics ofFisher’s type, f = m [1 − m ] . (41)The domain is comprised of square patches with contrasting di ff usivities D = d ( x ) I with a checkerboard pattern, asshown in Fig. 10 (a). For the blue patches d = . d = . m = . Γ (i.e., ¯ h = t = ∆ t n = .
1. Both the standard and mixed solutions werecomputed and compared. Because of the heterogeneous di ff usivity, the solution develops kinks along the interfacesof the patches and singularities at the corners. It is important to note the well-known fact that such irregularities(singularities) cannot be resolved by increasing the polynomial order. A feasible way of resolving such features isusing local h-adaptivity. In fact, numerical experimentation (not presented here) showed that p-adaptivity caused16 igure 9: Distribution of the error (top row) and polynomial order (bottom row) in a p-adaptive sequence. artificial oscillation near the corners as the polynomial order increases locally. For example, Fig. 11 (b) shows thedistribution of the flux magnitude, computed using the mixed method with the a priori adaptively refined mesh asshown in Fig. 11 (a). Figs. 11 (c) and (d) show the mass distribution computed using the standard and mixed methods,respectively. Even though the reference mesh Fig. 11 (a) has been used in both cases, the di ff erence in their respectivesolutions is apparent. This is due to the fact that the standard formulation uses a H -conforming space and is unableto approximate solutions with features such as discontinuities and singularities. By contrast, the mixed formulationuses a non-conforming L ( Ω ) space for m h and the exact, physically motivated conformity for the flux h h , i.e., H (div ).This allows the mixed method to capture such features.Fig. 12 shows the superiority of the mixed method over the standard single-field formulation. Here, the massprofiles along a line segment, coloured in red in Fig. 11 (a), that connects the centre and the right top corner of thesquare domain Ω are displayed. Along this line segment, there are two interior corners of patches where discontinuitiesin m h are expected. The solutions were obtained using various meshes at di ff erent levels of refinement with parameter h = / , / , /
20, and 1 /
40, including a reference mesh (denoted by Ref. Mesh), which is obtained by refiningalong the interfaces and corners of the patches. The discontinuities at the interior corners are captured almost exactlyusing the mixed method regardless of the refinement level. By contrast, using the standard formulation, none ofthe meshes resulted in a reasonable approximation of the discontinuities. Another interesting observation is that theapproximation of the mixed formulation converges from below. This is opposite to the standard formulation in whichthe approximation overestimates the solution.
The one-species Fisher-type equation (41) supports travelling wave solutions. The wave nature of the solutiondepends on the relative size of the reaction and di ff usion terms. When the reaction term is dominant, the wave frontsteepens and the wave travels with a finite speed. By contrast, when di ff usion is dominant, the influence of the reactionterm becomes less and the solution exhibits typical di ff usion behaviour, that is, it decays exponentially.Consider the planar domain as shown in Fig. 13, composed of rectangular patches that are arranged horizontallywith increasing di ff usivity between successive patches. The rectangular domain has height 0 . .
6, 0 .
4, 0 .
4, 0 .
4, and 2 .
2. The maximum di ff usivity d = . m = m is set to zero at t =
0. Since the di ff usivity is smallest in the firstleft patch, the solution starts to evolve slowly with a sharp wave front. As the wave passes each interface its speed17a) (b) Figure 10: The heterogeneous domain
Ω = [ − , with a checkerboard pattern where in (a) the blue regions have d = × − and the rest have d = × − , and in (b) an initial value m = . increases while the sharpness of the wavefront decreases. The wavefront is identified using a levelset method basedon a mass concentration value of 0 .
6. Uniform time steps of length ∆ t = .
05 have been used. The position of thewavefront against time is presented in Fig. 14 for various levels of mesh refinement. Importantly, the approximationsbased on the mixed method are accurate and converge to the correct wave speed. However, the approximation usingthe standard method initially overestimates the true speed of the travelling wave solution, and only slowly convergesto the correct speed. The p-adaptive mixed formulation is used with the second coarsest mesh, and it is also obtainedthat the speed of the wave is in good agreement with that of the converged results of the mixed or standard methods.
The class of problems considered here have application in various important areas including biological patternformation, morphogenesis [2] and electrophysiology [50].
Segregation pattern
A competition-di ff usion model involving three interacting species is considered. The level and mode of interactionbetween the species is the same. This, in e ff ect, means that the magnitude of each species that is consumed by theothers is the same as the other species that consumes it. The reaction term, for each i = , ,
3, is given by f i = m i [1 − a i m − a i m − a i m ] , (42)where the parameters in the model are represented in Table 1. It is assumed that all the three species have the samemobility rates, i.e., D = d I , where d = .
01. In cases where the dynamics is largely influenced by the reactionterm, it is important to analyse the local stability of the spatially homogeneous problem (i.e., ignoring the di ff usionterms). Such an analysis provides important insights into the range of parameter values for various possible spatio-temporal interaction patterns. A local stability analysis of the problem described by equations (42) reveals eightequilibrium points of which only, namely [ m , m , m ] = [1 / a , , m , m , m ] = [0 , / a , m , m , m ] = [0 , , / a ] are quasi-stable. In the case of the segregation problem, these equilibrium pointsrepresent regions which are exclusively occupied by one of the species.It is assumed initially that all the three species are distributed randomly over the domain Ω = [ − , as shownin Fig. 15. Such problems have been studied previously [25, 51, 52]. In those studies the numerical approaches wereeither the finite di ff erence method or the standard finite element method.18 igure 11: Comparison of approximations by standard and mixed methods, (a) an a priori adaptively refined mesh on which the mass concentrationapproximation at t = t =
6, where singularities at the corners of the patches are shown.Table 1: List of parameters for three species segregation problem d d d a a a a a a a a a .
01 0 .
01 0 .
01 1 3 3 3 1 3 3 3 119 igure 12: Mass concentration profile along a diagonal line from the centre of the square to the right corner (as indicated in Fig. 11(a) with a redline segment) at t =
6. Various uniform meshes at di ff erent refinement levels described by the mesh parameter h along with a reference mesh(designated as Ref. Mesh), which was adaptively refined along the interfaces and corners of the patches, were used.Figure 13: Domain for the travelling wave problem. Color indicates the di ff usivity, which is piecewise constant and increasing towards the right. igure 14: Comparison of speed of travelling wave solution with time. The broken lines correspond to the computation using the standardformulation (SFM) while the solid lines correspond to the mixed formulation (MFM). igure 15: Development of segregation pattern as a result of interactions of the three species (Equation (42), Table 1) at various times. The bluecolor surface plot represents the region dominated by m (i.e., m > m and m > m ), red color surface plot by m , and yellow surface plot by m . able 2: List of parameters for three species cyclic interaction problem d d d a a a a a a a a a .
01 0 .
01 0 .
01 1 2 7 7 1 2 2 7 1A relatively coarse mesh for such a problem is used, but to su ffi ciently represent the random initial condition andto captured the fast dynamics at the beginning stage order 6 polynomial approximation are used as shown in the lowerleft corner of Fig. 15. During the early stages of evolution, as show in Fig. 15 at t = t =
10, the dynamicsappears to be reasonably fast. Eventually, as shown in the second row of Fig. 15, as the species start to establishthemselves into well defined regions each occupied by one of the species, the interaction starts to proceed in a slowermanner. During this time the error distribution becomes more concentrated in the vicinity of the boundaries of theseregions. It is shown that the polynomial adaptation follows the error distribution very closely. These regions tend toa quasi-stable configuration, that is patches of convex shapes with triple junctions with angle of separation given by2 π/ Cyclic interaction
Here a three-species competition-di ff usion system of equations with the reaction term given by equation (42) isconsidered. The species react with each other in a cyclic way (based on parameters in Table 2) resulting in variouscomplex spatio-temporal patterns such as spiral-like, and band-like structures depending on the topology of the habitatand the initial configuration.A square habitat Ω = [ − , is considered with initial configuration as shown in Fig. 16 at t = ∆ t = .
2. Initially the polynomial order is set to 2 and adaptively increases to 6 (with parameter θ max = . θ min = .
1) throughout the simulation. As shown in Fig. 16, a spiral pattern starts to form turning in a clockwisedirection. The spiral shape consists of stripes of each species lying side-by-side, and eventually fills the region andcontinues with the same spiral feature. It is also shown that the error is high in the vicinity of the interface betweenthe species which led to the adaptivity taking place only on elements around such interfaces. This clearly shows thee ffi ciency of the p-adaptive mixed method. The propagation of ionic current in the cardiac muscle can be simulated using a monodomain model, which ismathematically equivalent to the reaction-di ff usion equation. The transmembrane electric potential can be viewedas a di ff using species, which “reacts” locally with the cellular ion channel densities. The reaction term depends onthe ion channel densities through a set of highly nonlinear and coupled ordinary di ff erential equations. Thus, theion channel densities can be viewed as non-di ff using species and treated simply as internal state variables. Thereare a large number of models available for the reaction term (called cardiac electrophysiology models) with varyingdegrees of complexity in terms of the number of (internal) variables. Using the proposed mixed method, we simulatethe phenomenon of spiral wave re-entry — the cause of several cardiac arrhythmias, such as ventricular tachycardia,atrial flutter, and atrial and ventricular fibrillation.A square block of cardiac tissue of dimension 100 mm is considered. The domain is subdivided into a relativelycoarse triangular mesh. The propagation of the transmembrane electric potential (more commonly known as the actionpotential) is governed by dd τ m − div( D ∇ m ) = f ( m , r ) + I Ω (cid:48) ( τ ) , (43)where the non-dimensional variable m ∈ [0 ,
1] is related to the transmembrane action potential E [mV] through therelation E = [100 m −
80] mV , igure 16: Development of spiral pattern as a result of cyclic interactions of three species (Equation (42), Table 2) at various stages. Region withblue depicts m , red m , and green m . able 3: List of parameters for spiral wave re-entry problem d α γ b c µ µ [mm ] [ − ] [ − ] [ − ] [ − ] [ − ] [ − ]0 .
01 0 .
01 0 .
01 1 2 7 7 r is a single internal variable representing the density of ionic channels, and I Ω (cid:48) is the external stimulus. The time t [ms] is non-dimensionalised as t = . τ ms . One of the simplest models capable of reproducing the spiral wave re-entry, the Aliev-Panfilov model [12], is adoptedfor the reaction term: f ( m , r ) = cm [ m − α ][1 − m ] − rm . (44)Equation (44) is supplemented by an ordinary di ff erential equation for the recovery (internal state) variables r :dd τ r = (cid:20) γ + µ r µ + m (cid:21) [ − r − cm [ m − b − , (45)where the parameters appearing in equations (44) and (45) are given in Table 3. We assume the conductivity to beisotropic, i.e., D = d I . The simulation is carried out using the IMEX mixed formulation with order k = k is the polynomial order used for the approximation of m and k + θ max = . θ min = .
03 starts with uniformly order 2 and increases locally to order 6.A horizontal planar wave is initiated by setting the action potential to E = −
40 mV on the region between y = y =
3. The wave form continues to propagate upwards as seen from the snapshot at t =
160 ms.Before the depolarising tail disappears, an external stimulus I Ω (cid:48) is applied to the strip of region, defined by Ω (cid:48) = { ( x , y ) : 50 < x < , < y < } , in order to initiate the spiral wave re-entry. The stimulus has magnitude 40 andis applied at t =
565 ms for a duration of 10 ms. This results in the development of the wavebreak (shown in Fig. 17at t = t =
924 ms,and continues afterwards.The computational aspects of this problem have been considered by several researchers in the electrophysiologyand electromechanics community. One notable work is by G¨oktepe and Kuhl [50] in which they used the standardfinite element approach with an implicit time integration scheme on a structured quad mesh. The proposed p-adaptive,IEMX mixed, as compared to their work, method is seen to capture the spiral wave re-entry dynamics more e ffi ciently.
8. Conclusion
A family of p-adaptive implicit-explicit, mixed finite element formulations has been proposed for a general classof reaction-di ff usion based problems. In contrast to single-field, standard finite element formulations, this class ofmethods provides accurate approximations of a wider class of solutions, including those with less regularity. Astandard formulation was shown to converge poorly, if at all, for such problems. The IMEX approach has beenshown to be e ffi cient and eliminates the dependence of algorithmic stability on the size of the spatial mesh size byhandling the non-local di ff usion part implicitly. This advantageous feature allows for mesh refinement, for examplein an adaptive strategy, without the need for changing the time step size, ∆ t . The explicit treatment of the localreaction term makes the implementation generic and modular for various classes of reaction kinetics as demonstratedby the wide range of problems that have been analysed in this paper. The finite element spaces are built using ahierarchical construction which, in addition to o ff ering optimal conditioning of the resulting linear system, makes theuse of the p -adaptivity strategy a natural choice [39]. The mixed formulation introduces additional DoFs. However,the computational complexity due to this increase in DoFs can be handled e ffi ciently using static condensation as themass concentration field (which is in L ) can be inverted locally since the local contributions are decoupled from oneanother. Moreover, this local inversion can also be used in block iterative schemes that involves computation of the25 igure 17: Evolution of a planar wave into a rotating spiral wave re-entry as a result of external stimulation I Ω (cid:48) =
40 applied on the red shadedregion (top left) during the time interval t =
570 ms to t =
580 ms. The planar wave is initiated with an initial excitation of E ( t = = −
40 mV onthe region shaded in blue.
Schur complement as an intermediate step. The Schur complement can then be computed exactly, resulting in a sparseglobal structure, rather than reverting to the common practice of approximating it.A distinguishing feature of the mixed formulation is that it leads to straightforward derivation and implemetation ofresidual based a posteriori error estimations without the need for computationally demanding postprocessing e ff ort asit is usually the case in literature, see, for example, [53–55]. This feature together with the hierarchical approximationof the mixed finite element method is exploited in formulating the p-adaptive strategy. It has been demonstrated by arange of examples that the p-adaptive algorithm performs very well in e ffi ciently resolving fine features.The performance of the proposed formulation is demonstrated by a number of challenging examples. The advan-tages of the proposed method over the standard techniques are showcased by the following two examples: a problemthat has singularities (see Section 7.1); and one that supports travelling wave solutions (see Section 7.2). The capabil-ity of this general p-adaptive framework is demonstrated by applying it to problems arising from di ff erent applicationssuch as electrophysiology [12, 50], and spatial pattern formation in theoretical ecology [6, 51, 52].Through a generalisation straightforwardly the proposed mixed method can be coupled to the mechanical defor-mation field for applications in cardiac electromechanics and chemo-mechanics. Due to the explicit treatment of thereaction term, the approach can be easily linked with, for example, electrophysiology models in the CellML repository[55]. Similarly, due to this explicit treatment of the reaction term, our computational approach can be easily used todrive the form of reaction kinetics models from experimental data following the approach proposed in [56]. Acknowledgements
The authors gratefully acknowledge the support provided by the EPSRC Strategic Support Package: Engineeringof Active Materials by Multiscale / Multiphysics Computational Mechanics - EP / R008531 / eferencesReferences [1] M. Chaplain, M. Ganesh, I. Graham, Spatio-temporal pattern formation on spherical surfaces: Numerical simulation and application to solidtumour growth, Journal of Mathematical Biology 42 (5) (2001) 387–423.[2] K. Garikipati, Perspectives on the mathematics of biological patterning and morphogenesis, Journal of the Mechanics and Physics of Solids99 (2017) 192 – 210.[3] P. K. Tapaswi, A. K. Saha, Pattern formation and morphogenesis: A reaction-di ff usion model, Bulletin of Mathematical Biology 48 (2) (1986)213 – 228.[4] S. F. Gilbert, Mathematical Modeling of Development, developmental biology. 6th edition Edition, Sinauer Associates, 2000.[5] D. Ambrosi, G. A. Ateshian, E. M. Arruda, S. C. Cowin, J. Dumais, A. Goriely, G. A. Holzapfel, J. D. Humphrey, R. Kemkemer, E. Kuhl,J. E. Olberding, L. A. Taber, K. Garikipati, Perspectives on biological growth and remodeling, Journal of the Mechanics and Physics of Solids59 (4) (2011) 863 – 883.[6] Y. Morishita, Y. Iwasa, Growth based morphogenesis of vertebrate limb bud, Bulletin of Mathematical Biology 70 (7) (2008) 1957–1978.[7] M. Tewary, J. Ostblom, L. Prochazka, T. Zulueta-Coarasa, N. Shakiba, R. Fernandez-Gonzalez, P. W. Zandstra, A stepwise model of reaction-di ff usion and positional information governs self-organized human peri-gastrulation-like patterning, Development 144 (23) (2017) 4298–4312.[8] A. M. Sebastian, M. Adrian, A two-scale reaction–di ff usion system with micro-cell reaction concentrated on a free boundary, ComptesRendus M´ecanique 336 (6) (2008) 481–486.[9] M. D. Ryser, S. V. Komarova, N. Nigam, The cellular dynamics of bone remodeling: A mathematical model, SIAM Journal on AppliedMathematics 70 (6) (2010) 1899–1921.[10] J. S. Matthew, A. L. Kerry, F. N. Donald, Chemotactic and di ff usive migration on a nonuniformly growing domain: numerical algorithmdevelopment and applications, Journal of Computational and Applied Mathematics 192 (2) (2006) 282 – 300.[11] R. C. P. Kerckho ff s, S. N. Healy, T. P. Usyk, A. D. McCulloch, Computational methods for cardiac electromechanics, Proceedings of theIEEE 94 (4) (2006) 769–783.[12] R. Aliev, A. Panfilov, A simple two-variable model of cardiac excitation, Chaos Solitons and Fractals 7 (3) (1996) 293–301.[13] Y. Zhang, X.-Q. Zhao, A reaction-di ff usion lyme disease model with seasonality, SIAM Journal on Applied Mathematics 73 (6) (2013)2077–2099.[14] W. Wang, X.-Q. Zhao, A nonlocal and time-delayed reaction-di ff usion model of dengue transmission, SIAM Journal on Applied Mathematics71 (1) (2011) 147–168.[15] R. E. Wilson, V. Capasso, Analysis of a reaction-di ff usion system modeling man–environment–man epidemics, SIAM Journal on AppliedMathematics 57 (2) (1997) 327–346.[16] M. Hemami, K. Parand, J. A. Rad, Numerical simulation of reaction-di ff usion neural dynamics models and their synchroniza-tion / desynchronization: Application to epileptic seizures, Computers and Mathematics with Applications 78 (11) (2019) 3644–3677.[17] D. Olmos, B. D. Shizgal, Pseudospectral method of solution of the fitzhugh–nagumo equation, Mathematics and Computers in Simulation79 (7) (2009) 2258–2278.[18] J. Ramos, A review of some numerical methods for reaction-di ff usion equations, Mathematics and Computers in Simulation 25 (6) (1983)538–548.[19] D. Barkley, A model for fast computer simulation of waves in excitable media, Physica D: Nonlinear Phenomena 49 (1) (1991) 61–70.[20] R. Ruiz-Baier, Primal-mixed formulations for reaction-di ff usion systems on deforming domains, Journal of Computational Physics 299(2015) 320–338.[21] J. Lang, Adaptive fem for reaction-di ff usion equations, Applied Numerical Mathematics 26 (1) (1998) 105–116.[22] N. Tuncer, A. Madzvamuse, A. Meir, Projected finite elements for reaction–di ff usion systems on stationary closed surfaces, Applied Numer-ical Mathematics 96 (2015) 45–71.[23] G. MacDonald, J. A. Mackenzie, M. Nolan, R. H. Insall, A computational method for the coupled solution of reaction–di ff usion equationson evolving domains and manifolds: Application to a model of cell migration and chemotaxis, Journal of Computational Physics 309 (2016)207–226.[24] C. Landsberg, A. Voigt, A multigrid finite element method for reaction-di ff usion systems on surfaces, Computing and Visualization in Science13 (2010) 177–185.[25] W. D. Mergia, K. C. Patidar, High-order semi-implicit linear multistep lg scheme for a three species competition-di ff usion system in two-dimensional spatial domain arising in ecology, Communications in Nonlinear Science and Numerical Simulation 84 (2020) 1–16.[26] D. Bo ffi , F. Brezzi, M. Fortin, Mixed finite element methods and applications, Springer Series in Computational Mathematics, Springer,Berlin, 2013.[27] L. P. Franca, T. J. Hughes, Two classes of mixed finite element methods, Computer Methods in Applied Mechanics and Engineering 69 (1)(1988) 89–129.[28] D. N. Arnold, Mixed finite element methods for elliptic problems, Computer Methods in Applied Mechanics and Engineering 82 (1) (1990)281–300.[29] H. Fu, H. Guo, J. Hou, J. Zhao, A stabilized mixed finite element method for steady and unsteady reaction–di ff usion equations, ComputerMethods in Applied Mechanics and Engineering 304 (2016) 102 – 117.[30] S. Liu, Y. Chen, A new two-grid method for expanded mixed finite element solution of nonlinear reaction di ff usion equations, Advances inApplied Mathematics and Mechanics 9 (3) (2017) 757–774.[31] S. J. Ruuth, Implicit-explicit methods for reaction-di ff usion problems in pattern formation, Journal of Mathematical Biology 34 (2) (1995)148–176.
32] U. M. Ascher, S. J. Ruuth, B. T. R. Wetton, Implicit-explicit methods for time-dependent partial di ff erential equations, SIAM Journal onNumerical Analysis 32 (3) (1995) 797–823.[33] H. Zhang, A. Sandu, Application of implicit-explicit general linear methods to reaction-di ff usion problems, AIP Conference Proceedings1648 (1) (2015) 150017.[34] I. Farag´o, F. Izs´ak, T. Szab´o, A. Kriston, An imex scheme for reaction-di ff usion equations: application for a pem fuel cell model, OpenMathematics 11 (4) (2013) 746–759.[35] O. Lakkis, A. Madzvamuse, C. Venkataraman, Implicit-explicit timestepping with finite element approximation of reaction-di ff usion systemson evolving domains, SIAM Journal on Numerical Analysis 51 (4) (2013) 2309–2330.[36] J. Lin, S. Reutskiy, A cubic b-spline semi-analytical algorithm for simulation of 3d steady-state convection-di ff usion-reaction problems,Applied Mathematics and Computation 371 (2020) 124944.[37] C. A. Kennedy, M. H. Carpenter, Additive runge–kutta schemes for convection–di ff usion–reaction equations, Applied Numerical Mathematics44 (1–2) (2003) 139–181.[38] C. A. Kennedy, M. H. Carpenter, Higher-order additive runge–kutta schemes for ordinary di ff erential equations, Applied Numerical Mathe-matics 136 (2019) 183–205.[39] M. Ainsworth, J. Coyle, Hierarchic finite element bases on unstructured tetrahedral meshes, International Journal for Numerical Methods inEngineering 58 (14) (2003) 2103–2130.[40] B. Guo, I. Babuska, S. N. Atluri, The h-p version of the finite element method. i: The basic approximation results, Computational Mechanics1 (1986) 203–220.[41] I. Babuska, W. Gui, The h, p and h-p versions of the finite element method in 1 dimension. part iii. the adaptive h-p version, NumerischeMathematik 49 (1986) 659–684.[42] F. Fuentes, B. Keith, L. Demkowicz, S. Nagaraj, Orientation embedded high order shape functions for the exact sequence elements of allshapes, Computers & Mathematics with Applications 70 (4) (2015) 353–458.[43] W. D¨orfler, V. Heuveline, Convergence of an adaptive hp finite element strategy in one space dimension, Applied Numerical Mathematics57 (10) (2007) 1108–1124.[44] Ł. Kaczmarczyk, Z. Ullah, K. Lewandowski, X. Meng, X.-Y. Zhou, I. Athanasiadis, H. Nguyen, C.-A. Chalons-Mouriesse, E. Richardson,E. Miur, A. Shvarts, M. Wakeni, C. Pearce, MoFEM: an open source, parallel finite element library, The Journal of Open Source Software.[45] T. J. Tautges, C. Ernst, C. Stimpson, R. J. Meyers, K. Merkley, Moab : a mesh-oriented database (2004).URL [46] T. Tautges, Moab-sd: integrated structured and unstructured mesh representation, Engineering with Computers 20 (2004) 286–293.[47] S. Balay, W. D. Gropp, L. C. McInnes, B. F. Smith, E ffi cient management of parallelism in object oriented numerical software libraries, in:E. Arge, A. M. Bruaset, H. P. Langtangen (Eds.), Modern Software Tools in Scientific Computing, Birkh¨auser Press, 1997, pp. 163–202.[48] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev,D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang,H. Zhang, PETSc Web page (2019).URL [49] S. Abhyankar, J. Brown, E. M. Constantinescu, D. Ghosh, B. F. Smith, H. Zhang, Petsc / ts: A modern scalable ode / dae solver library, arXivpreprint arXiv:1806.01437.[50] S. G¨oktepe, E. Kuhl, Computational modeling of cardiac electrophysiology: A novel finite element approach, International Journal forNumerical Methods in Engineering 79 (2) (2009) 156–178.[51] M. Mimura, M. Tohma, Dynamic coexistence in a three-species competition–di ff usion system, Ecological Complexity 21 (2015) 215–232.[52] M. Mimura, Y. Kan-on, Dynamic coexistence in a three-species competition–di ff usion system, Predation-Mediated Coexistence and Segre-gation Structures 18 (1986) 129–155.[53] M. Ainsworth, A posteriori error estimation for lowest order raviart–thomas mixed finite elements, SIAM Journal on Scientific Computing30 (1) (2008) 189–204.[54] D. Braess, R. Verfurth, A posteriori error estimators for the raviart-thomas element, SIAM Journal on Numerical Analysis 33 (6) (1996)2431–2444.[55] C. M. Lloyd, T. Yu, Cellml model repository, in: W. Dubitzky, O. Wolkenhauer, K.-H. Cho, H. Yokota (Eds.), Encyclopedia of SystemsBiology, Springer New York, New York, NY, 2013, pp. 376–378.[56] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems,Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937.usion system, Predation-Mediated Coexistence and Segre-gation Structures 18 (1986) 129–155.[53] M. Ainsworth, A posteriori error estimation for lowest order raviart–thomas mixed finite elements, SIAM Journal on Scientific Computing30 (1) (2008) 189–204.[54] D. Braess, R. Verfurth, A posteriori error estimators for the raviart-thomas element, SIAM Journal on Numerical Analysis 33 (6) (1996)2431–2444.[55] C. M. Lloyd, T. Yu, Cellml model repository, in: W. Dubitzky, O. Wolkenhauer, K.-H. Cho, H. Yokota (Eds.), Encyclopedia of SystemsBiology, Springer New York, New York, NY, 2013, pp. 376–378.[56] S. L. Brunton, J. L. Proctor, J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems,Proceedings of the National Academy of Sciences 113 (15) (2016) 3932–3937.