Speaking territories in the cultural landscape: Challenges of transformation in the Central Europe

Jana Pecnikova; Jana Pecnikova

doi:10.3934/urs.2023016

Urban Resilience and Sustainability

2023, Volume 1, Issue 4: 251-259. doi: 10.3934/urs.2023016

Previous Article Next Article

Research article Special Issues

Speaking territories in the cultural landscape: Challenges of transformation in the Central Europe

Jana Pecnikova ^,

Department of Romance Studies, Matej Bel University, Tajovskeho 40, Banska Bystrica, 97401 Slovakia

Received: 26 July 2023 Revised: 11 October 2023 Accepted: 16 October 2023 Published: 27 October 2023

Cultural landscape is a concept that characterizes the diversity of forms, urban structure, historical value and symbolic function, reflecting the fact that a cultural landscape is a materialized system that carries a variety of verbal, visual and physical aspects of human existence. This scientific paper is based on a definition of cultural landscape that consists of cultural layers which are often associated with the memory of place, and our research is based on Taylor's definition and exploration of cultural landscapes. In terms of the linear dimension of human history, cultures interact, sometimes destroy and forget, but at the same time, the cultural layers are stacked, layered and acquired in various forms and features. The methodology is based on a qualitative case study with focus on challenges connected with transformation and recovery of area that carries traumatic legacy of the past. In the conclusion, the state of historical revitalization is presented.

Keywords:

Citation: Jana Pecnikova. Speaking territories in the cultural landscape: Challenges of transformation in the Central Europe[J]. Urban Resilience and Sustainability, 2023, 1(4): 251-259. doi: 10.3934/urs.2023016

Related Papers:

[1]	Sani Aji, Poom Kumam, Aliyu Muhammed Awwal, Mahmoud Muhammad Yahaya, Kanokwan Sitthithakerngkiet . An efficient DY-type spectral conjugate gradient method for system of nonlinear monotone equations with application in signal recovery. AIMS Mathematics, 2021, 6(8): 8078-8106. doi: 10.3934/math.2021469
[2]	Khudija Bibi . Solutions of initial and boundary value problems using invariant curves. AIMS Mathematics, 2024, 9(8): 22057-22066. doi: 10.3934/math.20241072
[3]	Yasir Nadeem Anjam . The qualitative analysis of solution of the Stokes and Navier-Stokes system in non-smooth domains with weighted Sobolev spaces. AIMS Mathematics, 2021, 6(6): 5647-5674. doi: 10.3934/math.2021334
[4]	Song Lunji . A High-Order Symmetric Interior Penalty Discontinuous Galerkin Scheme to Simulate Vortex Dominated Incompressible Fluid Flow. AIMS Mathematics, 2016, 1(1): 43-63. doi: 10.3934/Math.2016.1.43
[5]	Hu Li . The modified quadrature method for Laplace equation with nonlinear boundary conditions. AIMS Mathematics, 2020, 5(6): 6211-6220. doi: 10.3934/math.2020399
[6]	Abdissalam Sarsenbi, Abdizhahan Sarsenbi . Boundary value problems for a second-order differential equation with involution in the second derivative and their solvability. AIMS Mathematics, 2023, 8(11): 26275-26289. doi: 10.3934/math.20231340
[7]	Yang Cao, Quan Shi, Sen-Lai Zhu . A relaxed generalized Newton iteration method for generalized absolute value equations. AIMS Mathematics, 2021, 6(2): 1258-1275. doi: 10.3934/math.2021078
[8]	Yasir Nadeem Anjam . Singularities and regularity of stationary Stokes and Navier-Stokes equations on polygonal domains and their treatments. AIMS Mathematics, 2020, 5(1): 440-466. doi: 10.3934/math.2020030
[9]	Yaojun Ye, Lanlan Li . Global existence and blow-up of solutions for logarithmic Klein-Gordon equation. AIMS Mathematics, 2021, 6(7): 6898-6914. doi: 10.3934/math.2021404
[10]	Fugeng Zeng, Peng Shi, Min Jiang . Global existence and finite time blow-up for a class of fractional $ p $-Laplacian Kirchhoff type equations with logarithmic nonlinearity. AIMS Mathematics, 2021, 6(3): 2559-2578. doi: 10.3934/math.2021155

Abstract

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{F}}}( \mathit{\boldsymbol{\phi}}) = 0 \end{equation}$

(1.1)

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

$\begin{array}{*{20}{c}} {\begin{array}{*{20}{r}} {\tilde{I}} \vec{\mathit {w}}_t + \mathcal{L}( \vec{\mathit {w}}) = 0 \\ \mathcal{H} \vec{\mathit {w}} = \vec{\mathit {g}}\\ {\tilde{I}} \vec{\mathit {w}} = {\tilde{I}} \vec{\mathit {f}} \end{array}}&{\begin{array}{*{20}{l}} (x, y) \in \Omega\\ (x, y) \in \partial\Omega\\ (x, y) \in \Omega \end{array}}&{\begin{array}{*{20}{l}} t > 0\\ t > 0\\ t = 0 \, . \end{array}} \end{array}$

(2.1)

Figure 1. Illustration of the computational domain

$\Omega$ and the specific boundary conditions.

DownLoad: Full-Size Img PowerPoint

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]

$\begin{equation} \mathcal{L}( \vec{\mathit {w}}) = \frac{1}{2}\left[A \vec{\mathit {w}}_x + (A \vec{\mathit {w}})_x + B \vec{\mathit {w}}_y + (B \vec{\mathit {w}})_y\right] - \epsilon\tilde{I}[ \vec{\mathit {w}}_{xx} + \vec{\mathit {w}}_{yy}] \, . \end{equation}$

(2.2)

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

$A = \begin{pmatrix} u & 0 & 1 \\ 0 & u & 0 \\ 1 & 0 & 0 \end{pmatrix} , \quad B = \begin{pmatrix} v & 0 & 0 \\ 0 & v & 1 \\ 0 & 1 & 0 \end{pmatrix} , \quad \tilde{I} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}$

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

$\begin{equation} \begin{aligned} u & = g_1 \quad & v & = g_2 & \quad & \text{at} \quad x = 0 \quad && \text{(West)} \\ p - \epsilon u_x & = g_3 \quad & -\epsilon v_x & = g_4 & \quad & \text{at} \quad x = 1 \quad && \text{(East)} \\ u & = 0 \quad & v & = 0 & \quad & \text{at} \quad y = 0 \quad && \text{(South)} \\ p-\epsilon v_y & = g_5 \quad & -\epsilon u_y & = g_6 & \quad & \text{at} \quad y = 1 \quad && \text{(North)} \, . \end{aligned} \end{equation}$

(2.3)

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

$\langle \vec{\mathit {\phi}}, \vec{\mathit {\psi}} \rangle = \int_\Omega \vec{\mathit {\phi}}^\top \vec{\mathit {\psi}} d\Omega , \quad \| \vec{\mathit {\phi}}\|^2 = \langle \vec{\mathit {\phi}}, \vec{\mathit {\phi}}\rangle \, .$

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

$\begin{equation} \frac{d}{dt}\| \vec{\mathit {w}}\|^2_{ \tilde{I}} + 2\epsilon \|\nabla \vec{\mathit {w}}\|^2_{ \tilde{I}} = BT \, , \end{equation}$

