A Jacobian-Free Newton–Krylov Method to Solve Tumor Growth Problems with Effective Preconditioning Strategies

Abstract

A Jacobian-free Newton–Krylov (JFNK) method with effective preconditioning strategies is introduced to solve a diffusion-based tumor growth model, also known as the Fisher–Kolmogorov partial differential equation (PDE). The time discretization of the PDE is based on the backward Euler and the Crank–Nicolson methods. Second-order centered finite differencing is used for the spatial derivatives. We introduce two physics-based preconditioners associated with the first- and second-order temporal discretizations. The theoretical time and spatial accuracies of the numerical scheme are verified through convergence tables and graphs that correspond to different computational settings. We present efficiency studies with and without using the preconditioners. Our numerical findings indicate the excellent performance of the newly proposed preconditioning strategies. In other words, when we turn the preconditioners on, the average number of GMRES and the Newton iterations are significantly reduced.

1. Introduction

Gliomas are the most common and very aggressive type of brain tumor. They have been known to be incurable due to their large invasion capacity. It is important to have good mathematical models and their numerical solutions to understand the true nature and the underlying physical mechanism of this phenomenon. Over the years, researchers have paid a lot of attention to the modeling part and solution techniques. A set of several review articles on this topic can be found in the following references [1,2,3,4,5]. Widely used mathematical models are based on the nonlinear reaction–diffusion type partial differential equations (PDEs) to describe the time evolution of the tumor growth. It is difficult to solve these models analytically because of the existence of nonlinearity and discontinuity in the source and diffusion coefficient terms. On the other hand, numerical solutions and computer simulations have been good alternatives. In the last two decades, there have been several attempts such as variational iteration method (VIM) [6], stochastic models [7,8], finite element methods (FEM) [9,10,11,12], finite difference schemes (FDSs) [6,11,13,14], cubic B-spline collocation methods [12], hybrid models [7,10], and Bayesian approaches [15] to approximately solve the mathematical model and simulate the time behavior of the tumor concentrations in the glioma cells.

In 2015, Wu et al. [6] introduced a variational iteration method (VIM) to study this problem by transforming the reaction–diffusion equation into an equivalent integral equation. This method enables one to avoid the treatment of the time derivatives as in classical finite difference sense. However, the method involves the so-called Picard iteration procedure that can be computationally very inefficient.

Later on, Mach et al. [9] investigated the finite-element nonlinear Galerkin method (FEM) in one spatial dimension and its application to the selected reaction–diffusion systems. Another FEM approach was carried out by Dehghan and Narimani [16] in 2018 and called the element-free Galerkin method to simulate the effects of cancer cell invasion on the surrounding tissue. We note that even though the FEM is preferable for irregular geometries compared to the FDS, the FEM is more prone to stability issues especially when the discontinuous diffusion coefficients and sharp gradients exist in the system.

In 2020, Dehghan and Narimani [17] used a collocation meshless technique called the radial basis function-generated finite difference (RBF-FD) scheme to solve the reaction–diffusion model in multi-dimensional settings. They employed a first-order semi-implicit scheme to approximate the time variable and solved the resulting linear system of algebraic equations using the biconjugate gradient stabilized (BiCGSTAB) iterative scheme with a zero-fill incomplete lower–upper (ILU) preconditioner. An obvious drawback of this approach is that, despite being incomplete, it involves an expensive LU decomposition, let alone the positive-definite requirement of the CG method.

Another study was carried out by Kapoor and Joshi in [12]. They used a modified cubic uniform algebraic hyperbolic tension B-spline differential quadrature method to discretize the spatial partial derivatives in 1D and 2D reaction–diffusion systems. They utilized a strong stability preserving Runge–Kutta (SSP-RK4) time scheme to solve the resulting ordinary differential Equation (ODE). One disadvantage of this approach is that it can be expensive to implement due to the high number of function evaluations in the cubic spline procedure and due to the stability restrictions imposed by the explicit SSP-RK4 method.

In this study, we adopt the proliferation–invasion (PI) model that is assumed to capture the most fundamental features of cancer cells (e.g., cell growth cannot be restrained and can occupy surrounding tissues). In this model, the rate of change of cell density can be written as a sum of the invasion of cells into nearby tissue (heterogeneous diffusion coefficient times the cell density) and the proliferation of cells (in terms of logistic growth). The resulting equations are similar to the reaction–diffusion equation with radiotherapy and chemotherapy effects [18]. We note that we neglect the effects of radiotherapy and chemotherapy in this paper because, for the time being, we aim to introduce an accurate and effective solution methodology and test its convergence and efficiency on a simplified model. However, we will consider a more realistic model with all additional effects included in a follow-up study [17]. We also note that a one-dimensional simplified model was numerically studied by the co-author of this work in [19]. The major drawback of [19] is that it relies on the construction, storage, and inversion of the Jacobian matrix of the nonlinear system. This approach can be computationally very inefficient and very difficult to implement in a situation where discontinuous or solution-dependent diffusion coefficient exists together with nonlinear source terms, which is the case herein.

In this paper, we utilize a Jacobian-free Newton–Krylov (JFNK) method with new effective preconditioning strategies to solve the growth model. The JFNK methods are synergic combinations of Newton type-methods for the super linearly convergent solution of nonlinear equations and Krylov subspace methods for solving the linear equations [20]. The implementation of a JFNK method requires the solution of a linear system that comes from the Newton correction step. Generally, the linear system is solved based on a particular Krylov subspace method and the generalized minimal residual (GMRES) iterative method [20,21]. The main advantage of a Jacobian-free Newton–Krylov method is that it does not require the formation, storage, and inversion of the Jacobian matrix. Perhaps, this itself is the main reason why the JFNK methods became very popular in the scientific community. They have been successfully used to solve boundary or initial valued problems coming from diverse areas; compressible or incompressible fluid dynamics problems, heat transfer engineering and thermal hydraulic problems, neutronics applications, solid mechanics problems, astrophysical applications, and optimization problems from industry to name a few [20,22,23,24,25,26,27,28].

As part of the JFNK method, we need to include a preconditioning step to accelerate the convergence of the Krylov block. We note that the preconditioning step became almost standard for any type of iterative scheme [21,29,30]. Here, we introduce very effective preconditioning techniques that share similarities to those presented in [31,32,33,34,35] in such a way that they all rely on the underlying physics. In the literature, these kinds of methods are referred to as physics-based or physical preconditioning. The main benefit of this preconditioning strategy is that one can target a particular part of the model that creates stiffness or causes slow convergence and only apply the procedure to such a part rather than to the entire system. The numerical results in Section 4 indicate that the new preconditioners perform very well, decreasing the average number of GMRES and the Newton iterations to very low numbers.

