An explicit Jacobian for Newton's method applied to nonlinear initial boundary value problems in summation-by-parts form

1.   Introduction

Nonlinear systems of partial differential equations are common in computational science and engineering, and present multiple challenges. Stability is needed for reliability and high accuracy for fine solution details. For fast turnaround and timely result delivery, generic systems of nonlinear equations from discretizations of the form

must be solved efficiently ^[1]. Several techniques exist to solve (1.1); for example, dual-time stepping ^[2,3], optimization algorithms ^[4], or iterative methods ^[1]. Among the classical iterative methods, Newton's method is an effective choice due to its quadratic convergence order. The obvious drawback with Newton's method is that the Jacobian must be known. Methods that bypass this requirement and instead approximate the Jacobian lead to lower convergence orders. A typical example is the secant method ^[1]. Alternatively, by using Newton-Krylov methodologies ^[5], only the action of the Jacobian is required and can be approximated by $J_{ \mathit{\boldsymbol{\mathcal{F}}}}(\mathit{\boldsymbol{\phi}}^k) \delta \mathit{\boldsymbol{u}} \approx (\mathit{\boldsymbol{\mathcal{F}}}(\mathit{\boldsymbol{\phi}}^k + \delta \mathit{\boldsymbol{u}}) - \mathit{\boldsymbol{\mathcal{F}}}(\mathit{\boldsymbol{\phi}}^k))/\delta$, where $\delta$ is small and $\mathit{\boldsymbol{u}}$ depends on the subspaces in the Krylov iterations. The advantage of Newton-Krylov methods is that an explicit Jacobian is never required, but sophisticated preconditioners becomes necessary ^[6] instead.

The focus in this paper is to facilitate the use of Newton's method where the key component is an exact explicit form of the Jacobian of (1.1). To exemplify our technique, we will use finite-difference operators in summation-by-parts (SBP) form ^[7,8] to discretize the incompressible Navier-Stokes (INS) equations in space. The boundary conditions will be weakly imposed via the simultaneous approximation term (SAT) technique ^[9]. In ^[10,11], a discretization based on the SBP-SAT technique of the nonlinear INS equations was proven to be stable, which is the key prerequisite.

Based on the formulation in ^[10], we show how the Jacobian can be explicitly calculated. It is also shown that the Jacobian has a block structure, where several blocks can be precomputed and reused when forming $\mathit{\boldsymbol{\mathcal{F}}}$, making the procedure very efficient.

To keep the paper focused on the derivation of the Jacobian, we follow ^[10] and consider a Cartesian grid. Exact Jacobians for numerical discretizations have recently been developed in ^[12] for so-called entropy stable numerical discretizations in SBP form in a periodic setting. Our new technique is not restricted to such specific discretizations, and we include the specific Jacobian related to the boundary conditions. The technique demonstrated in this paper can be used in a straightforward way on any numerical method for IBVPs that can be formulated in matrix-vector form. In addition, it can be readily extended to curvilinear and unstructured grids, arbitrary dimensions, and other sets of linear and nonlinear equations. In principle, all that is needed for the existence of the Jacobian is that $\mathit{\boldsymbol{\mathcal{F}}}$ is differentiable with respect to $\mathit{\boldsymbol{\phi}}$. This covers the various nonlinearities that arise in discretizations of the Navier-Stokes equations (both compressible and incompressible). However, the feasibility of explicitly deriving the Jacobian is highly dependent on the way in which $\mathit{\boldsymbol{\mathcal{F}}}$ is presented. As we shall see, discretizations in SBP-SAT form are particularly straightforward to differentiate, making Newton's method an attractive solution method.

The rest of the paper proceeds as follows. We introduce the continuous problem in Section 2 and present the semi-discrete formulation in Section 3. The Jacobian of the discretization is derived in Section 4. Implicit time integration is discussed in Section 5 and numerical experiments are performed in Section 6. Finally, conclusions are drawn in Section 7.

2.   Problem formulation

As an illustrative example of our technique, we consider the scenario illustrated in Figure 1. An incompressible fluid is moving from left to right. Hence, the left side is an inflow boundary, where Dirichlet conditions are imposed, and the right side is an outflow boundary, where natural boundary conditions ^[13] are imposed. The lower part of the domain is a no-slip wall and the upper side is an outflow boundary, where again natural conditions are imposed. The initial-boundary value problem for the INS equations that we consider is