(2.4)

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

$BT = \int_{\text{South}} \vec{\mathit {w}}^\top (B \vec{\mathit {w}} - 2\epsilon \tilde{I} \vec{\mathit {w}}_y) d x$

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

$\begin{equation} \| \vec{\mathit {w}}\|^2_{\tilde{I}}(T) + 2\epsilon \int_0^T\|\nabla \vec{\mathit {w}}\|^2_{\tilde{I}} dt \le \|f\|_{\tilde{I}}^2 \, , \end{equation}$

(2.5)

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

$\mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} \approx \mathit{\boldsymbol{u}}_x, \quad \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{u}} \approx \mathit{\boldsymbol{u}}_y \, .$

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

$\begin{equation} Q_{x} + Q_{x}^\top = E_N - E_{0x} \quad Q_{y} + Q_{y}^\top = E_M - E_{0y} \, , \end{equation}$

(3.1)

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]

$\begin{equation} \mathit{\boldsymbol{\tilde{I}}} \vec{\mathit {\mathit{\boldsymbol{w}}}}_t + \mathit{\boldsymbol{\mathcal{L}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \mathit{\boldsymbol{\mathcal{S}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) \, . \end{equation}$

(3.2)

The discrete spatial operator is given by

$\begin{equation*} \begin{aligned} \mathit{\boldsymbol{\mathcal{L}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = \frac{1}{2}\left[ \mathit{\boldsymbol{A}} (I_3\otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}} + (I_3\otimes \mathit{\boldsymbol{D_x}}) \mathit{\boldsymbol{A}} \vec{\mathit {\mathit{\boldsymbol{w}}}} + \mathit{\boldsymbol{B}} (I_3\otimes \mathit{\boldsymbol{D_y}}) \vec{\mathit {\mathit{\boldsymbol{w}}}} + (I_3\otimes \mathit{\boldsymbol{D_y}}) \mathit{\boldsymbol{B}} \vec{\mathit {\mathit{\boldsymbol{w}}}}\right] \\ & - \epsilon \mathit{\boldsymbol{\tilde{I}}} [(I_3\otimes \mathit{\boldsymbol{D_x}})^2 + (I_3\otimes \mathit{\boldsymbol{D_y}})^2] \vec{\mathit {\mathit{\boldsymbol{w}}}} \, , \end{aligned} \end{equation*}$

and the block matrices are

$\mathit{\boldsymbol{A}} = \begin{pmatrix} \mathit{\boldsymbol{U}} & {\bf{0}} & \mathit{\boldsymbol{I}} \\ {\bf{0}} & \mathit{\boldsymbol{U}} & {\bf{0}} \\ \mathit{\boldsymbol{I}} & {\bf{0}} & {\bf{0}} \end{pmatrix} , \quad \mathit{\boldsymbol{B}} = \begin{pmatrix} \mathit{\boldsymbol{V}} & {\bf{0}} & {\bf{0}} \\ {\bf{0}} & \mathit{\boldsymbol{V}} & \mathit{\boldsymbol{I}} \\ {\bf{0}} & \mathit{\boldsymbol{I}} & {\bf{0}} \end{pmatrix} , \quad \mathit{\boldsymbol{\tilde{I}}} = \begin{pmatrix} \mathit{\boldsymbol{I}} & {\bf{0}} & {\bf{0}} \\ {\bf{0}} & \mathit{\boldsymbol{I}} & {\bf{0}} \\ {\bf{0}} & {\bf{0}} & {\bf{0}} \end{pmatrix} \, ,$

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{S}}}^k( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = (I_3\otimes \mathit{\boldsymbol{P}}^{-1})\Sigma^k(I_3\otimes \mathit{\boldsymbol{P}}^{k})( \mathit{\boldsymbol{\mathcal{H}}}^k \vec{\mathit {\mathit{\boldsymbol{w}}}} - \vec{ \mathit{\boldsymbol{g}}}) \, . \end{equation}$

(3.3)

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

$\begin{equation*} \mathit{\boldsymbol{P}}^k = \begin{cases} E_{0x}\otimes P_y \quad & \text{on the west boundary } (k = W) \\ E_N\otimes P_y \quad & \text{on the east boundary } (k = E) \\ P_x\otimes E_{0y} \quad & \text{on the south boundary } (k = S) \\ P_x\otimes E_M \quad & \text{on the north boundary } (k = N) \, . \end{cases} \end{equation*}$

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

$\begin{equation} \begin{aligned} \mathit{\boldsymbol{\mathcal{S}}}^W( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1}) \Sigma^W (I_3\otimes \mathit{\boldsymbol{P}}^W) \underbrace{ \begin{pmatrix} \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1 \\ \mathit{\boldsymbol{v}} - \mathit{\boldsymbol{g}}_2 \\ \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1 \end{pmatrix} }_{ \mathit{\boldsymbol{\mathcal{H}}}^W \vec{\mathit {\mathit{\boldsymbol{w}}}} - \vec{\mathit { \mathit{\boldsymbol{g}}}}} \\ \mathit{\boldsymbol{\mathcal{S}}}^E( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1}) \Sigma^E (I_3\otimes \mathit{\boldsymbol{P}}^E) \underbrace{ \begin{pmatrix} \mathit{\boldsymbol{p}} - \epsilon \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_3 \\ -\epsilon \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{v}} - \mathit{\boldsymbol{g}}_4 \\ {\bf{0}} \end{pmatrix} }_{ \mathit{\boldsymbol{\mathcal{H}}}^E \vec{\mathit {\mathit{\boldsymbol{w}}}} - \vec{\mathit { \mathit{\boldsymbol{g}}}}} \\ \mathit{\boldsymbol{\mathcal{S}}}^S( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1}) \Sigma^S(I_3\otimes \mathit{\boldsymbol{P}}^S) \underbrace{ \begin{pmatrix} \mathit{\boldsymbol{u}} - {\bf{0}} \\ \mathit{\boldsymbol{v}} - {\bf{0}} \\ \mathit{\boldsymbol{v}} - {\bf{0}} \end{pmatrix} }_{ \mathit{\boldsymbol{\mathcal{H}}}^S \vec{\mathit {\mathit{\boldsymbol{w}}}} - {\bf{0}}} \\ \mathit{\boldsymbol{\mathcal{S}}}^N( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1}) \Sigma^N (I_3\otimes \mathit{\boldsymbol{P}}^N) \underbrace{ \begin{pmatrix} - \epsilon \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_6 \\ \mathit{\boldsymbol{p}} -\epsilon \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{v}} - \mathit{\boldsymbol{g}}_5 \\ {\bf{0}} \end{pmatrix} }_{ \mathit{\boldsymbol{\mathcal{H}}}^N \vec{\mathit {\mathit{\boldsymbol{w}}}} - \vec{\mathit { \mathit{\boldsymbol{g}}}}} \, , \end{aligned} \end{equation}$

(3.4)

where