The organization of this paper is as follows: In Section 2, we define the governing equations. In Section 3, we describe the Jacobian-free Newton–Krylov (JFNK) method together with the new physics-based preconditioning techniques. In Section 4, we present our numerical results to verify the numerical convergence of the proposed schemes and to show the performance and capabilities of the new preconditioning strategies. In Section 5, we summarize our concluding remarks and hint at possible future works.

2. Governing Equations

The governing PDEs are based on the proliferation–invasion (PI) model [17,18,19,36] to describe glioma cell invasion throughout the brain as a reaction–diffusion process;

$$\frac{\partial c(x,y,t)}{\partial t}=\frac{\partial }{\partial x}\left(D\frac{\partial c(x,y,t)}{\partial x}\right)+\frac{\partial }{\partial y}\left(D\frac{\partial c(x,y,t)}{\partial y}\right)+\rho c(x,y,t)\left(1-\frac{c(x,y,t)}{c_{\max }}\right),$$

(1)

where $c(x,y,t)$ is the tumor concentration of glioma cells at spatial locations $x,y$ and time t, $\rho$ is the net proliferation rate in units $da{y}^{-1}$, and $D(x,y)$ is the diffusion coefficient. Notice that we can solve Equation (1) as a one- or two- dimensional model by simply turning the y derivatives off and on. The birth and death rates are included in $\rho$, assuming a Verhulst (logistic growth) law including a tissue-carrying capacity $c_{\max }$ [37]. The diffusion coefficient can depend on the values of the concentration and may vary from patient to patient [38]. In this study, we use discontinuous diffusion coefficient profiles that contain low (gray region) and high (white region) diffusion zones [38]. In one- dimensional test cases, we use