In (2.1), $\vec{\mathit {w}} = (u, v, p)^\top$, where $u, v$ are the velocities in the $x, y$ direction, respectively, and $p$ is the pressure. Furthermore, $\Omega = [0, 1]^2$ is the domain and $\partial \Omega$ its boundary. The initial data $\vec{\mathit {f}}$ and boundary data $\vec{\mathit {g}}$ are sufficiently smooth and compatible functions and the spatial operator is given by ^[10]

The matrices in (2.1) and $(2.2)$ are

and $\epsilon > 0$ is the viscosity coefficient. To obtain a stable discretization with the SBP-SAT technique, the nonlinear convective terms in (2.1) are split into a skew-symmetric form in (2.2) by taking an average between the conservative and non-conservative formulations ^[10]. This formulation is possible due to the divergence constraint. Lastly, the explicit form of the boundary conditions $\mathcal{H} \vec{\mathit {w}} = \vec{\mathit {g}}$ reads

For completeness, we will show how to bound the solution. Only the south side of the domain is discussed explicitly for simplicity. Details of the upcoming analysis are found in ^[10].

For two vector functions $\vec{\mathit {\phi}}, \vec{\mathit {\psi}}$ defined on $\Omega$, we introduce the inner product and norm

By multiplying (2.1) by $2 \vec{\mathit {w}}^\top$ from the left and integrating over $\Omega$, we get

where $\nabla \vec{\mathit {w}} = (\nabla u, \nabla v, \nabla p)^\top$, $\|\nabla \vec{\mathit {w}}\|^2_{ \tilde{I}}$ is a dissipative volume term, and

contains the boundary terms evaluated at the south boundary. The other boundary terms are assumed dissipative and ignored. Imposing $u = v = 0$ results in $BT = 0$. Integrating (2.4) in time (assuming homogeneous boundary conditions on all sides) leads to

which bounds the semi-norm of the solution ($\| \vec{\mathit {w}}\|^2_{\tilde{I}}$) and its gradients ($\|\nabla \vec{\mathit {w}}\|^2_{\tilde{I}}$) for any time.

3.   The semi-discrete scheme

A brief introduction of the SBP-SAT technique is provided below, see ^[7,8] for extensive reviews. We discretize the domain $\Omega = [0, 1]^2$ with $N+1$ and $M+1$ grid points; $x_i = i/N$, $i = 0, \dots, N$ and $y_j = j/M$, $j = 0, \dots, M$ and let $n = (N+1)(M+1)$ denote the total number of grid points. A scalar function $q = q(x, y)$ defined on $\Omega$ is thereby represented on the grid by $\mathit{\boldsymbol{q}} = (q_{00}, \dots, q_{0M}, \dots q_{N0}, \dots q_{NM})^\top$ where $q_{ij} = q(x_i, y_j)$. For the vector-valued function $\vec{\mathit {w}} = (u, v, p)^\top$, the approximation is arranged as $\vec{\mathit {\mathit{\boldsymbol{w}}}} = (\mathit{\boldsymbol{u}}^\top, \mathit{\boldsymbol{v}}^\top, \mathit{\boldsymbol{p}}^\top)^\top$. Let $\mathit{\boldsymbol{D_x}} = (P^{-1}_x Q_x) \otimes I_{M+1}$ and $\mathit{\boldsymbol{D_y}} = I_{N+1} \otimes (P^{-1}_y Q_y)$, where $\otimes$ denotes the Kronecker product. Then the approximations of the spatial derivatives are given by

The matrices $P_{x, y}$ are diagonal and positive definite, so that $\mathit{\boldsymbol{P}} = P_x\otimes P_y$ forms a quadrature rule that defines the norm $\| \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{I_3\otimes \mathit{\boldsymbol{P}}} = \vec{\mathit {\mathit{\boldsymbol{w}}}}^\top (I_3\otimes \mathit{\boldsymbol{P}}) \vec{\mathit {\mathit{\boldsymbol{w}}}} \approx \iint_\Omega \vec{\mathit {w}}^\top \vec{\mathit {w}} d\Omega$. We have also introduced $I_{k}$, which is the identity matrix of size $k$. Moreover, the matrices $Q_{x, y}$ satisfy the SBP property

where $E_{0x, y} = \text{diag}(1, 0, 0, \dots, 0)$ and $E_{N, M} = \text{diag}(0, 0, 0, \dots, 1)$ are matrices of appropriate sizes.

By using the notation above, the semi-discrete approximation of (2.1) becomes ^[10]

The discrete spatial operator is given by

and the block matrices are