$\begin{align*} \Sigma^W & = \begin{pmatrix} - \mathit{\boldsymbol{U}}/2 + \epsilon \mathit{\boldsymbol{D_x}}^\top & {\bf{0}} & {\bf{0}} \\ {\bf{0}} & - \mathit{\boldsymbol{U}}/2 + \epsilon \mathit{\boldsymbol{D_x}}^\top & {\bf{0}} \\ {\bf{0}} & {\bf{0}} & - \mathit{\boldsymbol{I}} \end{pmatrix} , \quad &\Sigma^E & = (I_3\otimes \mathit{\boldsymbol{I}}) \\ \Sigma^S & = \begin{pmatrix} - \mathit{\boldsymbol{V}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top & {\bf{0}} & {\bf{0}} \\ {\bf{0}} & - \mathit{\boldsymbol{V}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top & {\bf{0}} \\ {\bf{0}} & {\bf{0}} & - \mathit{\boldsymbol{I}} \end{pmatrix} , \quad &\Sigma^N & = (I_3\otimes \mathit{\boldsymbol{I}}) \, . \end{align*}$

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

$\begin{equation} \begin{aligned} \mathit{\boldsymbol{\mathcal{S}}}^S( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3 \otimes \mathit{\boldsymbol{P}}^{-1}) \begin{pmatrix} - \mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}} \\ - \mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}} \\ - \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}} \end{pmatrix} \, , \end{aligned} \end{equation}$

(3.5)

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{S}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \sum\limits_{k\in\{W, E, S, N\}} \mathit{\boldsymbol{\mathcal{S}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}})^k \, . \end{equation}$

(3.6)

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

$\begin{equation} \frac{d}{dt}\| \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}} + 2\epsilon \|\nabla \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}} = \mathit{\boldsymbol{BT}} \, , \end{equation}$

(3.7)

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

$\begin{equation} \begin{aligned} \mathit{\boldsymbol{BT}} = & \underbrace{ \vec{\mathit {\mathit{\boldsymbol{w}}}}^\top(I_3 \otimes \mathit{\boldsymbol{P}}^S) \mathit{\boldsymbol{B}} \vec{\mathit {\mathit{\boldsymbol{w}}}} - 2\epsilon \vec{\mathit {\mathit{\boldsymbol{w}}}}^\top (I_3\otimes \mathit{\boldsymbol{P}}^S) \mathit{\boldsymbol{\tilde{I}}} (I_3 \otimes \mathit{\boldsymbol{D_y}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}}_I \\ & \underbrace{ + 2 \vec{\mathit {\mathit{\boldsymbol{w}}}} (I_3\otimes \mathit{\boldsymbol{P}}) \mathit{\boldsymbol{\mathcal{S}}}^S( \vec{\mathit {\mathit{\boldsymbol{w}}}})}_{II} \end{aligned} \end{equation}$

(3.8)

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

$\begin{equation*} \begin{aligned} \mathit{\boldsymbol{BT}} & = && \underbrace{ \mathit{\boldsymbol{v}}^\top \mathit{\boldsymbol{P}}^S( \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{V}} \mathit{\boldsymbol{v}} + 2 \mathit{\boldsymbol{p}} - 2 \epsilon \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{v}} ) -2 \epsilon \mathit{\boldsymbol{u}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{u}}}_I \\ &&& \underbrace{ - 2 \mathit{\boldsymbol{v}}^\top \mathit{\boldsymbol{P}}^S( \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}}/2 + \mathit{\boldsymbol{V}} \mathit{\boldsymbol{v}}/2 + \mathit{\boldsymbol{p}}^\top \mathit{\boldsymbol{v}} - \epsilon \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{v}}) + 2\epsilon \mathit{\boldsymbol{u}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{u}}}_{II} = 0 \, , \end{aligned} \end{equation*}$

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

$\begin{equation*} \| \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}} (T) + 2\epsilon \int_0^T\|\nabla \vec{\mathit {\mathit{\boldsymbol{w}}}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}} dt \le \| \mathit{\boldsymbol{f}}\|^2_{ \tilde{I}\otimes \mathit{\boldsymbol{P}}} \, , \end{equation*}$

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

$J_{ \mathit{\boldsymbol{h}}} = \begin{pmatrix} \frac{\partial h_{00}}{\partial u_{00}} & \dots & \frac{\partial h_{00}}{\partial u_{NM}} \\ \vdots &\ddots & \vdots \\ \frac{\partial h_{NM}}{\partial u_{00}} & \dots & \frac{\partial h_{NM}}{\partial u_{NM}} \end{pmatrix} \, .$

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

$\mathit{\boldsymbol{h}}( \mathit{\boldsymbol{u}}) = \begin{pmatrix} h_{00} \\ \vdots \\ h_{NM} \end{pmatrix} = \begin{pmatrix} u_{00} \\ \vdots \\ u_{NM} \end{pmatrix} = \mathit{\boldsymbol{u}} \, .$

Since

$\begin{align*} &\frac{\partial h_{00}}{\partial u_{00}} = 1 &\frac{\partial h_{00}}{\partial u_{01}} = 0 &&\frac{\partial h_{00}}{\partial u_{02}} = 0 &&\dots &&\frac{\partial h_{00}}{\partial u_{NM}} = 0 \\ &\frac{\partial h_{01}}{\partial u_{00}} = 0 &\frac{\partial h_{01}}{\partial u_{01}} = 1 &&\frac{\partial h_{01}}{\partial u_{02}} = 0 &&\dots &&\frac{\partial h_{01}}{\partial u_{NM}} = 0 \\ & & & \vdots \\ &\frac{\partial h_{NM}}{\partial u_{00}} = 0 &\frac{\partial h_{NM}}{\partial u_{01}} = 0 &&\frac{\partial h_{NM}}{\partial u_{02}} = 0 &&\dots &&\frac{\partial h_{NM}}{\partial u_{NM}} = 1 \, , \end{align*}$

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

$\begin{equation*} \begin{aligned} \mathit{\boldsymbol{h}}( \mathit{\boldsymbol{u}}) & = \begin{pmatrix} h_{00} \\ \vdots \\ h_{NM} \end{pmatrix} = \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} = \begin{pmatrix} D_{0, 0} & \dots & D_{0, NM} \\ \vdots &\ddots & \vdots \\ D_{NM, 0} & \dots & D_{NM, MN} \end{pmatrix} \begin{pmatrix} u_{00} \\ \vdots \\u_{NM} \end{pmatrix} \\ & = \begin{pmatrix} D_{0, 0}u_{00} + \dots + D_{0, NM}u_{NM} \\ \ \vdots \quad \quad \\ D_{NM, 0} u_{00} + \dots + D_{NM, MN}u_{NM} \end{pmatrix} \, . \end{aligned} \end{equation*}$

Then, in the same way

$\begin{align*} &\frac{\partial h_{00}}{\partial u_{00}} = D_{0, 0} &\frac{\partial h_{00}}{\partial u_{01}} = D_{0, 1} &&\frac{\partial h_{00}}{\partial u_{02}} = D_{0, 2} &&\dots &&\frac{\partial h_{00}}{\partial u_{NM}} = D_{0, NM} \\ &\frac{\partial h_{01}}{\partial u_{00}} = D_{1, 0} &\frac{\partial h_{01}}{\partial u_{01}} = D_{1, 1} &&\frac{\partial h_{01}}{\partial u_{02}} = D_{1, 2} &&\dots &&\frac{\partial h_{01}}{\partial u_{NM}} = D_{1, NM} \\ & & & \vdots \\ &\frac{\partial h_{NM}}{\partial u_{00}} = D_{NM, 0} &\frac{\partial h_{NM}}{\partial u_{01}} = D_{NM, 1} &&\frac{\partial h_{NM}}{\partial u_{02}} = D_{NM, 2} &&\dots &&\frac{\partial h_{NM}}{\partial u_{NM}} = D_{NM, NM} \, . \end{align*}$

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

