## 3 Numerical Algorithm
### 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}(c)\approx 0\) over the sequence of a linear system defined by  
\[\mathbf{J}(c^{k})\delta c^{k} \approx -\mathbf{F}(c^{k}),\] \[c^{k+1} \approx c^{k}+\delta c^{k},\hskip 14.226378 ptk=0,1,\cdots \tag{9}\]  
where \(\mathbf{J}(c^{k})=\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  
\[\|\mathbf{F}(c^{k})\|_{2}<\text{tol}_{res}\|\mathbf{F}(c^{0})\|_{2}, \tag{10}\]  
where \(\text{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}(c)-\mathbf{J}\delta c_{0}. \tag{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 \(j\)th Krylov iteration minimizes the quantity \(\|\mathbf{J}\delta c_{j}+\mathbf{F}(c)\|_{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}(\mathbf{J})^{i} \mathbf{r}_{0}, \tag{12}\]where the term \(\beta_{i}\) minimizes the residual. We terminate the Krylov iteration using the following inexact Newton criteria [40]  
\[\|f\delta c_{j}+\mathbf{F}(c)\|_{2}<\gamma\|\mathbf{F}(c)\|_{2}. \tag{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{Jv}\), where \(\mathbf{v}\in\{\tau_{0},\mathbf{Jv}_{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{Jv}\approx\frac{\mathbf{F}(c+\epsilon\mathbf{v})-\mathbf{F}(c)}{ \epsilon}, \tag{14}\]  
where \(\epsilon=\frac{1}{n\|\mathbf{v}\|_{2}}\sum_{i=1}^{n}b|c_{i}|+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). 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.