where $\mathit{\boldsymbol{U}}, \mathit{\boldsymbol{V}} \in \mathbb{R}^{n\times n}$ are diagonal matrices holding $\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{v}}$, respectively. The matrices $\mathit{\boldsymbol{I}}$ and ${\bf{0}}$ are the identity and the zero matrix of size $n\times n$. Furthermore, $\mathit{\boldsymbol{\mathcal{S}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ contains penalty terms that enforce the boundary conditions.

The purpose of the SAT $\mathit{\boldsymbol{\mathcal{S}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ is $i)$ to enforce the boundary conditions in (2.3) and $ii)$ to stabilize the solution. One penalty term for each of the boundary conditions in (2.3) will be constructed. Let $k\in\{W, E, S, N\}$. The SAT at boundary $k$ that enforces the boundary condition $H^k\vec{\mathit {w}} = \vec{\mathit {g}}$ has the general form

In (3.3), $\Sigma^k$ is the penalty matrix to be determined for stability at boundary $k$. The quadratures are

For the boundary conditions listed in (2.3), the penalty terms are

where

As an example, the south penalty term can be written as

which is a more convenient notation for the derivation of the Jacobian in Section 4.

We will show in the following section that this specific choice of penalty matrices leads to nonlinear stability. The total penalty term in (3.2) becomes

For completeness, we also show schematically how to obtain an energy estimate (again, all details can be found in ^[10]). Similarly to the continuous analysis, we omit all boundaries except for the south one. By mimicking the continuous path ^[14], we multiply (3.2) by $2 \vec{\mathit {\mathit{\boldsymbol{w}}}}^\top (I_3\otimes \mathit{\boldsymbol{P}})$ from the left and use the SBP property (3.1) to get

where $\|\nabla \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}} = (I_3\otimes \mathit{\boldsymbol{D_x}} \vec{\mathit {\mathit{\boldsymbol{w}}}})^\top (I_3\otimes \mathit{\boldsymbol{P}}) \mathit{\boldsymbol{\tilde{I}}}(I_3\otimes \mathit{\boldsymbol{D_x}} \vec{\mathit {\mathit{\boldsymbol{w}}}}) + (I_3\otimes \mathit{\boldsymbol{D_y}} \vec{\mathit {\mathit{\boldsymbol{w}}}})^\top (I_3\otimes \mathit{\boldsymbol{P}}) \mathit{\boldsymbol{\tilde{I}}}(I_3\otimes \mathit{\boldsymbol{D_y}} \vec{\mathit {\mathit{\boldsymbol{w}}}})$ is the dissipative volume term corresponding to the continuous one and

contains all terms evaluated at the boundary.

The semi-norm of the solution ($\| \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}}$) is bounded if the right-hand side of (3.7) is non-positive. By expanding (3.8) and using the explicit form of $\mathit{\boldsymbol{\mathcal{S}}}^S(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ stated in (3.4), we find

where term $I$ is obtained from the governing equation and term $II$ from the penalty term. As in the continuous setting, the boundary terms vanish. Integrating (3.7) in time (assuming homogeneous dissipative boundary conditions at all boundaries) leads to

which is the semi-discrete version of the estimate in (2.5).

4.   Exact computation of the Jacobian

In this section, we will explicitly compute the Jacobian of $\mathit{\boldsymbol{\mathcal{L}}}$ and $\mathit{\boldsymbol{\mathcal{S}}}$ in (3.2). Let $\mathit{\boldsymbol{h}} : \mathbb{R}^{n} \to \mathbb{R}^{n}$, where $n = (N+1)(M+1)$ is the total number of grid points, be a differentiable vector function. For a given vector $\mathit{\boldsymbol{u}} = (u_{00}, \dots u_{NM})^\top \in \mathbb{R}^{n}$, $\mathit{\boldsymbol{h}}$ outputs the vector $\mathit{\boldsymbol{h}}(\mathit{\boldsymbol{u}}) = (h_{00}, \dots h_{NM})^\top \in \mathbb{R}^{n}$. The Jacobian matrix $J_{ \mathit{\boldsymbol{h}}} \in \mathbb{R}^{n\times n}$ of $\mathit{\boldsymbol{h}}$ is given by

We will first derive the Jacobian of the different terms in $\mathit{\boldsymbol{\mathcal{L}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ and, at the end, add the terms and state the complete result. To start, consider the vector function

Since

the Jacobian of $\mathit{\boldsymbol{h}}(\mathit{\boldsymbol{u}}) = \mathit{\boldsymbol{u}}$ becomes $J_{ \mathit{\boldsymbol{u}}} = \mathit{\boldsymbol{I}}$. Now let

Then, in the same way

Thus, $J_{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}} = \mathit{\boldsymbol{D_x}}$.

To derive the Jacobian of the nonlinear term, $\mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}$, we let

By using the product rule, we get that

Hence,

where $\underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}} = \text{diag}(\mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})$.