$\mathit{\boldsymbol{h}}( \mathit{\boldsymbol{u}}) = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} = \begin{pmatrix} u_{00}[D_{0, 0}u_{00} + \dots + D_{0, NM}u_{NM}] \\ \vdots \\ u_{NM}[D_{NM, 0} u_{00} + \dots + D_{NM, MN}u_{NM}] \end{pmatrix} = \begin{pmatrix} u_{00}( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{00} \\ \vdots \\ u_{NM}( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{NM} \end{pmatrix} \, .$

By using the product rule, we get that

$\begin{align*} &\frac{\partial h_{00}}{\partial u_{00}} = u_{00} D_{0, 0} + ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{00} &&\frac{\partial h_{00}}{\partial u_{01}} = u_{00} D_{0, 1} & \dots\quad &\frac{\partial h_{00}}{\partial u_{NM}} = u_{00}D_{0, NM} \\ &\frac{\partial h_{01}}{\partial u_{00}} = u_{01} D_{1, 0} &&\frac{\partial h_{01}}{\partial u_{01}} = u_{01} D_{1, 1} + ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{01} & \dots \quad &\frac{\partial h_{01}}{\partial u_{NM}} = u_{00}D_{0, NM} \\ & \quad \quad \vdots \\ &\frac{\partial h_{NM}}{\partial u_{00}} = u_{NM} D_{NM, 0} &&\frac{\partial h_{NM}}{\partial u_{01}} = u_{NM} D_{NM, 1} & \dots\quad &\frac{\partial h_{NM}}{\partial u_{NM}} = u_{NM} D_{NM, NM} + ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{NM} \, . \end{align*}$

Hence,

$\begin{equation*} \begin{aligned} J_{ \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}} & = \begin{pmatrix} u_{00} D_{0, 0} + ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{00}& u_{00} D_{0, 1} & \dots & u_{00}D_{0, NM} \\ u_{01} D_{1, 0} & u_{01} D_{1, 1} + ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{01} & \dots & u_{01}D_{1, NM} \\ && \vdots \\ u_{NM} D_{NM, 0} & u_{NM} D_{NM, 1} & \dots & u_{NM}D_{NM, NM} + ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{NM} \end{pmatrix} \\ & = \underbrace{ \begin{pmatrix} u_{00} D_{0, 0} & u_{00} D_{0, 1} & \dots & u_{00}D_{0, NM} \\ u_{01} D_{1, 0} & u_{01} D_{1, 1} & \dots & u_{01}D_{1, NM} \\ && \vdots \\ u_{NM} D_{NM, 0} & u_{NM} D_{NM, 1} & \dots & u_{NM}D_{NM, NM} \end{pmatrix} }_{ \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}}} \\ & + \underbrace{ \begin{pmatrix} ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{00} & 0 & \dots & 0 \\ 0 & ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{01} & \dots & 0 \\ && \vdots \\ 0 & 0 & \dots & ( \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}})_{NM} \end{pmatrix} }_{\underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}}} = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} + \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}} \, , \end{aligned} \end{equation*}$

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

In a similar manner, for

$\mathit{\boldsymbol{h}}( \mathit{\boldsymbol{u}}) = \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}} = \begin{pmatrix} D_{0, 0}u_{00}^2 + \dots + D_{0, NM}u_{NM}^2 \\ \vdots \\ D_{NM, 0} u_{00}^2 + \dots + D_{NM, MN}u_{NM}^2 \end{pmatrix}$

we get that

$\begin{align*} &\frac{\partial h_{00}}{\partial u_{00}} = 2 D_{0, 0}u_{00} &\frac{\partial h_{00}}{\partial u_{01}} = 2 D_{0, 1}u_{01} \quad \dots \quad &\frac{\partial h_{00}}{\partial u_{NM}} = 2 D_{0, NM}u_{NM} \\ &\frac{\partial h_{01}}{\partial u_{00}} = 2 D_{1, 0}u_{00} &\frac{\partial h_{01}}{\partial u_{01}} = 2 D_{1, 1}u_{01} \quad \dots \quad &\frac{\partial h_{01}}{\partial u_{NM}} = 2 D_{1, NM}u_{NM} \\ & \quad \quad \vdots \\ &\frac{\partial h_{NM}}{\partial u_{00}} = 2 D_{NM, 0}u_{00} &\frac{\partial h_{NM}}{\partial u_{01}} = 2 D_{NM, 1}u_{01} \quad \dots \quad &\frac{\partial h_{01}}{\partial u_{NM}} = 2 D_{NM, NM}u_{NM} \, . \end{align*}$

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

$\begin{equation} J_{ \mathit{\boldsymbol{u}}} = \mathit{\boldsymbol{I}}, \quad J_{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} } = \mathit{\boldsymbol{D_x}}, \quad J_{ \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}} = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} + \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}}, \quad J_{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}}} = 2 \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \, . \end{equation}$

(4.1)

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

$\mathit{\boldsymbol{h}}^1 = \begin{pmatrix} h^1_{00} \\ h^1_{01} \\ \vdots \\ h^1_{NM} \end{pmatrix} \quad \text{and} \quad \vec{\mathit {\mathit{\boldsymbol{w}}}} = \begin{pmatrix} \mathit{\boldsymbol{u}} \\ \mathit{\boldsymbol{v}} \\ \mathit{\boldsymbol{p}} \end{pmatrix}$

it follows that

$\begin{align*} J_{h^1_{00}} = \frac{\partial h^1_{00}}{\partial \vec{\mathit {\mathit{\boldsymbol{w}}}}} & = \begin{pmatrix} \frac{\partial h^1_{00}}{\partial w_{00}} & \dots & \frac{\partial h^1_{00}}{\partial w_{3NM}} \end{pmatrix} = \begin{pmatrix} \frac{\partial h^1_{00}}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial h^1_{00}}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial h^1_{00}}{\partial \mathit{\boldsymbol{p}}} \end{pmatrix} \in \mathbb{R}^{1\times 3n} \end{align*}$

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

$J_{ \mathit{\boldsymbol{h}}^1} = \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \vec{\mathit {\mathit{\boldsymbol{w}}}}} = \begin{pmatrix} \frac{\partial h^1_{00}}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial h^1_{00}}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial h^1_{00}}{\partial \mathit{\boldsymbol{p}}} \\ \frac{\partial h^1_{01}}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial h^1_{01}}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial h^1_{01}}{\partial \mathit{\boldsymbol{p}}} \\ \vdots & \vdots & \vdots \\ \frac{\partial h^1_{NM}}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial h^1_{NM}}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial h^1_{NM}}{\partial \mathit{\boldsymbol{p}}} \end{pmatrix} = \begin{pmatrix} \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \mathit{\boldsymbol{p}}} \end{pmatrix} \in \mathbb{R}^{n\times 3n} \, .$

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

