A bounded numerical solver for a fractional FitzHugh–Nagumo equation and its high-performance implementation

Abstract

In this paper, we propose a high-performance implementation of a space-fractional FitzHugh–Nagumo model. Our implementation is based on a positivity- and boundedness-preserving finite-difference model to approximate the solutions of a Riesz space-fractional reaction-diffusion equation. The model generalizes the FitzHugh–Nagumo model. The stability and convergence of the difference scheme are thoroughly discussed. Moreover, we prove the existence and uniqueness of numerical solutions, positivity, boundedness and consistency of the model. The scheme is based on weighted and shifted Grünwald differences. The conjugate gradient method is used then to solve the sparse matrix system. The MPI and PETSc libraries are used for the computational implementation. We investigate the influence of some computer factors on the performance of our implementation and scalability. More precisely, we consider the number of cores, the size of the computation mesh and the orders of the fractional derivatives. Tests are evaluated on a ccNUMA architecture with two CPUs.

Implicit formulas

The purpose of this section is to provide explicit forms of the finite-difference method (29) under all the possible cases. More precisely, note that the right-hand side of that equation can be written as

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {{\mathfrak {a}} _{0 , 0} ^k v _{0 , 0} ^{k + 1} - {\mathfrak {b}} _1 v _{1 , 0} ^{k + 1} - {\mathfrak {b}} _2 v _{0 , 1} ^{k + 1} - \sum _{l = 2} ^M {\widetilde{\omega }} _{l + 1} ^1 v _{l , 0} ^{k + 1} - \sum _{l = 2} ^N {\widetilde{\omega }} _{l + 1} ^2 v _{0 , l} ^{k + 1}}, \quad \forall (m , n ) = (0 , 0), \\ \displaystyle {{\mathfrak {a}} _{0 , n} ^k v _{0 , n} ^{k + 1} - {\mathfrak {b}} _1 v _{1 , n} ^{k + 1} - {\mathfrak {b}} _2 \left( v _{0 , n - 1} ^{k + 1} + v _{0 , n + 1} ^{k + 1} \right) - \sum _{l = 2} ^M {\widetilde{\omega }} _{l + 1} ^1 v _{l , n} ^{k + 1} - \sum _{\begin{array}{c} l = 0 \\ l \ne n - 1 , \\ n , n + 1 \end{array}} ^N {\widetilde{\omega }} _{|n - l| + 1} ^2 v _{0 , l} ^{k + 1}}, \quad \forall (m , n) \in 0 \times I _{N - 1}, \\ \displaystyle {{\mathfrak {a}} _{0 , N} ^k v _{0 , N} ^{k + 1} - {\mathfrak {b}} _1 v _{1 , N} ^{k + 1} - {\mathfrak {b}} _2 v _{0 , N - 1} ^{k + 1} - \sum _{l = 2} ^M {\widetilde{\omega }} _{l + 1} ^1 v _{l , N} ^{k + 1} - \sum _{l = 0} ^{N - 2} {\widetilde{\omega }} _{N - l + 1} ^2 v _{0 , l} ^{k + 1}}, \quad \forall (m , n) = (0 , N), \\ \displaystyle {{\mathfrak {a}} _{m , 0} ^k v _{m , 0} ^{k + 1} - {\mathfrak {b}} _1 \left( v _{m - 1 , 0} ^{k + 1} + v _{m + 1 , 0} ^{k + 1} \right) - {\mathfrak {b}} _2 v _{m , 1} ^{k + 1} - \sum _{\begin{array}{c} l = 0 \\ l \ne m - 1 , \\ m , m + 1 \end{array}} ^M {\widetilde{\omega }} _{|m - l| + 1} ^1 v _{l , 0} ^{k + 1} - \sum _{l = 2} ^N {\widetilde{\omega }} _{l + 1} ^2 v _{m , l} ^{k + 1}}, \quad \forall (m , n) \in I _{M - 1} \times 0, \\ \displaystyle {{\mathfrak {a}} _{m , n} ^k v _{m , n} ^{k + 1} - {\mathfrak {b}} _1 \left( v _{m - 1 , n} ^{k + 1} + v _{m + 1 , n} ^{k + 1} \right) - {\mathfrak {b}} _2 \left( v _{m , n - 1} ^{k + 1} + v _{m , n + 1} ^{k + 1} \right) } \\ \displaystyle {\qquad - \sum _{\begin{array}{c} l = 0 \\ l \ne m - 1 , \\ m , m + 1 \end{array}} ^M {\widetilde{\omega }} _{|m - l| + 1} ^1 v _{l , n} ^{k + 1} - \sum _{\begin{array}{c} l = 0 \\ l \ne n - 1 , \\ n , n + 1 \end{array}} ^N {\widetilde{\omega }} _{|n - l| + 1} ^2 v _{m , l} ^{k + 1}}, \quad \forall (m , n) \in J, \\ \displaystyle {{\mathfrak {a}} _{m , N} ^k v _{m , N} ^{k + 1} - {\mathfrak {b}} _1 \left( v _{m - 1 , N} ^{k + 1} + v _{m + 1 , N} ^{k + 1} \right) - {\mathfrak {b}} _2 v _{m , N - 1} ^{k + 1} - \sum _{\begin{array}{c} l = 0 \\ l \ne m - 1 , \\ m , m + 1 \end{array}} ^M {\widetilde{\omega }} _{|m - l| + 1} ^1 v _{l , N} ^{k + 1}} \\ \qquad \displaystyle {- \sum _{l = 0} ^{N - 2} {\widetilde{\omega }} _{N - l + 1} ^2 v _{m , l} ^{k + 1}}, \quad \forall (m , n) \in I _{M - 1} \times N, \\ \displaystyle {{\mathfrak {a}} _{M , 0} ^k v _{M , 0} ^{k + 1} - {\mathfrak {b}} _1 v _{M - 1 , 0} ^{k + 1} - {\mathfrak {b}} _2 v _{M , 1} ^{k + 1} - \sum _{l = 0} ^{M - 2} {\widetilde{\omega }} _{M - l + 1} ^1 v _{l , 0} ^{k + 1} - \sum _{l = 2} ^N {\widetilde{\omega }} _{l + 1} ^2 v _{M , l} ^{k + 1}}, \quad \forall (m , n) = (M , 0), \\ \displaystyle {{\mathfrak {a}} _{M , n} ^k v _{M , n} ^{k + 1} - {\mathfrak {b}} _1 v _{M - 1 , n} ^{k + 1} - {\mathfrak {b}} _2 \left( v _{M , n - 1} ^{k + 1} + v _{M , n + 1} ^{k + 1} \right) - \sum _{l = 0} ^{M - 2} {\widetilde{\omega }} _{M - l + 1} ^1 v _{l , n} ^{k + 1}} \\ \quad \displaystyle {- \sum _{\begin{array}{c} l = 0 \\ l \ne n - 1 , \\ n , n + 1 \end{array}} ^N {\widetilde{\omega }} _{|n - l| + 1} ^2 v _{M , l} ^{k + 1}}, \quad \forall (m , n) \in M \times I _{N - 1}, \\ \displaystyle {{\mathfrak {a}} _{M , N} ^k v _{M , N} ^{k + 1} - {\mathfrak {b}} _1 v _{M - 1 , N} ^{k + 1} - {\mathfrak {b}} _2 v _{M , N - 1} ^{k + 1} - \sum _{l = 0} ^{M - 2} {\widetilde{\omega }} _{M - l + 1} ^1 v _{l , N} ^{k + 1} - \sum _{l = 0} ^{N - 2} \widetilde{\omega } _{N - l + 1} ^2 v _{M , l} ^{k + 1}}, \quad \forall (m , n) = (M , N). \end{array} \right. \end{aligned}$$

(A.1)