In a similar manner, for

we get that

Hence, $J_{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}}} = 2 \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}}$. To summarize, we have shown that

4.1.   The Jacobian of the spatial operator

Having established these building blocks, we next consider the terms in $\mathit{\boldsymbol{\mathcal{L}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$. Since these terms have a block structure, so will their Jacobians. Let $\mathit{\boldsymbol{h}}^1, \mathit{\boldsymbol{h}}^2, \mathit{\boldsymbol{h}}^3: \mathbb{R}^{3n} \to \mathbb{R}^n$ be differentiable functions of $\vec{\mathit {\mathit{\boldsymbol{w}}}}$ and define $\tilde{ \mathit{\boldsymbol{h}}} : \mathbb{R}^{3n}\to \mathbb{R}^{3n}$ given by

Since

it follows that

and similarly for every element in $\mathit{\boldsymbol{h}}^1$. Therefore, the Jacobian of $\mathit{\boldsymbol{h}}^1$ can be expressed as

The same holds for $\mathit{\boldsymbol{h}}^2$ and $\mathit{\boldsymbol{h}}^3$. Thus, the Jacobian of $\tilde{ \mathit{\boldsymbol{h}}}$ is given by

and each block in $J_{\tilde{ \mathit{\boldsymbol{h}}}}$ is of size $n\times n$.

For the first term in $\mathit{\boldsymbol{\mathcal{L}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$, $\mathit{\boldsymbol{A}} (I_3 \otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}$, we get

The last identities are useful when deriving $J_{ \mathit{\boldsymbol{A}} (I_3 \otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}}$. By using (4.1), we get that

Thus,

Likewise for the second term in $\mathit{\boldsymbol{\mathcal{L}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$, $(I_3 \otimes \mathit{\boldsymbol{D_x}}) \mathit{\boldsymbol{A}} \vec{\mathit {\mathit{\boldsymbol{w}}}}$, note that

where we have used that $\mathit{\boldsymbol{V}} \mathit{\boldsymbol{u}} = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{v}}$. Hence,

The next two terms in $\mathit{\boldsymbol{\mathcal{L}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ are treated in a similar manner and we get that

Finally, the contribution to the Jacobian of the linear viscous terms simply becomes

Adding all of the terms proves the following proposition, which is the first of the two main results of this paper.

Proposition 1. The Jacobian $J_{ \mathit{\boldsymbol{\mathcal{L}}}}$ of the discrete operator $\mathit{\boldsymbol{\mathcal{L}}}$ in (3.2) is

where

4.2.   The Jacobian of the penalty terms

By following the procedure presented above, we next derive the Jacobian for $\mathit{\boldsymbol{\mathcal{S}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$. To start, we rewrite $\mathit{\boldsymbol{\mathcal{S}}}^S(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ as

The Jacobian of $\mathit{\boldsymbol{\mathcal{S}}}^S(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ is

The first block in $J_{ \mathit{\boldsymbol{\mathcal{S}}}^S}$ becomes

Since $\mathit{\boldsymbol{P}}^S$ is diagonal, we have $\mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^s \mathit{\boldsymbol{u}} = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}}$ and the second block is

Note that $\mathit{\boldsymbol{\mathcal{S}}}^S$ does not depend on $\mathit{\boldsymbol{p}}$ and also that $\mathit{\boldsymbol{\mathcal{S}}}_2^S$ and $\mathit{\boldsymbol{\mathcal{S}}}_3^S$ are both independent of $\mathit{\boldsymbol{u}}$. Hence, the remaining non-zero blocks of $J_{ \mathit{\boldsymbol{\mathcal{S}}}^S}$ are

where we have used that $\mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^s \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{P}}^s \mathit{\boldsymbol{V}} \mathit{\boldsymbol{v}}$. Therefore,

For non-homogeneous boundary conditions, the boundary data $\mathit{\boldsymbol{g}}$ will affect the Jacobian if the SATs are nonlinear with respect to $\vec{\mathit {\mathit{\boldsymbol{w}}}}$. We illustrate this by considering $\mathit{\boldsymbol{\mathcal{S}}}^W_1$ (the first block in $\mathit{\boldsymbol{\mathcal{S}}}^W$), which we rewrite in a similar manner as we did for $\mathit{\boldsymbol{\mathcal{S}}}^S_1$ in (4.3) and get

Note that the terms $- \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^W(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1)/2$ and $\epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^W(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1)$ are nonlinear and linear with respect to $\vec{\mathit {\mathit{\boldsymbol{w}}}}$ (via $\mathit{\boldsymbol{u}}$), respectively. The Jacobian to the linear term simply becomes

For the nonlinear term, we use that $\mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^W = \mathit{\boldsymbol{P}}^W \mathit{\boldsymbol{U}}$ and $\mathit{\boldsymbol{U}} \mathit{\boldsymbol{g}}_1 = \underline{{ \mathit{\boldsymbol{g}}_1}} \mathit{\boldsymbol{u}}$, which yield

Since $\mathit{\boldsymbol{\mathcal{S}}}_1^W$ is independent of both $\mathit{\boldsymbol{v}}$ and $\mathit{\boldsymbol{p}}$, its Jacobian becomes

The Jacobian of the other penalty terms are derived in a similar manner and we have therefore proved the second main result of this paper.

Proposition 2. The Jacobian of the total penalty term (3.6) is

where

Remark 1. We see from Proposition 4.1 and Proposition 2 that parts of the blocks in the Jacobian of both $J_{ \mathit{\boldsymbol{\mathcal{L}}}}$ and $J_{ \mathit{\boldsymbol{\mathcal{S}}}}$ are obtained directly from the construction of $\mathit{\boldsymbol{\mathcal{L}}}$. The few remaining parts are obtained by $i)$ matrix multiplications between a diagonal matrix and a non-diagonal one (for example, $\mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}}$) and $ii)$ matrix additions. This leads to few new additional operations and hence efficiency.

5.   The fully discrete scheme

To evolve the system (3.2) in time, we will, for simplicity and ease of explanation, use the implicit backward Euler method. More accurate and efficient methods could be used in the same manner in practice. For an ordinary differential system of equations of the form

where $\mathit{\boldsymbol{\phi}}$ is a function defined on the grid and $\mathcal{M}$ is a constant matrix, the backward Euler scheme becomes

In (5.1), $\Delta t$ is the size of the time step and the superindices $i$ and $i+1$ are the solution at time level $i$ and $i+1$, respectively.

In order to obtain $\mathit{\boldsymbol{\phi}}^{i+1}$, the system of nonlinear equations in (5.1) must be solved. One strategy is to first form the function in (1.1), which results in

If we find a vector $\mathit{\boldsymbol{\phi}}^*$ such that $\mathit{\boldsymbol{\mathcal{F}}}(\mathit{\boldsymbol{\phi}}^*) = 0$, then $\mathit{\boldsymbol{\phi}}^{i+1} = \mathit{\boldsymbol{\phi}}^*$. To solve (5.2), we employ Newton's method ^[1], which is described in Algorithm 1. This allows us to solve a sequence of linear systems of equations and arrive at an approximation of $\mathit{\boldsymbol{\phi}}^{i+1}$.

For the INS equations, $\mathit{\boldsymbol{\phi}} = \vec{\mathit {\mathit{\boldsymbol{w}}}}$, $\mathit{\boldsymbol{\mathcal{H}}}(\mathit{\boldsymbol{\phi}}) = \mathit{\boldsymbol{\mathcal{L}}}(\mathit{\boldsymbol{\phi}}) - \mathit{\boldsymbol{\mathcal{S}}}(\mathit{\boldsymbol{\phi}})$, and $\mathcal{M} = \mathit{\boldsymbol{\tilde{I}}}$. Hence, (5.2) becomes

Furthermore, $J_{ \mathit{\boldsymbol{\mathcal{H}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}}) = J_{ \mathit{\boldsymbol{\mathcal{L}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}}) - J_{ \mathit{\boldsymbol{\mathcal{S}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$, which yields

to be used in the Newton iterations. In (5.4), $J_{ \mathit{\boldsymbol{\mathcal{L}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ and $J_{ \mathit{\boldsymbol{\mathcal{S}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$ are given in Proposition 4.1 and Proposition 4.2, respectively.

6.   Numerical experiments

A simple finite-difference approximation of the Jacobian is given by ^[5]

The approximation in (6.1) was used during the implementation of the analytical expression of $J_{ \mathit{\boldsymbol{\mathcal{F}}}}$ since we expected $\|J - \hat{J}\|_{\infty}$ to be small. This allowed us to write unit tests ensuring that the Jacobian had been correctly implemented, by comparing it to the approximation. In (6.1), a small $\delta$ leads to a good approximation. However, note that if $\delta$ is chosen too small, the approximation will be contaminated by floating-point roundoff errors, which limits the practically achievable accuracy of $J$ ^[5].

Computing difference approximations of the Jacobian also allowed us to compare the efficiency of Newton's method using approximate versus analytical Jacobians. Note that computing the approximation (6.1) requires $n$ evaluations of $\mathit{\boldsymbol{\mathcal{F}}}$, resulting in $O(n^2)$ complexity, compared to the $O(n)$ complexity of evaluating the exact Jacobian. Table 1 shows the execution times for evaluating the analytical Jacobian versus computing the approximation (6.1) at increasing resolutions.

As expected, due to the large number of evaluations of $\mathit{\boldsymbol{\mathcal{F}}}$ needed to compute the approximation, such a strategy quickly becomes infeasible.

It is readily seen that the number of floating point operations needed to evaluate the discrete spatial operator $\mathit{\boldsymbol{\mathcal{L}}}$ grows linearly with the degrees of freedom $n$. Consider, for example, the term $\mathit{\boldsymbol{A}}(I_3 \otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}$. The first product, $(I_3 \otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}$, are finite difference approximations at each point in the grid, resulting in $Cn$ operations, where $C$ depends on the width of the difference stencil. The matrix $\mathit{\boldsymbol{A}}$ is a $3$-by-$3$ block matrix with diagonal blocks, and so results in another $O(n)$ number of operations. Analogously, the remaining terms in $\mathit{\boldsymbol{\mathcal{L}}}$ each contribute $O(n)$ operations. The arithmetic complexity of evaluating a penalty term $\mathit{\boldsymbol{\mathcal{S}}}$ is $O(\sqrt{n})$ (assuming equal resolution in the horizontal and vertical directions), since $\mathit{\boldsymbol{\mathcal{S}}}$ acts only on the grid boundary. Hence, the arithmetic complexity of evaluating $\mathit{\boldsymbol{\mathcal{F}}}$ is $O(n)$.

Let us study the arithmetic complexity of evaluating the Jacobian $J_{ \mathit{\boldsymbol{\mathcal{F}}}}$ of $\mathit{\boldsymbol{\mathcal{F}}}$. Inspecting the form of the Jacobian $J_{ \mathit{\boldsymbol{\mathcal{L}}}}$ in Proposition 4.1 we see a number of terms that need to be evaluated. The partial derivatives $\underline{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}$, $\underline{ \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{v}}}$, etc, have already been computed as part of the evaluation of $\mathit{\boldsymbol{\mathcal{L}}}$, and hence can be disregarded. Similarly, terms that do not depend on the solution, such as $\mathit{\boldsymbol{D_x}}$, $\mathit{\boldsymbol{D_x}}^2$, etc, can be disregarded since they remain constant throughout the simulation. Finally, we have terms of the type $\mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}}$, $\mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{V}}$, etc. These are all products of a diagonal matrix and a banded difference stencil matrix, and each contribute with $O(n)$ operations. Summing the terms uses $O(n)$ operations. Therefore, the arithmetic complexity of evaluating $J_{ \mathit{\boldsymbol{\mathcal{L}}}}$ is $O(n)$. In fact, the number of operations needed to evaluate products like $\mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}}$ or $\mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{V}}$ do not exceed the number of operations needed to compute the discrete partial derivatives involved in $\mathit{\boldsymbol{\mathcal{L}}}$. Hence, the cost ratio of evaluating $J_{ \mathit{\boldsymbol{\mathcal{L}}}}$ and evaluating $\mathit{\boldsymbol{\mathcal{L}}}$ is less than $1$ (i.e., the additional cost of evaluating $J_{ \mathit{\boldsymbol{\mathcal{L}}}}$ is small). As before, the arithmetic complexity of evaluating the Jacobian with respect to a boundary penalty $\mathit{\boldsymbol{\mathcal{S}}}$ is $O(\sqrt{n})$ since it acts only on the boundary of the grid. Thus, the total arithmetic complexity of evaluating $J_{ \mathit{\boldsymbol{\mathcal{F}}}}$ is less than the cost of evaluating $\mathit{\boldsymbol{\mathcal{F}}}$.

6.1.   The order of accuracy

The method of manufactured solution ^[15] is used to verify the implementation. In all computations in this subsection, the initial guess is the solution from the previous time step and the tolerance $tol$ in Algorithm 1 is set to $10^{-12}$. For the SBP operators SBP21 and SBP42, the expected orders of accuracy for the system (3.2) are 2 and 3, respectively ^[16,17].

Remark 2. The SBP21 and SBP42 operators (in general, SBPpq) are difference operators that satisfy the SBP property (3.1) and all derivations in this paper hold for these difference operators. The numbers pq correspond to the interior accuracy and accuracy on the boundaries, respectively. For a diagonal norm based SBPpq operator, $D \mathit{\boldsymbol{f}} = \mathit{\boldsymbol{f}}_x + \mathcal{O}(h^l)$, where $l = p$ in the interior and $l = p/2$ at the boundaries. For stable approximations, this leads, in general, to $p +1$ order accurate solutions ^[16,17].

The manufactured solution we have used is

Inserting (6.2) into (2.1) leads to a non-zero right-hand side $\vec{\mathit {k}}(t, x, y)$, which is evaluated on the grid and added to the right-hand side of (3.2) by the vector $\vec{\mathit { \mathit{\boldsymbol{k}}}}(t)$. Since $\vec{\mathit { \mathit{\boldsymbol{k}}}}$ is independent of $\vec{\mathit {\mathit{\boldsymbol{w}}}}$, it does not affect the Jacobian. The initial and boundary data are also taken from (6.2). The step size is chosen to be $\Delta t = 10^{-5}$ and the computations are terminated at $t = 1$. Next, we compute the pointwise error vector $\vec{\mathit { \mathit{\boldsymbol{e}}}}$ and its $L_2$-norm $\|\vec{\mathit { \mathit{\boldsymbol{e}}}}\|_{I_3\otimes \mathit{\boldsymbol{P}}}$. The spatial convergence rate for the SBP operators is given by $r = \log(\| \mathit{\boldsymbol{e}}\|_{i}/\| \mathit{\boldsymbol{e}}\|_{j})/\log((j-1)/(i-1))$, where $i$ and $j$ refer to the number of grid points in both spatial dimensions. The order of accuracy in space are presented in Table 2 and agree well with the theory.

Next, we consider the steady-state problem of (2.1) and (3.2), which means that the goal is to find $\vec{\mathit {\mathit{\boldsymbol{w}}}}^*$ such that

As before, we want to find an approximation to the vector $\vec{\mathit {\mathit{\boldsymbol{w}}}}^*$ which satisfies

The Jacobian of $\mathit{\boldsymbol{\mathcal{F}}}$ is $J_{ \mathit{\boldsymbol{\mathcal{F}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}}) = J_{ \mathit{\boldsymbol{\mathcal{L}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}}) - J_{ \mathit{\boldsymbol{\mathcal{S}}}}(\vec{\mathit {\mathit{\boldsymbol{w}}}})$. When the iterate $\vec{\mathit {\mathit{\boldsymbol{w}}}}^k$ is far away from $\vec{\mathit {\mathit{\boldsymbol{w}}}}^*$, Newton's method may not converge and other techniques must initially be applied. We choose the SOR method ^[1] until $\|F(\vec{\mathit {\mathit{\boldsymbol{w}}}}^k)\|_\infty$ is sufficiently small. For SOR, the next iterate is given by $\vec{\mathit {\mathit{\boldsymbol{w}}}}^{k+1} = \vec{\mathit {\mathit{\boldsymbol{w}}}}^k(1-\alpha) + (\vec{\mathit {\mathit{\boldsymbol{w}}}}^k - \mathit{\boldsymbol{h}}^k)\alpha$, where $\mathit{\boldsymbol{h}}^k$ is the Newton step from Algorithm 1 and $\alpha \in (0, 1]$.

To verify our procedure, we choose the steady manufactured solution to be ^[18]

and the computational domain is changed to $\Omega = [-0.5, 1]\times [-1, 1]$ for $\epsilon = 1/20$. Inserting (6.5) into the time-independent version of (2.1) leads to $\vec{\mathit {k}}(t, x, y) = 0$. The initial guess is $\vec{\mathit {\mathit{\boldsymbol{w}}}}^0 = (1, 1, \dots, 1)^\top$ and the tolerance $tol$ in Algorithm 1 is again set to $10^{-12}$. Table 3 shows the error and convergence rates, which again agree well with the theory. In Figure 2, the streamlines and the velocity field are illustrated for the converged solution on the grid containing $100\times 100$ points. They agree well with previous results ^[18].

6.2.   The convergence rate of the Newton iteration

Next, we will test the main development in this paper. For $\vec{\mathit {\mathit{\boldsymbol{w}}}}^k$ sufficiently close to $\vec{\mathit {\mathit{\boldsymbol{w}}}}^*$, Newton's method converges quadratically in any norm ^[1], which means that $e_{k+1} = C e_k^2$, where $C$ varies marginally between iterations and $e_k = \| \vec{\mathit {\mathit{\boldsymbol{w}}}}^k - \vec{\mathit {\mathit{\boldsymbol{w}}}}^*\|$. To verify that, we consider a grid of size $100\times 100$ with the SBP42 operator. The exact solution $\vec{\mathit {\mathit{\boldsymbol{w}}}}^*$ is approximated by the last iterate. By the assumption that $C$ is constant, the relation

is obtained for a general convergence rate $p$, which yields

The error $e_k = \| \vec{\mathit {\mathit{\boldsymbol{w}}}}^k- \vec{\mathit {\mathit{\boldsymbol{w}}}}^*\|_\infty$ is presented in Table 4 together with the estimations of $p$. The convergence rate agrees well with the expected theoretical one, which verifies that the Jacobian of $\mathit{\boldsymbol{\mathcal{F}}}$ is correct.

Next, we move on to a more realistic case where the boundary data is set to $g_1 = 1$, $g_2 = g_3 = g_4 = g_5 = g_6 = 0,$ and $\epsilon = 0.01$, which will lead to a boundary layer. The computations are performed on $\Omega = [0, 1]^2$ with $200\times 200$ grid points with the SBP42 operator. Figure 3 illustrates $\mathit{\boldsymbol{u}}$ for the converged solution and the iterative convergence order, $p$, is presented in Table 5. The estimated iterative convergence order agrees well with what is theoretically expected.

In the last experiment, we consider the incompressible Euler equation using a curvilinear grid. Both the south and north sides are solid surfaces, where the normal velocity is zero. The west side is an inflow boundary where $u = 1$ and $v = 0$ are specified, and at the east side, $p = 0$ is imposed. We change the domain to $\Omega = [-1.5, 1.5]\times [0, 0.8]$ and include a smooth bump at the south boundary given by $y(x) = 0.0625e^{-25x^2}$ ^[19]. In Figure 4, the converged solution is illustrated and the estimated iterative convergence rate $p$ is presented in Table 6 for the initial guess $(\mathit{\boldsymbol{u}}^0; \mathit{\boldsymbol{v}}^0; \mathit{\boldsymbol{p}}^0) = (1, \dots, 1; 0, \dots, 0; 1, \dots 1)$. Again, the results agree well with the theoretical value.

Remark 3. In the example above containing the curvilinear grid, the geometry is included in the definition of the derivative operators ^[20]. This simplifies the derivation of the corresponding Jacobian since no special consideration of the metric terms are required. The work in ^[20] was extended in ^[21] to included even more complex geometries with several sub-domains and non-conforming interfaces. In general, as long as the discretization can be written in a matrix-vector form (including unstructured mesh discretizations ^[22]), the additional effort to derive the Jacobian is small.

7.   Conclusions

We derived an explicit expression for the Jacobian of a finite-difference discretization of the incompressible Navier-Stokes equations. Both the Jacobian of the system of equations and the Jacobian of the related boundary condition was computed exactly. By using the block structure of the discretization, we showed that the Jacobian had a block structure as well, which leads to a compact and clean expression. We also showed that large parts of the Jacobian were computed by evaluating the discretization. We showed that the Jacobian could be used both in steady-state and time-dependent simulations. The numerical discretization was verified by manufactured solutions and the spatial convergence rates agreed well with the theoretical expectations. Furthermore, the computed estimates of the iterative convergence rates for Newton's method was two, which verified that the Jacobian was correctly computed. The methodology used in this paper is general and can be used in a straightforward manner for any numerical discretization of initial boundary value problems that can be written in matrix-vector form. The method is particularly straightforward for approximations in SBP-SAT form.

Author contributions

J. N.: Conceptualization, Methodology, Funding acquisition, Project administration, Supervision, Writing–review and editing; F. L., O. Å.: Conceptualization, Formal analysis, Methodology, Software, Validation, Visualization, Writing–original draft.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

J. N. was supported by Vetenskapsrådet, Sweden [award 2021-05484 VR] and the University of Johannesburg.

Conflict of interest

The authors declare no conflict of interest.