$J_{\tilde{ \mathit{\boldsymbol{h}}}} = \begin{pmatrix} \frac{\partial \mathit{\boldsymbol{h}}^1 }{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{h}}^1 }{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{h}}^1 }{\partial \mathit{\boldsymbol{p}}} \\ \frac{\partial \mathit{\boldsymbol{h}}^2 }{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{h}}^2 }{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{h}}^2 }{\partial \mathit{\boldsymbol{p}}} \\ \frac{\partial \mathit{\boldsymbol{h}}^3 }{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{h}}^3 }{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{h}}^3 }{\partial \mathit{\boldsymbol{p}}} \end{pmatrix} \in \mathbb{R}^{3n\times 3n}$

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

$\tilde{ \mathit{\boldsymbol{h}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \begin{pmatrix} \mathit{\boldsymbol{h}}^1( \vec{\mathit {\mathit{\boldsymbol{w}}}}) \\ \mathit{\boldsymbol{h}}^2( \vec{\mathit {\mathit{\boldsymbol{w}}}}) \\ \mathit{\boldsymbol{h}}^3( \vec{\mathit {\mathit{\boldsymbol{w}}}}) \end{pmatrix} = \mathit{\boldsymbol{A}} (I_3 \otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}} = \begin{pmatrix} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{p}} \\ \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{v}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} \end{pmatrix} = \begin{pmatrix} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{p}} \\ \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{v}}}} \mathit{\boldsymbol{u}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} \end{pmatrix} \, .$

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

$\begin{align*} \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \mathit{\boldsymbol{u}}} & = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} + \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}} \quad & \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \mathit{\boldsymbol{v}}} & = {\bf{0}} \quad & \frac{\partial \mathit{\boldsymbol{h}}^1}{\partial \mathit{\boldsymbol{p}}} & = \mathit{\boldsymbol{D_x}} \\ \frac{\partial \mathit{\boldsymbol{h}}^2}{\partial \mathit{\boldsymbol{u}}} & = \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{v}}}} \quad &\frac{\partial \mathit{\boldsymbol{h}}^2}{\partial \mathit{\boldsymbol{v}}} & = \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} \quad &\frac{\partial \mathit{\boldsymbol{h}}^2}{\partial \mathit{\boldsymbol{p}}} & = {\bf{0}} \\ \frac{\partial \mathit{\boldsymbol{h}}^3}{\partial \mathit{\boldsymbol{u}}} & = \mathit{\boldsymbol{D_x}} \quad &\frac{\partial \mathit{\boldsymbol{h}}^3}{\partial \mathit{\boldsymbol{v}}} & = {\bf{0}} \quad &\frac{\partial \mathit{\boldsymbol{h}}^3}{\partial \mathit{\boldsymbol{p}}} & = {\bf{0}} \, . \end{align*}$

Thus,

$J_{ \mathit{\boldsymbol{A}} (I_3 \otimes \mathit{\boldsymbol{D_x}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}} = \begin{pmatrix} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} + \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}} & {\bf{0}} & \mathit{\boldsymbol{D_x}} \\ \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{v}}}} & \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} & {\bf{0}} \\ \mathit{\boldsymbol{D_x}} & {\bf{0}} & {\bf{0}} \end{pmatrix} \, .$

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

$(I_3 \otimes \mathit{\boldsymbol{D_x}}) \mathit{\boldsymbol{A}} \vec{\mathit {\mathit{\boldsymbol{w}}}} = \begin{pmatrix} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{p}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{v}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} \end{pmatrix} = \begin{pmatrix} \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}} + \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{p}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{V}} \mathit{\boldsymbol{u}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}} \end{pmatrix} \, ,$

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

$J_{(I_3 \otimes \mathit{\boldsymbol{D_x}}) \mathit{\boldsymbol{A}} \vec{\mathit {\mathit{\boldsymbol{w}}}}} = \begin{pmatrix} 2 \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} & {\bf{0}} & \mathit{\boldsymbol{D_x}} \\ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{V}} & \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} & {\bf{0}} \\ \mathit{\boldsymbol{D_x}} & {\bf{0}} & {\bf{0}} \end{pmatrix} \, .$

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

$\begin{align*} J_{ \mathit{\boldsymbol{B}} (I_3 \otimes \mathit{\boldsymbol{D_y}}) \vec{\mathit {\mathit{\boldsymbol{w}}}}} = & \begin{pmatrix} \mathit{\boldsymbol{V}} \mathit{\boldsymbol{D_y}} & \underline{{ \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{u}}}} & {\bf{0}} \\ {\bf{0}} & \mathit{\boldsymbol{V}} \mathit{\boldsymbol{D_y}} + \underline{{ \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{v}}}} & \mathit{\boldsymbol{D_y}} \\ {\bf{0}} & \mathit{\boldsymbol{D_y}} & {\bf{0}} \end{pmatrix} \\ J_{(I_3 \otimes \mathit{\boldsymbol{D_y}}) \mathit{\boldsymbol{B}} \vec{\mathit {\mathit{\boldsymbol{w}}}}} = & \begin{pmatrix} \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{V}} & \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{U}} & {\bf{0}} \\ {\bf{0}} & 2 \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{V}} & \mathit{\boldsymbol{D_y}} \\ {\bf{0}} & \mathit{\boldsymbol{D_y}} & {\bf{0}} \end{pmatrix} \, . \end{align*}$

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

$J_{-\epsilon \mathit{\boldsymbol{\tilde{I}}} [(I_3\otimes \mathit{\boldsymbol{D_x}})^2 +(I_3\otimes \mathit{\boldsymbol{D_y}})^2] \vec{\mathit {\mathit{\boldsymbol{w}}}}} = - \begin{pmatrix} \epsilon ( \mathit{\boldsymbol{D_x}}^2 + \mathit{\boldsymbol{D_y}}^2) & {\bf{0}} & {\bf{0}} \\ {\bf{0}} & \epsilon ( \mathit{\boldsymbol{D_x}}^2 + \mathit{\boldsymbol{D_y}}^2) & {\bf{0}} \\ {\bf{0}} & {\bf{0}} & {\bf{0}} \end{pmatrix} \, .$

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

$\begin{equation} J_{ \mathit{\boldsymbol{\mathcal{L}}}} = \begin{pmatrix} J_{11} & \frac{1}{2} \left(\underline{{ \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{u}}}} + \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{U}}\right)& \mathit{\boldsymbol{D_x}} \\ \frac{1}{2} \left(\underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{v}}}} + \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{V}}\right) & J_{22} & \mathit{\boldsymbol{D_y}} \\ \mathit{\boldsymbol{D_x}} & \mathit{\boldsymbol{D_y}} & {\bf{0}} \\ \end{pmatrix} \end{equation}$

(4.2)

where