$$\begin{array}{c}\hfill D\left(x\right)=\left\{\begin{array}{cc}0.13\hfill & \mathrm{if} 0\le x\le 7.5 (\mathrm{gray} \mathrm{region})\hfill \\ 0.65\hfill & \mathrm{if} 7.5<x\le 42.5 (\mathrm{white} \mathrm{region})\hfill \\ 0.13\hfill & \mathrm{if} 42.5<x\le 50 (\mathrm{gray} \mathrm{region}),\hfill \end{array}\right.\end{array}$$

(2)

where x belongs to the interval $I=[0,50]$, and in two -dimensional cases, as suggested by Shim et al. [39], we set

$$\begin{array}{c}\hfill D(x,y)=\left\{\begin{array}{cc}0.65\hfill & \mathrm{if} 0\le d\le 17.5 (\mathrm{white} \mathrm{region})\hfill \\ 0.13\hfill & \mathrm{otherwise} (\mathrm{gray} \mathrm{region}),\hfill \end{array}\right.\end{array}$$

(3)

where $d=\sqrt{{(x-x_0)}^2+{(x-y_0)}^2}$ and $\left(x_0{y}_0\right)=(25,25)$ is the center of the two-dimensional computational domain $D=[0,50]\times [0,50]$. The zero-flux boundary conditions and one and two- dimensional Gaussian initial tumor profiles are defined as

$$c_x(0,y,t)=c_x(50,y,t)=c_y(x,0,t)=c_y(x,50,t)=0,$$

(4)

$$c(0,t)=g\left(x\right)=\frac{1}{\sqrt{2\pi }ϵ}{e}^{-\frac{1}{2}{\left(\frac{x-x_0}{ϵ}\right)}^2}c(x,y,0)=g(x,y)=\frac{1}{\sqrt{2\pi }ϵ}{e}^{-\frac{d^2}{2{ϵ}^2}},$$

(5)

where we use $ϵ=0.01$ in our calculations. Figure 1 illustrates one- and two- dimensional Gaussian initial tumors and discontinuous diffusion coefficient profiles.

3. Numerical Algorithm

We present two implicit time discretization methods: one is based on the first-order backward Euler and the other is the second-order Crank–Nicolson scheme. Here, we purposely considered two different time-convergent schemes. The first reason is the fact that the local time error can accumulate differently depending on the time accuracy of the method. Furthermore, the second reason is that the time accuracy can affect the overall performance of the JFNK method. In Section 4, we present a head on comparison study for the two time methods and the aforementioned numerically verified behaviors. The spatial derivatives are discretized using the second-order centered finite differencing procedure on a staggered grid configuration. An obvious advantage of using a staggered grid is that it enables us to position the discontinuous diffusion coefficients exactly at the cell edges, rather than extrapolating them from the cell centers.

3.1. The Backward Euler (BE) Method

The first-order backward Euler (BE) time scheme together with the second-order space discretization of Equation (1) can be given in the following form

$$\begin{array}{c}\hfill \frac{c_{i,j}^{n+1}-c_{i,j}^n}{\Delta t}-\frac{D_{i+\frac{1}{2},j}{c}_{xi+\frac{1}{2},j}^{n+1}-D_{i-\frac{1}{2},j}{c}_{xi-\frac{1}{2},j}^{n+1}}{\Delta x}-\frac{D_{i,j+\frac{1}{2}}{c}_{yi,j+\frac{1}{2}}^{n+1}-D_{i,j-\frac{1}{2}}{c}_{yi,j-\frac{1}{2}}^{n+1}}{\Delta y}\\ \hfill -\rho c_{i,j}^{n+1}\left(1-\frac{c_{i,j}^{n+1}}{c_{\max }}\right)\approx 0,\end{array}$$

(6)

where $c_{i,j}^n$ approximates the true solution at the grid location $(x_i,y_j)$ and time $t_n$ (e.g., $c_{i,j}^n\approx c(x_i,y_j,t_n)=c(i\Delta x,j\Delta y,n\Delta t)$ and the following derivatives are approximated by

$$c_{xi+\frac{1}{2},j}^{n+1}\approx \frac{c_{i+1,j}^{n+1}-c_{i,j}^{n+1}}{\Delta x}{c}_{xi-\frac{1}{2},j}^{n+1}\approx \frac{c_{i,j}^{n+1}-c_{i-1,j}^{n+1}}{\Delta x}$$

$$c_{yi,j+\frac{1}{2}}^{n+1}\approx \frac{c_{i,j+1}^{n+1}-c_{i,j}^{n+1}}{\Delta y},c_{yi,j-\frac{1}{2}}^{n+1}\approx \frac{c_{i,j}^{n+1}-c_{i,j-1}^{n+1}}{\Delta y}.$$

Here, $\Delta x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ and $\Delta y=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$ are the size of the $(i,j)$th cell in x and y directions, $(x_i,y_j)$ represents the center of the $(i,j)$th cell, and $(x_{i\pm \frac{1}{2}},y_{j\pm \frac{1}{2}})$ denote the left, right, top, and bottom edges of the ${(i,j)}^{th}$ cell (Figure 2). The spatial domain $[0,50]\times [0,50]$ is divided into $M_x$ and $M_y$ sections in the x and y directions. We use the uniform grids in this study, implying that $(\Delta x_i,\Delta y_j)=(\Delta x,\Delta y)=(\frac{50}{M_x},\frac{50}{M_y})$ where $i=0,\dots M_x-1$ and $j=0,\dots M_y-1$.

We utilize the following discrete Neumann boundary conditions for Equation (4)

$$c_{x-\frac{1}{2},j}^{n+1}\approx 0,c_{x{M}_x-\frac{1}{2},j}^{n+1}\approx 0,c_{yi,-\frac{1}{2}}^{n+1}\approx 0,c_{yi,M_y-\frac{1}{2}}^{n+1}\approx 0.$$

(7)

Equation (7) can be used to set the ghost cell values in the x direction as

$$0\approx c_{x-\frac{1}{2},j}^{n+1}\approx (c_{0,j}^{n+1}-c_{-1,j}^{n+1})/\Delta x\Rightarrow c_{-1,j}^{n+1}\approx c_{0,j}^{n+1},$$

$$0\approx c_{x{M}_x-\frac{1}{2},j}^{n+1}\approx (c_{M_x,j}^{n+1}-c_{M_x-1,j}^{n+1})/\Delta x\Rightarrow c_{M_x,j}^{n+1}\approx c_{M_x-1,j}^{n+1},$$

and, similarly, in the y direction as

$$0\approx c_{yi,-\frac{1}{2}}^{n+1}\approx (c_{i,0}^{n+1}-c_{i,-1}^{n+1})/\Delta y\Rightarrow c_{i,-1}^{n+1}\approx c_{i,0}^{n+1},$$

$$0\approx c_{yi,M_y-\frac{1}{2}}^{n+1}\approx (c_{i,M_y}^{n+1}-c_{i,M_y-1}^{n+1})/\Delta y\Rightarrow c_{i,M_y}^{n+1}\approx c_{i,M_y-1}^{n+1}.$$

3.2. The Crank–Nicolson (CN) Method

The Crank–Nicolson method is a second-order accurate scheme in both time and space. It utilizes the time average of the flux and source terms together with the centered spatial discretizations. Along with these principles, the governing equation yields

$$\begin{array}{ccc}\hfill \frac{c_{i,j}^{n+1}-c_{i,j}^n}{\Delta t}& -& \frac{1}{2}\left(\frac{D_{i+\frac{1}{2},j}{c}_{xi+\frac{1}{2},j}^{n+1}-D_{i-\frac{1}{2},j}{c}_{xi-\frac{1}{2},j}^{n+1}}{\Delta x}+\frac{D_{i+\frac{1}{2},j}{c}_{xi+\frac{1}{2},j}^n-D_{i-\frac{1}{2},j}{c}_{xi-\frac{1}{2},j}^n}{\Delta x}\right)\hfill \\ & -& \frac{1}{2}\left(\frac{D_{i,j+\frac{1}{2}}{c}_{yi,j+\frac{1}{2}}^{n+1}-D_{i,j-\frac{1}{2}}{c}_{yi,j-\frac{1}{2}}^{n+1}}{\Delta y}+\frac{D_{i,j+\frac{1}{2}}{c}_{y,i,j+\frac{1}{2}}^n-D_{i,j-\frac{1}{2}}{c}_{yi,j-\frac{1}{2}}^n}{\Delta y}\right)\hfill \\ & -& \frac{1}{2}\left(\rho c_{i,j}^{n+1}\left(1-\frac{c_{i,j}^{n+1}}{c_{\max }}\right)+\rho c_{i,j}^n\left(1-\frac{c_{i,j}^n}{c_{\max }}\right)\right)\approx 0.\hfill \end{array}$$

(8)

Equations (6) and (8) can be viewed as a system of nonlinear algebraic equations (e.g., $\mathbf{F}\left(c\right)\approx 0$) that need to be solved for the unknown $c^{n+1}$. The resulting nonlinear equation system will be solved using the Jacobian-free Newton Krylov (JFNK) method [20,31,34]. After applying a desired number of nonlinear iterations, the JFNK method guarantees the convergence of all nonlinearities in the system. Of course, we have to keep in mind that, in order for the JFNK method—just like the other nonlinear solvers—to converge to a physically correct solution, we have to provide a good initial guess. This is not an issue here, because we are solving an initial value problem. In other words, at the beginning of each time update, we have a very good initial guess coming from the previous time step solution.

3.3. The Jacobian-Free Newton–Krylov Method

In the following section, we briefly describe the JFNK method. The JFNK method has two components, namely the Newton and Krylov blocks. The Newton block iteratively solves $\mathbf{F}\left(c\right)\approx 0$ over the sequence of a linear system defined by

$$\begin{array}{ccc}\hfill \mathbf{J}\left(c^k\right)\delta c^k& \approx & -\mathbf{F}\left(c^k\right),\hfill \\ \hfill c^{k+1}& \approx & c^k+\delta c^k,k=0,1,\cdots \hfill \end{array}$$

(9)

where $\mathbf{J}\left(c^k\right)=\frac{\partial \mathbf{F}}{\partial c}$ is the Jacobian matrix and $\delta c^k$ is the update vector. The Newton iteration is stopped if the following criterion is matched

$$\parallel \mathbf{F}\left(c^k\right){\parallel }_2<{\mathrm{tol}}_{res}{\parallel \mathbf{F}\left(c^0\right)\parallel }_2,$$

(10)

where ${\mathrm{tol}}_{res}$ is the user-defined tolerance, typically ${10}^{-8}$. The resulting linear system at the Newton correction step (Equation (9)) is solved utilizing a Krylov subspace method. In particular, we use the Arnoldi-based GMRES method [21]. At the beginning of the GMRES iteration, we need to define an initial residual, ${\mathbf{r}}_0$, for a given initial guess $\delta c_0$,

$${\mathbf{r}}_0\approx -\mathbf{F}\left(c\right)-\mathbf{J}\delta c_0.$$

(11)

Since the Krylov (GMRES) iteration is executed at a fixed k, we omit the index k for convenience in Equation (11). Let j denote the Krylov iteration index. The jth Krylov iteration minimizes the quantity $\parallel \mathbf{J}\delta c_j{+\mathbf{F}\left(c\right)\parallel }_2$ within a subspace of small dimensions relative to n (the number of unknown vectors) in a least-squares sense. $\delta c_j$ is drawn from the subspace spanned by the Krylov vectors, $\{{\mathbf{r}}_0,\mathbf{J}{\mathbf{r}}_0,{\mathbf{J}}^2{\mathbf{r}}_0,\cdots ,{\mathbf{J}}^{j-1}{\mathbf{r}}_0\}$, and can be written as the following combination

$$\delta c_j\approx \delta c_0+\sum _{i=0}^{j-1}{\beta }_i{\left(\mathbf{J}\right)}^i{\mathbf{r}}_0,$$

(12)

where the term ${\beta }_i$ minimizes the residual. We terminate the Krylov iteration using the following inexact Newton criteria [40]

$$\parallel J\delta c_j{+\mathbf{F}\left(c\right)\parallel }_2<\gamma {\parallel \mathbf{F}\left(c\right)\parallel }_2.$$

(13)

Here, parameter $\gamma$ is chosen to determine the convergence of the linear solver at each Newton iteration (we usually set $\gamma ={10}^{-3}$). We note that one main advantage of the Krylov subspace method that we use here is that it does not need the formation of the Jacobian matrix. It only requires the matrix-vector multiplication, $\mathbf{J}\mathbf{v}$, where $\mathbf{v}\in \{{\mathbf{r}}_0,\mathbf{J}{\mathbf{r}}_0,{\mathbf{J}}^2{\mathbf{r}}_0,\cdots \}$. This is why it is simply referred to as the Jacobian-free implementation, in which the operation of the Jacobian matrix times a vector can be estimated by

$$\mathbf{J}\mathbf{v}\approx \frac{\mathbf{F}(c+ϵ\mathbf{v})-\mathbf{F}\left(c\right)}{ϵ},$$

(14)

where $ϵ=\frac{1}{{n\parallel \mathbf{v}\parallel }_2}{\sum }_{i=1}^{n}b\left|c_i\right|+b$, n is the dimension of the linear system, and b is a constant whose magnitude is within a few orders of magnitude of the square root of machine round-off (it is typically set to ${10}^{-6}$ for 64-bit double precision).

Remark 1.

We remark that we do not employ any global convergence strategy such as the trust region method or line search method to improve the convergence of the nonlinear block [21,41] in this research. Typically, the line search or trust region of globalization methods are employed if the initial guess is far from the solution or the converged solution is not physically the correct one. The main reason why we do not use such methods is that, in both methods, one needs to calculate the direction search vectors which require the knowledge of a Jacobian matrix (exact or approximate), but in our matrix’s free implementation, we do not calculate the Jacobian matrix by any means. The second reason is that we are solving an initial value problem and, at the beginning of each nonlinear iteration, the previous time step solutions become a very good initial guess for the next set of solutions. With regard to the second point, we also note that, despite the fact they are time implicit, our integration techniques do not use arbitrarily large time steps for accuracy reasons. Furthermore, this prevents us drifting far away from the next time level solutions.

As we pointed out earlier, we utilize preconditioning techniques to accelerate the convergence of the JFNK’s Krylov step. Usually, the original linear system (9) is modified when a preconditioner is in use (the Jacobian matrix is multiplied by a preconditioning matrix from the left or right). When the linear system is multiplied by a preconditioning matrix from the left, the resulting method is called the left preconditioner (or the right preconditioner if the multiplication is performed from the right). We note that the left preconditioner also modifies the right-hand side of the linear system (9), meaning that the modified nonlinear residuals must be recovered. This brings additional coding complexity and computational costs [20,21,42]. On the other hand, the right preconditioner does not affect the right-hand side term and does not possess the aforementioned difficulties. Therefore, we prefer the right preconditioning (RP) strategy in this paper.

The RP procedure modifies and improves the linear system (9) as follows

$$\mathbf{J}{P}^{-1}\left(P\delta c\right)\approx -\mathbf{F}\left(c\right),$$

(15)

where P is the preconditioning matrix. The action of the right preconditioning can be defined similarly to (14), e.g.,

$$\mathbf{J}{P}^{-1}\mathbf{v}\approx \frac{\mathbf{F}(c+ϵP^{-1}\mathbf{v})-\mathbf{F}\left(c\right)}{ϵ},$$

(16)

where $\mathbf{v}$ is the Krylov vector. The execution of this procedure involves the following steps:

Step-1: 

Given $\mathbf{v}$, solve $Pz=\mathbf{v}$ for z (e.g., $z=P^{-1}\mathbf{v}$)

Step-2: 

Perform $\mathbf{J}z\approx \left(\mathbf{F}\right(c+ϵz)-\mathbf{F}(c\left)\right)/ϵ$,

where z is the preconditioned Krylov vector. Below, we describe two different ways of obtaining $z^{\prime }s$.

3.4. The Physics-Based Preconditioners

3.4.1. Based on the Backward Euler Method

We follow the steps of Kadioglu et al. [31,34] to derive the preconditioning equation. First, we recast the discretization (6) as

$$\frac{c_{i,j}^{k+1}-c_{i,j}^n}{\Delta t}-{\left[{\left(\left(D(x,y)c_x\right)\right)}_x\right]}_{i,j}^{k+1}-{\left[{\left(\left(D(x,y)c_y\right)\right)}_y\right]}_{i,j}^{k+1}\approx f_{i,j}^k,$$

(17)

where k represents the kth nonlinear iteration. Adding and subtracting $c_{i,j}^k$ to the time part and ${\left[{\left(\left(D(x,y)c_x\right)\right)}_x\right]}_{i,j}^k$, ${[{\left(\left(D(x,y)c_y\right)\right)}_y]}_{i,j}^k$ to both sides of Equation (17), we obtain

$$\begin{array}{c}\hfill \frac{\delta c_{i,j}}{\Delta t}-{[{\left(\left(D(x,y)c_x\right)\right)}_x]}_{i,j}^{k+1}+{[{\left(\left(D\left(x\right)c_x\right)\right)}_x]}_{i,j}^k-{[{\left(\left(D(x,y)c_y\right)\right)}_y]}_{i,j}^{k+1}+{[{\left(\left(D\left(x\right)c_y\right)\right)}_y]}_{i,j}^k\\ \hfill \approx -\frac{c_{i,j}^k-c_{i,j}^n}{\Delta t}+{[{\left(\left(D(x,y)c_x\right)\right)}_x]}_{i,j}^k+{[{\left(\left(D(x,y)c_y\right)\right)}_y]}_{i,j}^k+f_{i,j}^k,\end{array}$$

(18)

where $\delta c_{i,j}=c_{i,j}^{k+1}-c_{i,j}^k$. Assuming the equal diffusion coefficients at both iterates, we have

$$\begin{array}{c}\hfill \delta c_{i,j}-\frac{\Delta t}{{\Delta x}^2}{D}_{i+\frac{1}{2},j}(\delta c_{i+1,j}-\delta c_{i,j})+\frac{\Delta t}{{\Delta x}^2}{D}_{i-\frac{1}{2},j}(\delta c_{i,j}-\delta c_{i-1,j})\\ \hfill -\frac{\Delta t}{{\Delta y}^2}{D}_{i,j+\frac{1}{2}}(\delta c_{i,j+1}-\delta c_{i,j})+\frac{\Delta t}{{\Delta y}^2}{D}_{i,j-\frac{1}{2}}(\delta c_{i,j}-\delta c_{i,j-1})\approx -\mathrm{res}{c}_{i,j},\end{array}$$

(19)

where $\mathrm{res}{c}_{i,j}$ represents the nonlinear residual term at the $k^{th}$ iteration,

$$\mathrm{res}{c}_{i,j}\approx c_{i,j}^k-c_{i,j}^n-\Delta t{\left[{\left(\left(D(x,y)c_x\right)\right)}_x\right]}_{i,j}^k-\Delta t{[{\left(\left(D(x,y)c_y\right)\right)}_y]}_{i,j}^k-\Delta t{f}_{i,j}^k.$$

(20)

Combining the common terms yields

$$-c\delta c_{i-1,j}-b\delta c_{i+1,j}+a\delta c_{i,j}-d\delta c_{i,j+1}-e\delta c_{i,j-1}\approx -\mathrm{res}{c}_{i,j},$$

(21)

where

$$a=1+\frac{\Delta t}{{\Delta x}^2}{D}_{i+\frac{1}{2},j}+\frac{\Delta t}{{\Delta x}^2}{D}_{i-\frac{1}{2},j}+\frac{\Delta t}{{\Delta y}^2}{D}_{i,j+\frac{1}{2}}+\frac{\Delta t}{{\Delta y}^2}{D}_{i,j-\frac{1}{2}},$$

$$b=\frac{\Delta t}{{\Delta x}^2}{D}_{i+\frac{1}{2},j},c=\frac{\Delta t}{{\Delta x}^2}{D}_{i-\frac{1}{2},j},d=\frac{\Delta t}{{\Delta y}^2}{D}_{i,j+\frac{1}{2}},e=\frac{\Delta t}{{\Delta y}^2}{D}_{i,j-\frac{1}{2}}.$$

for $i=1,\cdots ,M_x-1$ and $j=1,\cdots ,M_y-1$.

3.4.2. Based on the Crank–Nicolson Method

To obtain the preconditioning equation, we recast the discretization (8) as

$$\begin{array}{c}\hfill \frac{c_{i,j}^{k+1}-c_{i,j}^n}{\Delta t}-\frac{{[{((D(x,y)c_x))}_x]}_{i,j}^{k+1}+{[{((D(x,y)c_x))}_x]}_{i,j}^n}{2}\\ \hfill -\frac{{[{((D(x,y)c_y))}_y]}_{i,j}^{k+1}+{[{((D(x,y)c_y))}_y]}_{i,j}^n}{2}\approx \frac{f_{i,j}^k+f_{i,j}^n}{2}.\end{array}$$

(22)

If we follow the similar steps as in Section 3.4.1, we obtain the following update formula

$$-c\delta c_{i-1,j}-b\delta c_{i+1,j}+a\delta c_{i,j}-d\delta c_{i,j+1}-e\delta c_{i,j-1}\approx -\mathrm{res}{c}_{i,j},$$

(23)

where

$$\begin{array}{ccc}\hfill \mathrm{res}{c}_{i,j}& \approx & c_{i,j}^k-c_{i,j}^n-\frac{\Delta t}{2}({[{((D(x,y)c_x))}_x]}_{i,j}^k+{[{((D(x,y)c_x))}_x]}_{i,j}^n)\hfill \\ & & -\frac{\Delta t}{2}({[{((D(x,y)c_y))}_y]}_{i,j}^k+{[{((D(x,y)c_y))}_y]}_{i,j}^n)-\frac{\Delta t}{2}(f_{i,j}^k+f_{i,j}^n),\hfill \end{array}$$

(24)

and

$$a=1+\frac{\Delta t}{{2\Delta x}^2}{D}_{i+\frac{1}{2},j}+\frac{\Delta t}{{2\Delta x}^2}{D}_{i-\frac{1}{2},j}+\frac{\Delta t}{{2\Delta y}^2}{D}_{i,j+\frac{1}{2}}+\frac{\Delta t}{{2\Delta y}^2}{D}_{i,j-\frac{1}{2}},$$

$$b=\frac{\Delta t}{{2\Delta x}^2}{D}_{i+\frac{1}{2},j},c=\frac{\Delta t}{{2\Delta x}^2}{D}_{i-\frac{1}{2},j},d=\frac{\Delta t}{{2\Delta y}^2}{D}_{i,j+\frac{1}{2}},e=\frac{\Delta t}{{2\Delta y}^2}{D}_{i,j-\frac{1}{2}}.$$

for $i=1,\cdots ,M_x-1$ and $j=1,\cdots ,M_y-1$.

Notice that we need boundary conditions when solving Equations (21) and (23) for $\delta c$. Essentially, we define few ghost cell values similarly to those described in Section 3.1 as $\delta c_{-1,j}\approx \delta c_{0,j}$ and $\delta c_{M_x,j}\approx \delta c_{M_x-1,j}$ for $j=1,\cdots ,M_y-1$ and $\delta c_{i,-1}\approx \delta c_{i,0}$ and $\delta c_{i,M_y}\approx \delta c_{i,M_y-1}$ for $i=1,\cdots ,M_x-1$. The update Equations (21) and (23) are solved iteratively to produce the preconditioned Krylov vector z.

The iterative procedure is based on a multi-grid V-cycle method specifically designed for our applications. A typical multi-grid V-cycle method consists of few recursively applied Gauss–Seidel (GS) or other basic such as Jacobi iterations from finest to coarsest and from coarsest to finest grid levels [41,43]. It is well known that the GS or other basic iterative methods (sometimes referred to as smoothers) are very good at smoothing the high frequency errors in just a few iterations. However, they are extremely slow to produce fully converged solutions. Often, the iteration count for convergence is related to the problem size (e.g., $M_x\times M_y$ in our 2D cases) which is impractical. We take advantage of the GS smoother to design our multi-grid V-cycle in the following way. At the beginning of each V-cycle, we perform two Gauss–Seidel iterations to remove or smooth out the high-frequency errors occurring at the finest grid. Then, we restrict ourselves to a coarser grid at which we form coarse grid solutions and the related coarse residuals to apply two more GS sweeps. We note that, at this coarse level, the previously low-frequency errors act like high-frequency ones and the application of a small number of GS sweeps are sufficient to smooth them out. This continues recursively until we reach the coarsest level (two grid cells in each direction) at which we perform hundred GS iterations. At the end of the coarsest level iterations, we recursively apply the prolongation step to interpolate the coarse level solutions to the finer levels at which we carry out two more GS sweeps. The completion of these steps defines one multi-grid V-cycle. In our test problems, we carry out five V-cycles that we found sufficient to decrease the GMRES iteration count to very small numbers. However, one can always increase or decrease the number of multi-grid V-cycles depending on the physical model and its numerical discretization.

Remark 2.

We remark that the preconditioning procedure that we are using in this paper is referred to as the physics-based or partial differential equation (PDE)-based method [20,22,29,41]. A different physical process usually corresponds to a different class of PDE’s. For example, the diffusion process is modeled by a parabolic type of PDE or a wave phenomenon described by a hyperbolic type of PDE. Some physical process, fluid flow for instance, can have a mixture of hyperbolic and parabolic characteristics. More importantly, some physical process can have multiple time-scale behavior resulting in less or more stiff terms in the governing equations. Depending on the stiffness measure, there exist well-established numerical techniques such as explicit algorithms (for non-stiff problems), implicit algorithms for stiff problems, or a semi-implicit algorithm for the mixed equation systems. We also note that stiff parts in the system are the source of the stringent time step stability restrictions. Therefore, having some physical insight regarding the model helps target specific stiff parts (diffusion terms in our case) to step over fine time scales resulting in the accelerated convergence of the linear algebra problem (Krylov step in our case). We remark that there exist matrix-based preconditioning technologies such as LU decomposition, incomplete LU and its variants, or other sorts of approaches [25,29]. They can be successfully implemented to accelerate the convergence. However, the main disadvantage of almost all of them is that they all require some information about the Jacobian matrix of the nonlinear system. As noted earlier, forming the Jacobian matrix can be complicated and storing it can be expensive. Our matrix-free method (JFNK) does not require the formation or storage of the Jacobian matrix. Our preferred physics-based preconditioner also has the same characteristics and avoids the explicit formation of the preconditioning matrix and its inverse, but it rather computes the effect of the inverse of the preconditioning matrix multiplied by a Krylov vector. This basic structure of our preconditioning procedure provides a unique opportunity for flexible and straight forward code executions. The reliability and efficiency of the preconditioners are demonstrated in the following section.

4. Results

4.1. 1D Case

4.1.1. Convergence Analysis

The theoretical order of accuracy of the numerical schemes that we use here is first -order in time and second order in space for the backward Euler and second order in time and second order in space for the Crank–Nicolson method. Of course, these theoretical orders exist for the linear models and linear methods. Whenever the model becomes nonlinear or sharp structures appear in the solution, it may become difficult to verify the expected order of accuracy. We have both situations in our model and solution profiles. Therefore, it makes sense to check whether we lose numerical accuracy in time or space. To do this, we will perform temporal and spatial convergence analysis, that subsequently described. Although there exists analytical solutions to the Fisher–Kolmogorov equation in certain situations [44], we prefer computational means to carry the convergence analysis. This is mainly because our numerical analysis can be easily applied to more general cases such as a coupled nonlinear fluid dynamical system with more physical effects. The time convergence analysis was performed in the following manner. We fix the number of grid points (e.g., $M=2048$) and run the code with sequentially refined time steps up to the final time $t_{final}=10$. For instance, we obtain the following approximate solution set $c^{\Delta t},c^{\frac{\Delta t}{2}},c^{\frac{\Delta t}{4}},\dots ,c^{\frac{\Delta t}{2^k}}$ where $k=0,1,\dots ,8$ and $\Delta t={10}^{-1}$. Then, we measure the $L_p$ norms of errors between two consecutive time step solutions, e.g.,

$${\left|\right|c^{\Delta t}-c^{\frac{\Delta t}{2}}\left|\right|}_p,{\left|\right|c^{\frac{\Delta t}{2}}-c^{\frac{\Delta t}{4}}\left|\right|}_p,\dots ,$$

where $p\in \{1,2,\infty \}$. These time errors can be superimposed (in log scale) on top of the reference first- and second-order slope indicator lines to provide insights into the convergence behavior. Alternatively, one can define the time rate of convergence by

$$r=ln\left(\left|\right|E(\Delta t/2^{k-1})\left|\right|/\left|\right|E(\Delta t/2^k)\left|\right|\right)/ln2,$$

(25)

where

$$E(\Delta t/2^{k-1})={\left|\right|c^{\Delta t/2^{k-1}}-c^{\Delta t/2^k}\left|\right|}_{p}k=1,\dots ,8.$$

Based on the second description, we made convergence tables to demonstrate and verify the theoretical time order of accuracy of the numerical schemes.

Table 1. Crank–Nicolson, $M=2048,\Delta t={10}^{-1}$.

Time Convergence Analysis with $\mathit{E}={\left|\right|{\mathit{c}}^{\Delta \mathit{t}/2^{\mathit{k}-1}}-{\mathit{c}}^{\Delta \mathit{t}/2^{\mathit{k}}}\left|\right|}_{\mathit{p}}$ in Columns 1, 3, 5
$\mathit{p}=\mathbf{1},\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate ($\mathit{r}$ in (25))	$\mathit{p}=\mathbf{2},\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate ($\mathit{r}$ in (25))	$\mathit{p}=\infty ,\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate ($\mathit{r}$ in (25))
$7.22\times {10}^{-5}$	2.0000	$4.28\times {10}^{-5}$	2.0004	$1.38\times {10}^{-4}$	2.0005
$5.55\times {10}^{-6}$	2.0000	$1.07\times {10}^{-5}$	2.0001	$3.44\times {10}^{-5}$	2.0001
$1.39\times {10}^{-6}$	2.0029	$2.67\times {10}^{-6}$	1.9985	$8.59\times {10}^{-6}$	1.9973
$3.46\times {10}^{-7}$	1.9980	$6.69\times {10}^{-7}$	1.9992	$2.15\times {10}^{-6}$	1.9998
$8.67\times {10}^{-8}$	1.9711	$1.67\times {10}^{-7}$	2.0105	$5.38\times {10}^{-7}$	2.0172
$2.21\times {10}^{-8}$	1.9506	$4.15\times {10}^{-8}$	2.0359	$1.33\times {10}^{-7}$	2.0687
$5.72\times {10}^{-9}$	1.9014	$1.01\times {10}^{-8}$	2.0615	$3.17\times {10}^{-8}$	2.1685
$1.53\times {10}^{-9}$		$2.42\times {10}^{-9}$		$7.05\times {10}^{-9}$

illustrates the time behavior of the Crank–Nicolson method using different norms. The convergence rates slightly change depending on the norm, however, they are all very close to the theoretical value which is two. In

Table 2. Backward Euler, $M=2048,\Delta t={10}^{-1}$.

Time Convergence Analysis with $\mathit{E}={\left|\right|{\mathit{c}}^{\Delta \mathit{t}/2^{\mathit{k}-1}}-{\mathit{c}}^{\Delta \mathit{t}/2^{\mathit{k}}}\left|\right|}_{\mathit{p}}$ in Columns 1, 3, 5
$\mathit{p}=\mathbf{1},\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate ($\mathit{r}$ in (25))	$\mathit{p}=\mathbf{2},\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate ($\mathit{r}$ in (25))	$\mathit{p}=\infty ,\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate ($\mathit{r}$ in (25))
$1.34\times {10}^{-3}$	0.9969	$2.17\times {10}^{-3}$	0.9984	$6.56\times {10}^{-3}$	0.9996
$6.72\times {10}^{-4}$	0.9983	$1.08\times {10}^{-3}$	0.9989	$3.28\times {10}^{-3}$	0.9989
$3.36\times {10}^{-4}$	0.9991	$5.43\times {10}^{-4}$	0.9993	$1.64\times {10}^{-3}$	0.9991
$1.68\times {10}^{-4}$	0.9995	$2.71\times {10}^{-4}$	0.9996	$8.21\times {10}^{-4}$	0.9995
$8.41\times {10}^{-5}$	0.9998	$1.36\times {10}^{-4}$	0.9998	$4.11\times {10}^{-4}$	0.9997
$4.21\times {10}^{-5}$	0.9999	$6.79\times {10}^{-5}$	0.9999	$2.05\times {10}^{-4}$	0.9998
$2.10\times {10}^{-5}$	0.9999	$3.39\times {10}^{-5}$	0.9999	$1.03\times {10}^{-4}$	0.9999
$1.05\times {10}^{-5}$		$1.70\times {10}^{-5}$		$5.14\times {10}^{-5}$

, we verify the time convergence behavior of the backward Euler method. Naturally, time errors are bigger compared to the Crank–Nicolson results, but again, we have a perfect match between the theoretical and numerical rate of convergence. This verification analysis is very important because it not only means that we have a reliable code implementation, but it also numerically proves that the JFNK method consistently converges all nonlinearities in the physical system as intended.

To numerically demonstrate the spatial convergence, we perform nine runs using $\frac{\Delta x}{2^k}$$(k=0,1,\dots ,8)$ and a fixed fine time step $\Delta t={10}^{-5}$ where $\Delta x=\frac{50}{2^6}$. In these runs, the time scheme is set to be the CN method that, by construction, contributes smaller time errors to the truncation term. With these settings and similar error norm measures described above, we formed

Table 3. Spatial errors and Conv. rate between sequentially refined grid solutions with $\Delta t={10}^{-5}$.

Spatial Convergence Analysis with $\mathit{E}={\left|\right|{\mathit{c}}^{\Delta \mathit{x}/2^{\mathit{k}-1}}-{\mathit{c}}^{\Delta \mathit{x}/2^{\mathit{k}}}\left|\right|}_{\mathit{p}}$ in Columns 1, 3, 5
$\mathit{p}=\mathbf{1},\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate	$\mathit{p}=\mathbf{2},\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate	$\mathit{p}=\infty ,\mathit{k}=\mathbf{1},\cdots ,\mathbf{8}$	Conv. Rate
$6.67\times {10}^{-4}$	1.9975	$1.93\times {10}^{-3}$	1.9987	$9.70\times {10}^{-3}$	1.9224
$1.67\times {10}^{-4}$	2.0006	$4.83\times {10}^{-4}$	1.9985	$2.56\times {10}^{-3}$	1.9972
$4.18\times {10}^{-5}$	2.0001	$1.21\times {10}^{-4}$	1.9999	$6.41\times {10}^{-4}$	1.9994
$1.04\times {10}^{-5}$	2.0002	$3.02\times {10}^{-5}$	1.9950	$1.60\times {10}^{-4}$	1.9847
$2.61\times {10}^{-6}$	2.0015	$7.57\times {10}^{-6}$	2.0064	$4.05\times {10}^{-5}$	2.0220
$6.52\times {10}^{-7}$	1.9639	$1.88\times {10}^{-6}$	2.0005	$9.97\times {10}^{-6}$	2.0647
$1.67\times {10}^{-7}$	1.9678	$4.71\times {10}^{-7}$	2.0402	$2.38\times {10}^{-6}$	2.2178
$4.27\times {10}^{-8}$		$1.14\times {10}^{-7}$		$5.12\times {10}^{-7}$

at the final time $t_{final}=10$. Table 3 clearly shows second-order spatial convergence. Again, this is a good indication of reliable and accurate code implementation. This also means that we can readily extend our work to multi-dimensional and more complex models. Here, we note that one could fix $\Delta t/\Delta x^2$ or $\Delta t^2/\Delta x^2$ and perform only one convergence analysis (space–time) for each numerical scheme to check the possible accuracy loss. However, in this way, it may be difficult to identify the source of inaccuracies (e.g., it could be related to the time or space terms in the truncation functional). Therefore, performing separate analysis as we did here is important to minimize the influence of the time error to the spatial one.

Figure 3 shows the time simulation of the tumor concentration of glioma cells from the initial to the final time $t_{final}=1000$ (days). We have two sets of figures. The figures in the left column correspond to the coarse grid ($M=128$) and the figures on the right column belongs to the fine grid ($M=4096$) runs. All runs are executed using the Crank–Nicolson time scheme with the same time step $\Delta t=0.01$. The simulation starts with a high concentration of glioma cells in a very narrow region. The concentration diffuses for a while (approximately 50 days), then it starts ramping up and spreading at the same time. The sub-figures show that the tumor cells spread the entire region in approximately 75 days. The last row of the figure indicates that the tumor concentration slowly reaches the steady state in approximately 1000 days. We see the effects of the discontinuous diffusion coefficient (e.g., the symmetric kinks in the solution profiles in the middle row). Clearly, the tumor cells in the high diffusion zone spread much faster than those in the low zone. We also notice the effects of the grid refinement. Fine grid runs (right column) capture better (more accurate) solution profiles that can be very significant in the sense that the amount of medical treatment (overdose/underdose) can depend on this.

In Figure 4, we compare our computational results with the results of [19] for a series of time runs. Both results are using $M=100$ grid points and fixed time steps $\Delta t={10}^{-2}$. We can see that the results are very similar. This is expected because both numerical schemes are second-order accurate in space and time. However, as far as efficiency is concerned, our algorithm is much faster to converge with two Newton and one GMRES iteration per Newton step for the aforementioned grid and time step options.

4.1.2. Performance Study of the Preconditioners

In this section, we include an efficiency study to monitor the performance of the new preconditioners. We ran the code with a fixed $\Delta t={10}^{-1}$ and $M=2048$ grid points until the final time $t_{final}=5.0$. As we pointed out earlier, the time accuracy of the numerical scheme affects the convergence speed of the JFNK method. In other words, more accurate time schemes result in better conditioned systems (the condition number is improved). This automatically accelerates the convergence of the linear solver (Krylov block) and consequently may decrease the number of nonlinear iterations (Newton block). Therefore, we expect that the JFNK method based on the Crank–Nicolson and the backward Euler time schemes behave differently. Indeed, we observe this behavior in Figure 5. Without the preconditioners, the average numbers of nonlinear and linear iterations are significantly different for the two time schemes. In fact, the JFNK method with the backward Euler scheme takes almost twice the number linear and nonlinear iterations to achieve the same final time run. A truly meaningful acceleration manifests when we turn the preconditioners on (Figure 5). The average number of Newton iterations are from four to six, which are very reasonable, but the average number of GMRES iterations (linear iterations) significantly drops whenever we activate the preconditioners. Notice that when the time progresses, the sharp structure smears and the solution profile becomes smoother (Figure 3). Then, the JFNK method with the preconditioners converges effortlessly, requiring only one GMRES iteration (second raw of Figure 5).

Another important aspect of our algorithm is the convergence rate of the nonlinear block. In the literature, it is reported that the JFNK method has the super-linearly convergence rate property [20,41]. We computationally verify this behavior in Figure 6. Figure 6 shows the convergence rate as a function of nonlinear iterations. To produce this figure, we use $M=2048$ grid points and a fixed time step $\Delta t={10}^{-1}$. The code runs until the final time $t_{final}=1.0$ with the second-order version (Crank–Nicolson-based) of the preconditioner (Section 3.4.2). Figure 6 clearly indicates that the rate of convergence of our algorithm with the new preconditioner is super linear.

4.2. Two-Dimensional Case

In Section 4.1, by solving one-dimensional test cases, we were able to demonstrate the accuracy and the effectiveness of the numerical schemes. We expect a similar performance from the algorithm when it is applied to multi-dimensional and more complicated models. This is essentially because any multi-dimensional problem can be put into a nonlinear equation form with an increasing number of components in the unknown vector. However, this is not an issue for the JFNK method, because it is known to self-consistently guarantee the convergence of all nonlinearities in the system. Indeed, just to support the above discussion, we performed a time convergence analysis by solving the two-dimensional model with the initial and boundary settings described in Section 2.

Table 4. 2D: Crank–Nicolson, $M_x=M_y=2048,\Delta t={10}^{-1}$.

Time Convergence Analysis with $\mathit{E}={\left|\right|{\mathit{c}}^{\Delta \mathit{t}/2^{\mathit{k}-1}}-{\mathit{c}}^{\Delta \mathit{t}/2^{\mathit{k}}}\left|\right|}_{\mathit{p}}$ in Columns 1, 3, 5
$\mathit{p}=\mathbf{1},\mathit{k}=\mathbf{1},\cdots ,\mathbf{5}$	Conv. Rate	$\mathit{p}=\mathbf{1},\mathit{k}=\mathbf{1},\cdots ,\mathbf{5}$	Conv. Rate	$\mathit{p}=\infty ,\mathit{k}=\mathbf{1},\cdots ,\mathbf{5}$	Conv. Rate
$2.04\times {10}^{-3}$	2.0007	$2.66\times {10}^{-4}$	2.0008	$1.10\times {10}^{-4}$	2.0001
$5.09\times {10}^{-4}$	2.0003	$6.64\times {10}^{-5}$	2.0002	$2.76\times {10}^{-5}$	2.0000
$1.27\times {10}^{-4}$	2.0057	$1.66\times {10}^{-5}$	2.0026	$6.90\times {10}^{-6}$	2.0030
$3.17\times {10}^{-5}$	2.0087	$4.14\times {10}^{-6}$	2.0071	$1.72\times {10}^{-6}$	2.0079
$7.87\times {10}^{-6}$		$1.03\times {10}^{-6}$		$4.28\times {10}^{-7}$

clearly verifies the second-order convergence. Figure 7 illustrates the two-dimensional performance study by employing the second-order space and time scheme. Sub-figures show the average number of Newton and GMRES iterations in time with and without using the preconditioner. With the preconditioner, the method runs extremely efficiently, needing only a few Newton and GMRES steps. This remarkable performance is due to the newly devised physics-based preconditioners and their internal multi-grid mechanism whose convergence is independent of the problem size. Here, we would like to note that we never intend to measure the CPU time performance of our algorithms. This is because such a study can depend on the computer technology and the code implementation that can vary from person to person.

Figure 8 shows the time history of the glioma cell invasion in a two-dimensional setting. Again, we start with a very high concentration of glioma cells (approximately 40 cells) in a very narrow region at the center of the computational domain ($50\times 50$ mm). The tumor concentration rapidly spreads into the high diffusion zones in approximately 50 days (third row in Figure 8). Upon hitting the low-diffusion zones, the spreading slows down for a while but the concentration starts building up and eventually reaches the steady state value (40 cells). Sub-figures in Figure 8 also indicate the sharp capturing of the tumor profiles with our second-order time and space scheme.

5. Conclusions

In this paper, we introduced a Jacobian-free Newton–Krylov method with very effective preconditioning strategies to study tumor growth dynamics. The JFNK method does not require the formation, storage, and inversion of the Jacobian matrix that make it very unique. We utilized two physics-based preconditioners which are associated with the time integration schemes to accelerate the JFNK’s convergence. Our preconditioners rely on physical models rather than prescribed matrices. This provides more practical applications and flexible coding. We solved test problems to demonstrate the numerical convergence (both in time and space) of the proposed schemes. Theoretical temporal and spatial accuracies were verified through tables. We performed efficiency studies in one- and two-dimensional settings and provided figures that illustrate the average number of linear and nonlinear iterations with and without using the preconditioners. Our study clearly indicates that the JFNK method with the new preconditioners works very well. The average number of GMRES iterations drops down to very low numbers. This is a very important conclusion because, when we apply the methodology to more complex tumor dynamics (coupled with fluid dynamics [11]), the efficiency will be crucial. We are currently working on more realistic coupled fluid dynamics and interface dynamics model with full radiotherapeutic and chemotherapeutic effects to simulate the tumor growth behavior and will report our findings in a separate journal paper.