$\begin{align*} J_{11} & = \frac{1}{2} \left( \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} + \underline{{ \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{u}}}} + 2 \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}} + \mathit{\boldsymbol{V}} \mathit{\boldsymbol{D_y}} + \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{V}}\right) - \epsilon ( \mathit{\boldsymbol{D_x}}^2 + \mathit{\boldsymbol{D_y}}^2)\\ J_{22} & = \frac{1}{2} \left( \mathit{\boldsymbol{V}} \mathit{\boldsymbol{D_y}} + \underline{{ \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{v}}}} + 2 \mathit{\boldsymbol{D_y}} \mathit{\boldsymbol{V}} + \mathit{\boldsymbol{U}} \mathit{\boldsymbol{D_x}} + \mathit{\boldsymbol{D_x}} \mathit{\boldsymbol{U}}\right) - \epsilon ( \mathit{\boldsymbol{D_x}}^2 + \mathit{\boldsymbol{D_y}}^2) \, . \end{align*}$

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{S}}}^S( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \begin{pmatrix} \mathit{\boldsymbol{\mathcal{S}}}_1^S \\ \mathit{\boldsymbol{\mathcal{S}}}_2^S \\ \mathit{\boldsymbol{\mathcal{S}}}_3^S \\ \end{pmatrix} = (I_3\otimes \mathit{\boldsymbol{P}}^{-1}) \begin{pmatrix} - \mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}} \\ - \mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}} /2+ \epsilon \mathit{\boldsymbol{D_y}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}} \\ - \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}} \end{pmatrix} \in \mathbb{R}^{3n} \, . \end{equation}$

(4.3)

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

$J_{ \mathit{\boldsymbol{\mathcal{S}}}^S} = \begin{pmatrix} \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_1^S}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_1^S}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_1^S}{\partial \mathit{\boldsymbol{p}}} \\ \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_2^S}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_2^S}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_2^S}{\partial \mathit{\boldsymbol{p}}} \\ \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_3^S}{\partial \mathit{\boldsymbol{u}}} & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_3^S}{\partial \mathit{\boldsymbol{v}}} & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_3^S}{\partial \mathit{\boldsymbol{p}}} \end{pmatrix} \in \mathbb{R}^{3n\times 3n} \, .$

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

$\begin{align*} \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_1^S}{\partial \mathit{\boldsymbol{u}}} & = \Big(- \underbrace{\frac{\partial}{\partial \mathit{\boldsymbol{u}}} \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}}/2}_{ = \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{V}} \mathit{\boldsymbol{P}}^S/2} + \underbrace{ \frac{\partial}{\partial \mathit{\boldsymbol{u}}} \epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_y}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}}}_{ = \epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_y}}^\top \mathit{\boldsymbol{P}}^S} \Big) = \mathit{\boldsymbol{P}}^{-1} \left(- \mathit{\boldsymbol{V}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top\right) \mathit{\boldsymbol{P}}^S \, . \end{align*}$

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

$\begin{align*} \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_1^S}{\partial \mathit{\boldsymbol{v}}} & = \Big(- \underbrace{\frac{\partial}{\partial \mathit{\boldsymbol{v}}} \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{v}}/2}_{ = \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^S /2} + \underbrace{ \frac{\partial}{\partial \mathit{\boldsymbol{v}}}\epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^S \mathit{\boldsymbol{u}}}_{ = {\bf{0}}} \Big) = - \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^S /2 \, . \end{align*}$

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

$\begin{align*} \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}^W_2}{\partial \mathit{\boldsymbol{v}}} & = \mathit{\boldsymbol{P}}^{-1}(- \mathit{\boldsymbol{V}} + \epsilon \mathit{\boldsymbol{D_y}}^\top) \mathit{\boldsymbol{P}}^S , \quad & \frac{\partial \mathit{\boldsymbol{\mathcal{S}}}_3^W}{\partial \mathit{\boldsymbol{v}}} & = - \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^S \, , \end{align*}$

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,

$J_{ \mathit{\boldsymbol{\mathcal{S}}}^S} = (I_3 \otimes \mathit{\boldsymbol{P}}^{-1}) \begin{pmatrix} - \mathit{\boldsymbol{V}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top & - \mathit{\boldsymbol{U}}/2 & {\bf{0}} \\ {\bf{0}} & - \mathit{\boldsymbol{V}} + \epsilon \mathit{\boldsymbol{D_y}}^\top & {\bf{0}} \\ {\bf{0}} & - \mathit{\boldsymbol{I}} & {\bf{0}} \end{pmatrix} (I_3 \otimes \mathit{\boldsymbol{P}}^S) \, .$

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

$\begin{align*} \mathit{\boldsymbol{\mathcal{S}}}^W_1( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = \mathit{\boldsymbol{P}}^{-1}(- \mathit{\boldsymbol{U}}/2 + \epsilon \mathit{\boldsymbol{D_x}}^\top) \mathit{\boldsymbol{P}}^W( \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1) \\ & = - \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{U}} \mathit{\boldsymbol{P}}^W( \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1)/2 + \epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^W( \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1) \, . \end{align*}$

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

$\begin{align*} \frac{\partial}{\partial \mathit{\boldsymbol{u}}}\left(\epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^W( \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1) \right) = \underbrace{ \frac{\partial}{\partial \mathit{\boldsymbol{u}}}\left(\epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^W \mathit{\boldsymbol{u}}\right) }_{ = \epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^W} - \underbrace{ \frac{\partial}{\partial \mathit{\boldsymbol{u}}}\left(\epsilon \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{D_x}}^\top \mathit{\boldsymbol{P}}^W \mathit{\boldsymbol{g}}_1\right) }_{ = 0} \, . \end{align*}$

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

$\begin{align*} -\frac{\partial}{\partial \mathit{\boldsymbol{u}}}\left( \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^W \mathit{\boldsymbol{U}}( \mathit{\boldsymbol{u}} - \mathit{\boldsymbol{g}}_1)/2\right) & = \underbrace{ -\frac{\partial}{\partial \mathit{\boldsymbol{u}}}\left( \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^W \mathit{\boldsymbol{U}} \mathit{\boldsymbol{u}})/2\right)}_ { = - \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^W \mathit{\boldsymbol{U}}} + \underbrace{ \frac{\partial}{\partial \mathit{\boldsymbol{u}}}\left( \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^W\underline{{ \mathit{\boldsymbol{g}}_1}} \mathit{\boldsymbol{u}})/2\right) }_{ = \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^W \underline{{ \mathit{\boldsymbol{g}}_1}}/2} \\ & = \mathit{\boldsymbol{P}}^{-1}(- \mathit{\boldsymbol{U}} + \underline{{ \mathit{\boldsymbol{g}}_1}}/2) \mathit{\boldsymbol{P}}^W \end{align*}$

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

$J_{ \mathit{\boldsymbol{\mathcal{S}}}^W_1}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \begin{pmatrix} \mathit{\boldsymbol{P}}^{-1}(-( \mathit{\boldsymbol{U}} - \underline{{ \mathit{\boldsymbol{g}}_1}}/2) + \epsilon \mathit{\boldsymbol{D_x}}^\top) \mathit{\boldsymbol{P}}^W & {\bf{0}} & {\bf{0}} \end{pmatrix} \in \mathbb{R}^{n\times 3n}$

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

$\begin{equation} J_{ \mathit{\boldsymbol{\mathcal{S}}}} ( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \sum\limits_{k\in\{W, E, S, N\}} J_{ \mathit{\boldsymbol{\mathcal{S}}}^k}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) \, , \end{equation}$

(4.4)

where

$\begin{equation*} \begin{aligned} J_{ \mathit{\boldsymbol{\mathcal{S}}}^W}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1}) \begin{pmatrix} -( \mathit{\boldsymbol{U}}-\underline{{ \mathit{\boldsymbol{g}}_1}}/2) + \epsilon \mathit{\boldsymbol{D_x}}^\top & {\bf{0}} & {\bf{0}} \\ -( \mathit{\boldsymbol{V}}-\underline{{ \mathit{\boldsymbol{g}}_2}})/2 & - \mathit{\boldsymbol{U}}/2 + \epsilon \mathit{\boldsymbol{D_x}}^\top & {\bf{0}} \\ - \mathit{\boldsymbol{I}} & {\bf{0}} & {\bf{0}} \end{pmatrix} (I_3\otimes \mathit{\boldsymbol{P}}^W) \\ J_{ \mathit{\boldsymbol{\mathcal{S}}}^E}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^E) \begin{pmatrix} - \epsilon \mathit{\boldsymbol{D_x}} & {\bf{0}} & \mathit{\boldsymbol{I}} \\ {\bf{0}} & - \epsilon \mathit{\boldsymbol{D_x}} & {\bf{0}} \\ {\bf{0}} & {\bf{0}} & {\bf{0}} \end{pmatrix} \\ J_{ \mathit{\boldsymbol{\mathcal{S}}}^S}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^S) \begin{pmatrix} - \mathit{\boldsymbol{V}}/2 + \epsilon \mathit{\boldsymbol{D_y}}^\top& - \mathit{\boldsymbol{U}}/2 & {\bf{0}} \\ {\bf{0}} & - \mathit{\boldsymbol{V}} + \epsilon \mathit{\boldsymbol{D_y}}^\top & {\bf{0}} \\ {\bf{0}} & - \mathit{\boldsymbol{I}} & {\bf{0}} \end{pmatrix} (I_3\otimes \mathit{\boldsymbol{P}}^S) \\ J_{ \mathit{\boldsymbol{\mathcal{S}}}^N}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) & = (I_3\otimes \mathit{\boldsymbol{P}}^{-1} \mathit{\boldsymbol{P}}^N) \begin{pmatrix} - \epsilon \mathit{\boldsymbol{D_y}} & {\bf{0}} & {\bf{0}} \\ {\bf{0}} & - \epsilon \mathit{\boldsymbol{D_y}} & \mathit{\boldsymbol{I}} \\ {\bf{0}} & {\bf{0}} & {\bf{0}} \end{pmatrix} \, . \end{aligned} \end{equation*}$

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

$\begin{equation*} \mathcal{M} \mathit{\boldsymbol{\phi}}_t + \mathit{\boldsymbol{\mathcal{H}}}( \mathit{\boldsymbol{\phi}}) = 0 \, , \end{equation*}$

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

$\begin{equation} \frac{ \mathcal{M} ( \mathit{\boldsymbol{\phi}}^{i+1} - \mathit{\boldsymbol{\phi}}^i)}{\Delta t} + \mathit{\boldsymbol{\mathcal{H}}}( \mathit{\boldsymbol{\phi}}^{i+1}) = 0 \, . \end{equation}$

(5.1)

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{F}}}( \mathit{\boldsymbol{\phi}}^{i+1}) = \frac{ \mathcal{M} ( \mathit{\boldsymbol{\phi}}^{i+1} - \mathit{\boldsymbol{\phi}}^i)}{\Delta t} + \mathit{\boldsymbol{\mathcal{H}}}( \mathit{\boldsymbol{\phi}}^{i+1}) \, . \end{equation}$

(5.2)

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}$ .

Algorithm 1 Newton's method
1: Input: $\mathit{\boldsymbol{\phi}}^0$ and tolerance $tol$
2: Output: An approximation of $\mathit{\boldsymbol{\phi}}^$ , where $\mathit{\boldsymbol{\mathcal{F}}}(\mathit{\boldsymbol{\phi}}^) = 0$
3: for $j = 0, 1, 2, \dots$ do
4: solve $J_{ \mathit{\boldsymbol{\mathcal{F}}}}(\mathit{\boldsymbol{\phi}}^j)\mathit{\boldsymbol{h}}^j = - \mathit{\boldsymbol{\mathcal{F}}}(\mathit{\boldsymbol{\phi}}^j)$
5: set $\mathit{\boldsymbol{\phi}}^{j+1} = \mathit{\boldsymbol{\phi}}^j + \mathit{\boldsymbol{h}}^j$
6: if $\\| \mathit{\boldsymbol{\mathcal{F}}}(\mathit{\boldsymbol{\phi}}^{j+1})\\| < tol$ then
7: return $\mathit{\boldsymbol{\phi}}^{j+1}$
8: end if
9: end for

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{F}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^{i+1}) = \frac{1}{\Delta t} \left( \begin{bmatrix} \mathit{\boldsymbol{u}}^{i+1} \\ \mathit{\boldsymbol{v}}^{i+1} \\ {\bf{0}} \end{bmatrix} - \begin{bmatrix} \mathit{\boldsymbol{u}}^{i} \\ \mathit{\boldsymbol{v}}^{i} \\ {\bf{0}} \end{bmatrix} \right) + \mathit{\boldsymbol{\mathcal{L}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^{i+1}) - \mathit{\boldsymbol{\mathcal{S}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^{i+1}) \, . \end{equation}$

(5.3)

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

$\begin{equation} J_{ \mathit{\boldsymbol{\mathcal{F}}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) = \frac{1}{\Delta t} \mathit{\boldsymbol{\tilde{I}}} + J_{ \mathit{\boldsymbol{\mathcal{L}}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) - J_{ \mathit{\boldsymbol{\mathcal{S}}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}) \, , \end{equation}$

(5.4)

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]

$\begin{equation} J_{i, j} \approx \hat{J}_{i, j} = \frac{ \mathit{\boldsymbol{\mathcal{F}}}_i( \vec{\mathit {\mathit{\boldsymbol{w}}}} + \delta_j \mathit{\boldsymbol{e}}_j) - \mathit{\boldsymbol{\mathcal{F}}}_i( \vec{\mathit {\mathit{\boldsymbol{w}}}})}{\delta_j} \, . \end{equation}$

(6.1)

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.

Table 1. Execution times for computing the exact Jacobian of

$\mathit{\boldsymbol{\mathcal{F}}}$ versus execution times computing an approximate Jacobian using the finite difference approximation (6.1). Even at low resolutions, using difference approximations to compute the Jacobian is clearly unrealistic.

Resolution	Exact	FD approximation
$5 \times 5$	0.005s	0.858s
$10 \times 10$	0.006s	13.85s
$15 \times 15$	0.006s	76.49s
$20 \times 20$	0.007s	261.6s

| Show Table

DownLoad: CSV

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

$\begin{equation} \begin{aligned} u & = 1 + 0.1 \sin(3\pi x-0.01t)\sin(3\pi y-0.01t) \\ v & = \sin(3\pi x-0.01t)\sin(3\pi y-0.01t) \\ p & = \cos(3\pi x-0.01t)\cos(3\pi y-0.01t) \, . \end{aligned} \end{equation}$

(6.2)

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.

Table 2. Error and convergence rate.

operator	SBP21		SBP42
$N$	$\\| \mathit{\boldsymbol{e}}\\|$	$r$	$\\| \mathit{\boldsymbol{e}}\\|$	$r$
21	4.13e-02	–	1.90e-02	–
41	9.73e-03	2.16	2.19e-03	3.23
61	4.17e-03	2.13	6.34e-04	3.12
81	2.28e-03	2.12	2.70e-04	3.01
Theoretical		2		3

| Show Table

DownLoad: CSV

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{L}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^*) = \mathit{\boldsymbol{\mathcal{S}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^*) \, . \end{equation}$

(6.3)

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

$\begin{equation} \mathit{\boldsymbol{\mathcal{F}}} ( \vec{\mathit {\mathit{\boldsymbol{w}}}}^*) = \mathit{\boldsymbol{\mathcal{L}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^*) - \mathit{\boldsymbol{\mathcal{S}}}( \vec{\mathit {\mathit{\boldsymbol{w}}}}^*) = 0 \, . \end{equation}$

(6.4)

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]

$\begin{equation} \begin{aligned} u & = 1-e^{\lambda x} \cos(2\pi y) & v & = \frac{1}{2 \pi}\lambda e^{\lambda x}\sin(2\pi y) \\ p & = \frac{1}{2}\left(1 - e^{2\lambda x}\right) & \lambda & = \frac{1}{2 \epsilon} - \sqrt{\frac{1}{4 \epsilon^2} + 4 \pi^2} \end{aligned} \end{equation}$

(6.5)

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}$ . shows the error and convergence rates, which again agree well with the theory. In , 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].

Table 3. Error and (accuracy) convergence rate of (6.5).

operator	SBP21		SBP42
$N$	$\\| \mathit{\boldsymbol{e}}\\|$	$r$	$\\| \mathit{\boldsymbol{e}}\\|$	$r$
21	2.04e-01	–	4.95e-02	–
41	4.56e-02	2.16	6.86e-03	2.85
61	2.04e-02	1.98	2.20e-03	2.80
81	1.16e-02	1.97	9.76e-04	2.83
101	7.46e-03	1.97	5.16e-04	2.85
Theoretical		2		3

| Show Table

DownLoad: CSV

Figure 2. Streamlines and the velocity field of (6.5).

DownLoad: Full-Size Img PowerPoint

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

$\frac{e_{k+1}}{e_{k}} \approx \left(\frac{e_{k}}{e_{k-1}}\right)^p$

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

$p \approx \frac{\log(e_{k+1}/e_k)}{\log(e_{k}/e_{k-1})} \, .$

The error $e_k = \| \vec{\mathit {\mathit{\boldsymbol{w}}}}^k- \vec{\mathit {\mathit{\boldsymbol{w}}}}^*\|_\infty$ is presented in 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.

Table 4. Errors and the estimated (iterative) convergence rates of (6.5).

$k$	$\\|e_k\\|_\infty$	$p$
1	3.56e+00	–
2	1.85e+01	–
3	1.89e+00	-1.38
4	1.21e+00	0.20
5	5.53e-01	1.74
6	1.10e-01	2.07
7	3.21e-03	2.19
8	3.31e-06	1.94
9	4.14e-12	1.98
Theoretical		2

| Show Table

DownLoad: CSV

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. 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.

Figure 3. Viscous flow over a solid surface. The plot illustrates

$\mathit{\boldsymbol{u}}$ for the converged solution.

DownLoad: Full-Size Img PowerPoint

Table 5. Errors and the estimated (iterative) convergence rates for the flow over a solid surface.

$k$	$\\|e_k\\|_\infty$	$p$
1	1.68e+00	–
2	6.17e-01	–
3	1.15e-01	1.67
4	4.37e-03	1.95
5	1.02e-05	1.85
6	5.51e-11	2.00
Theoretical		2

| Show Table

DownLoad: CSV

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 , the converged solution is illustrated and the estimated iterative convergence rate $p$ is presented in 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.

Figure 4. Inviscid flow over a smooth bump. The plot illustrates the velocity field (arrows) and

$\mathit{\boldsymbol{u}}$ (color figure) at the converged solution.

DownLoad: Full-Size Img PowerPoint

Table 6. Errors and the estimated (iterative) convergence rates for the bump.

$k$	$\\|e_k\\|_\infty$	$p$
1	3.56e+00	–
2	4.91e-01	–
3	1.74e-02	1.69
4	3.33e-05	1.87
5	1.55e-10	1.96
Theoretical		2

| Show Table

DownLoad: CSV

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.

References

[1]	Taylor K (2012) Landscape and meaning: Context for a global discourse on cultural landscape values. Managing Cultural Landscapes. Available from: http://hdl.handle.net/1885/59693.
[2]	Taylor K (2014) Cities as cultural landscapes. In: Bandarin, F., Van Oers, R., Reconnecting the City: The Historic Urban Landscape Approach and the Future of Urban Heritage, 172–202. https://doi.org/10.1002/9781118383940.ch7
[3]	Taylor K, Clair AS, Mitchell NJ (2015) Conserving Cultural Landscapes: Challenges and New Directions. New York: Routledge.
[4]	Oevermann H, Mieg HA (2015) Industrial Heritage Sites in Transformation. New York: Routledge.
[5]	Collavitti AM (2018) Urban Heritage Management. Switzerland: Springer, Cham. https://doi.org/10.1007/978-3-319-72338-9
[6]	Ogawa Y (2019) The Memory Police. New York: Pantheon.
[7]	Assman A (2011) Cultural Memory and Western Civilization. Cambridge: Cambridge University Press.
[8]	Macdonald S (2013) Memorylands. London: Routledge.
[9]	Czepczynski M (2012) Cultural Landscapes of Post-Socialist Countries. Dorchester: Ashgate.
[10]	Kubicki P, Mach E (2022) European Cities in the Process of Constructing and Transmitting European Cultural Heritage. Kraków: Księgarnia Akademicka Publishing. https://doi.org/10.12797/9788381386708
[11]	Kostialova K (2015) The main aspects of the transformation of town disctricts—Podborova housing estate. Kulturna a socialna diverzita na Slovensku/Cultural and Social Diversity in Slovakia. Banska Bystrica: Belianum.

Reader Comments

Your name:*

Email:*
© 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Urban Resilience and Sustainability

Metrics

Article views(1305) PDF downloads(34) Cited by(0)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(2)

Urban Resilience and Sustainability

Speaking territories in the cultural landscape: Challenges of transformation in the Central Europe

Related Papers:

Abstract

1. Introduction

2. Problem formulation

3. The semi-discrete scheme

4. Exact computation of the Jacobian

4.1. The Jacobian of the spatial operator

4.2. The Jacobian of the penalty terms

5. The fully discrete scheme

6. Numerical experiments

6.1. The order of accuracy

6.2. The convergence rate of the Newton iteration

7. Conclusions

Author contributions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Urban Resilience and Sustainability

Speaking territories in the cultural landscape: Challenges of transformation in the Central Europe

Related Papers:

Abstract

1. Introduction

2. Problem formulation

3. The semi-discrete scheme

4. Exact computation of the Jacobian

4.1. The Jacobian of the spatial operator

4.2. The Jacobian of the penalty terms

5. The fully discrete scheme

6. Numerical experiments

6.1. The order of accuracy

6.2. The convergence rate of the Newton iteration

7. Conclusions

Author contributions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog