1 Introduction

The compressible Navier–Stokes equations are the foundation of computational fluid dynamics (CFD) for modelling the flow of viscous compressible fluids. Consequently, numerical methods for approximating their solutions are vastly studied. For a numerical scheme to yield a convergent sequence of approximate solutions, it must be a stable discretisation of the well-posed continuous problem. For linear partial differential equations (PDEs), the energy method, which depends heavily on integration by parts (IBP), is often used to prove well-posedness. In the (semi-)discrete setting, analogous stability proofs can be obtained by using the discrete energy method, where IBP is mimicked using summation-by-parts (SBP). Numerical methods formulated to satisfy the SBP property are thus frequently used for various PDEs, including CFD problems (see e.g. [6, 9, 24, 31, 32]).

Different numerical methods can be formulated in the SBP framework. These include the finite-difference methods (see e.g. [20, 21, 28]), the finite-volume methods (see e.g. [7, 23, 29]) and the discontinuous Galerkin spectral element methods (see e.g. [12]). The latter two may be formulated on unstructured grids, that are sometimes preferred for domains with complex geometries.

Herein, we focus the attention on local finite-volume methods that only use nearest neighbours to approximate the derivatives. These are still the workhorse methods in production CFD, due to their simplicity and robustness, and since the local structure allows for easier parallelisation. A well-known drawback is however the difficulty of finding consistent second-derivative approximations, which hampers their usability for the compressible Navier–Stokes equations. It was, for example, shown in [29] that a commonly used edge-based approximation is inconsistent on general unstructured grids. Although proofs of convergence exist for finite-volume methods, they often rely on admissible meshes (in the sense of Def. 3.1 in [11], see e.g. [3, 13]), that require normal derivative approximations at volume faces to be orthogonal to the face. This severely constrains mesh generation. Hence, it is desirable to design a local finite-volume scheme that runs on standard unstructured grids such as Delaunay triangulations.

In the interest of accurately discretising the viscous terms of the compressible Navier–Stokes equations on such grids, we study the Laplacian approximation proposed by Chandrashekar in [7]. His approximation incorporated the Dirichlet boundary conditions weakly, and the resulting operator was shown to satisfy the SBP property. The approximation was then used to discretise the heat equation, and numerical experiments showed that the scheme converged with second-order rate on triangulated grids.

In this work, we slightly modify the Laplacian operator from [7], by not including any boundary conditions directly in the operator. We mimic the proof of Chandrashekar, and demonstrate that the modified operator maintain the SBP property proved in [7]. To study the consistency and convergence of the Laplacian approximation, we use the heat equation as a model equation. We propose a numerical scheme for this equation where the Dirichlet boundary conditions are imposed strongly by injection (see e.g. [16] for more information about this technique), and the Neumann and Robin boundary coditions are imposed weakly similar to [7]. This approach is analogous to the one used in [15] to prove both linear and non-linear stability for the compressible Navier–Stokes equations augmented with the no-slip adiabatic wall boundary conditions on structured grids.

The main goal herein is to mathematically prove the convergence of the proposed scheme for the heat equation, thus also proving the consistency, in a weak sense, of the second-derivative approximation. By utilising the SBP properties of the Laplacian operator, we find a priori \(H^1\) estimates for the numerical solution. These estimates guarantee that the numerical solution converges weakly (up to a subsequence) to a weak solution of the heat equation. Furthermore, we show that the numerical solution converges strongly by employing Aubin–Lions’ lemma, and subsequently show that the weak solution is unique. The present proof is valid for general triangular grids with Lipschitz boundaries, and does not require admissible meshes. By using the method of manufactured solutions, we verify by a numerical experiment that the scheme is convergent.

Remark 1.1

To the best of our knowledge, this is the first convergence proof for a local finite-volume method for the second-derivative that does not require admissible meshes. We note that some multi-point flux approximations (MPFA) finite-volume methods have been proven convergent by identifying them as mixed finite-element approximations (see e.g. [2, 18]).

The proof presented herein is also easily adapted to weakly imposed boundary conditions. Stability for such a scheme for the heat equation was established in [7], and herein we show that injected Dirichlet boundary conditions also yield a stable scheme. That is, both approaches are applicable, and we have chosen the strong imposition to provide an alternative.

The paper is further organised as follows. Section 2 defines the problem, whereas a priori estimates for the continuous solution are found in Sect. 3. In Sect. 4 we state the weak formulation of the problem. Section 5 concerns the spatial discretisation and provides the proof of the slightly altered Laplacian operator being SBP. In Sect. 6, the numerical scheme that approximates our problem is stated. Furthermore, the SBP properties of the Laplacian operator are utilised to obtain discrete a priori estimates similar to those found for the continuous solution. Using these estimates, we show in Sect. 7 that the approximate solution obtained by the proposed numerical scheme converges weakly to a weak solution of the original problem. Furthermore, we show in Sect. 8 that the solution indeed converges strongly by using Aubin–Lions’ lemma. The solution is subsequently shown to be unique in Sect. 9. Finally, Sect. 10 provides a numerical example that demonstrates the convergence of the scheme.

2 Problem Statement

Consider the heat equation on a two-dimensional open polygonal Lipschitz domain, \(\Omega \), with boundary \(\partial \Omega \):

$$\begin{aligned} \begin{aligned} v_t&= \nabla \cdot (\mu \nabla v),&\text {on } \Omega \times (0,T], \\ v&= g^{\scriptscriptstyle D},&\text {on } \partial \Omega ^D \times [0,T], \\ \mu \nabla v \cdot \varvec{n}&= g^{\scriptscriptstyle N},&\text {on } \partial \Omega ^N \times [0,T],\\ \mu \nabla v \cdot \varvec{n}+ \alpha v&= g^{\scriptscriptstyle R},&\text {on } \partial \Omega ^R \times [0,T], \\ v |_{t=0}&= f,&\text {on } \Omega . \end{aligned} \end{aligned}$$
(1)

The superscripts D, N, R indicate the Dirichlet, Neumann and Robin parts of the boundary with corresponding boundary data. We assume \(\partial \Omega ^D \cup \partial \Omega ^N \cup \partial \Omega ^R = \partial \Omega \), and \(\partial \Omega ^D \cap \partial \Omega ^N = \partial \Omega ^D \cap \partial \Omega ^R = \partial \Omega ^N \cap \partial \Omega ^R = \emptyset \). Furthermore, \(\varvec{n}\) denotes the outward unit normal vector, \(f \in L^2(\Omega )\) is the initial data, and \(\mu > 0\), \(\alpha \ge 0\) are constants. We take \(g^{\scriptscriptstyle D} \in H^1(0,T; H^{1/2}(\partial \Omega ^D))\) and \(g^{\scriptscriptstyle N,R} \in L^2(0,T;L^2(\partial \Omega ^{N,R}))\).

To simplify the forthcoming analysis, we define a function, w, such that \(w \in L^2(0,T;H^1(\Omega ))\) and \(w_t \in L^2(0,T;H^1(\Omega ))\), and \(w|_{\partial \Omega ^D} = g^{\scriptscriptstyle D}\) (in the sense of traces). By the trace theorem, we know there exists such a \(w \in H^1(\Omega )\) (see [1]). Lastly, we choose w to satisfy \(w|_{t=0} = f\), and

$$\begin{aligned} \begin{aligned} \mu \nabla w \cdot \varvec{n}&= 0&\text {on } \partial \Omega ^N, \\ \mu \nabla w \cdot \varvec{n}+ \alpha w&= 0&\text {on } \partial \Omega ^R. \\ \end{aligned} \end{aligned}$$
(2)

Then, by introducing \(u = v-w\) (see e.g. [1, 17]), (1) can be recast to

$$\begin{aligned} u_t&= \nabla \cdot (\mu \nabla u) + F,&\text {on } \Omega \times (0,T], \end{aligned}$$
(3a)
$$\begin{aligned} u&= 0,&\text {on } \partial \Omega ^D \times [0,T], \end{aligned}$$
(3b)
$$\begin{aligned} \mu \nabla u \cdot \varvec{n}&= g^{\scriptscriptstyle N},&\text {on } \partial \Omega ^N \times [0,T], \end{aligned}$$
(3c)
$$\begin{aligned} \mu \nabla u \cdot \varvec{n}+ \alpha u&= g^{\scriptscriptstyle R},&\text {on } \partial \Omega ^R \times [0,T], \end{aligned}$$
(3d)
$$\begin{aligned} u |_{t=0}&= 0,&\text {on } \Omega . \end{aligned}$$
(3e)

Here,

$$\begin{aligned} F = \nabla \cdot (\mu \nabla w) - w_t, \end{aligned}$$
(4)

is a forcing function.

Remark 2.1

We could have made all boundary conditions homogeneous by defining w differently. However, we choose non-zero Neumann and Robin data to keep the regularity assumptions on the boundary data to a minimum.

3 A Priori Estimates for the Continuous Problem

To obtain a priori estimates on u, we use the energy method (see e.g. [17]). By inserting F given in (4) into (3a) and integrating by parts, we obtain

$$\begin{aligned} \int _\Omega uu_t \, d\textbf{x}&= \int _\Omega u \left( \nabla \cdot (\mu \nabla u) + \nabla \cdot (\mu \nabla w) - w_t \right) \, d\textbf{x}, \\ \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&= - \int _\Omega \left( \mu \nabla u \cdot \left( \nabla u + \nabla w \right) + uw_t \right) \, d\textbf{x}\\ {}&\quad + \int _{\partial \Omega } u (\mu \nabla u \cdot \varvec{n}+ \mu \nabla w \cdot \varvec{n}) \, ds, \\ {}&= - \mu \Vert \nabla u\Vert ^2_{L^2(\Omega )} - \int _\Omega \left( \mu \nabla u \cdot \nabla w + uw_t \right) \, d\textbf{x}\\ {}&\quad + \int _{\partial \Omega } u (\mu \nabla u \cdot \varvec{n}+ \mu \nabla w \cdot \varvec{n}) \, ds. \end{aligned}$$

Using Cauchy–Schwarz’s and Young’s inequality on the first integral on the right-hand side, we obtain

$$\begin{aligned} \begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&\le - \mu \Vert \nabla u\Vert ^2_{L^2(\Omega )} + \tfrac{\epsilon }{2} \mu \Vert \nabla u\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\epsilon } \mu \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} \\ {}&\quad + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} \\&\quad + \int _{\partial \Omega } u \left( \mu \nabla u \cdot \varvec{n}+ \mu \nabla w \cdot \varvec{n}\right) \, ds. \end{aligned} \end{aligned}$$
(5)

By choosing \(\epsilon = 1\), the term \( - \mu \Vert \nabla u\Vert ^2_{L^2(\Omega )} + \tfrac{\epsilon }{2} \mu \Vert \nabla u\Vert ^2_{L^2(\Omega )} = -\tfrac{\mu }{2} \Vert \nabla u\Vert ^2_{L^2(\Omega )}\). Since \(\epsilon \) is now determined, \(\tfrac{1}{2 \epsilon } \mu \Vert \nabla w\Vert ^2_{L^2(\Omega )} = \tfrac{\mu }{2}\Vert \nabla w\Vert ^2_{L^2(\Omega )}\), which is bounded by definition. Hence, (5) reads

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&\le - \tfrac{\mu }{2}\Vert \nabla u\Vert ^2_{L^2(\Omega )} + \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} \\ {}&\quad + \int _{\partial \Omega } u \left( \mu \nabla u \cdot \varvec{n}+ \mu \nabla w \cdot \varvec{n}\right) \, ds. \end{aligned}$$

Inserting the boundary conditions for w and u given in (2) and (3b)–(3d), respectively, we obtain

$$\begin{aligned} \begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot ,t)\Vert ^2_{L^2(\Omega )}&\le - \tfrac{\mu }{2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} + \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} \\&\quad + \int _{\partial \Omega ^N} \underline{u g^{\scriptscriptstyle N}} \, ds + \int _{\partial \Omega ^R} \left( u \left( \underline{g^{\scriptscriptstyle R}} - \alpha u\right) \underline{-\alpha uw} \right) \, ds. \end{aligned} \end{aligned}$$
(6)

Consider the underlined boundary terms above. We follow [19], and bound these terms by first using the Cauchy–Schwarz inequality:

$$\begin{aligned} \begin{aligned} \int _{\partial \Omega ^N} \!\!\!\!\!\! u g^{\scriptscriptstyle N} \, ds + \int _{\partial \Omega ^R} \!\!\!\!\!\! u g^{\scriptscriptstyle R} - \alpha u w \, ds&\le \Vert u\Vert _{L^2(\partial \Omega ^N)} \Vert g^{\scriptscriptstyle N}\Vert _{L^2(\partial \Omega ^N)} + \Vert u\Vert _{L^2(\partial \Omega ^R)} \Vert g^{\scriptscriptstyle R}\Vert _{L^2(\partial \Omega ^R)} \\ {}&\quad + \alpha \Vert u\Vert _{L^2(\partial \Omega ^R)} \Vert w\Vert _{L^2(\partial \Omega ^R)}, \end{aligned} \end{aligned}$$
(7)

and then by using the trace theorem, which states that \(\Vert u\Vert _{L^2(\partial \Omega )} \le C \Vert u\Vert _{H^1(\Omega )}\), \(C>0\) (see e.g. [1]):

$$\begin{aligned}&\Vert u\Vert _{L^2(\partial \Omega ^N)} \Vert g^{\scriptscriptstyle N}\Vert _{L^2(\partial \Omega ^N)}+ \Vert u\Vert _{L^2(\partial \Omega ^R)} \left( \Vert g^{\scriptscriptstyle R}\Vert _{L^2(\partial \Omega ^R)} + \alpha \Vert w\Vert _{L^2(\partial \Omega ^R)} \right) \\&\quad \quad \lesssim \Vert u\Vert _{H^1(\Omega )} \left( \Vert g^{\scriptscriptstyle N}\Vert _{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert _{L^2(\partial \Omega ^R)} + \alpha \Vert w\Vert _{H^1(\Omega )} \right) . \end{aligned}$$

Here, we have introduced the notation \(a \lesssim b\) for \(a \le Cb\), where \(C > 0\) is a constant.

By employing Young’s inequality, the boundary terms (7) finally read

$$\begin{aligned}&\int _{\partial \Omega ^N} \!\!\!\!\!\! u g^{\scriptscriptstyle N} \, ds + \int _{\partial \Omega ^R} \!\!\!\!\!\! u g^{\scriptscriptstyle R} - \alpha u w \, ds \lesssim \tfrac{\beta }{2} \Vert u\Vert ^2_{H^1(\Omega )} \\ {}&\quad + \tfrac{1}{2\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )} \right) . \end{aligned}$$

The preliminary estimate (6), can then be stated as

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&\lesssim - \tfrac{\mu }{2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} + \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} \\ {}&\quad + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} + \tfrac{\beta }{2} \Vert u\Vert ^2_{H^1(\Omega )} \\&\quad +\tfrac{1}{2\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )}\right) - \int _{\partial \Omega ^R} \alpha u^2 \, ds. \end{aligned}$$

The last term on the right-hand side is negative semi-definite, since \(\alpha \ge 0\). We neglect it in the remaining analysis. Hence we have

$$\begin{aligned} \begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&+ \tfrac{\mu }{2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} - \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} - \tfrac{\beta }{2} \Vert u\Vert ^2_{H^1(\Omega )} \\ {}&\lesssim \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} \\ {}&\quad + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )} \right) . \end{aligned} \end{aligned}$$
(8)

Consider the three last terms on the left-hand side of the above inequality. By adding and subtracting \(\tfrac{\mu }{2} \Vert u\Vert ^2_{L^2(\Omega )}\), they can be rewritten as

$$\begin{aligned}&\tfrac{\mu }{2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} - \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} - \tfrac{\beta }{2} \Vert u\Vert ^2_{H^1(\Omega )} + \tfrac{\mu }{2} \Vert u\Vert ^2_{L^2(\Omega )} - \tfrac{\mu }{2} \Vert u\Vert ^2_{L^2(\Omega )} = \tfrac{\mu - \beta }{2}\Vert \nabla u\Vert ^2_{H^1(\Omega )} \nonumber \\ {}&\quad - \tfrac{\mu + \delta }{2} \Vert u\Vert ^2_{L^2(\Omega )}. \end{aligned}$$
(9)

By choosing \(\beta \) sufficiently small \(\tfrac{\mu - \beta }{2} \Vert u\Vert ^2_{H^1(\Omega )} \ge \tfrac{\mu - \beta }{2} \Vert u\Vert ^2_{L^2(\Omega )} \ge 0\). From (8) we then have

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&\lesssim \tfrac{\beta + \delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} + \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} \\ {}&\quad + \tfrac{1}{2\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )} \right) . \end{aligned}$$

Employing Grönwall’s inequality (see e.g. [10]), we obtain

$$\begin{aligned} \Vert u(\cdot , \cdot , t)\Vert ^2_{L^2(\Omega )}&\lesssim e^{(\beta + \delta )t} \left( \Vert u(\cdot , \cdot , 0)\Vert ^2_{L^2(\Omega )} + \int _0^T \left( \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{1}{\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} \right) \, dt \right. \\ {}&\left. \quad + \int _0^T \tfrac{1}{\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )} \right) \, dt \right) . \end{aligned}$$

The inequality holds for all \(0 \le t \le T\). In Sect. 2 we defined \(w \in L^2(0,T;H^1(\Omega ))\), \(w_t \in L^2(0,T;H^1(\Omega ))\) and \(g^{\scriptscriptstyle N,R} \in L^2(0,T;L^2(\partial \Omega ^{N,R}))\). Using this, we find that the right-hand side of the above inequality is bounded. Thus, \(u \in L^\infty (0,T;L^2(\Omega ))\). Lastly, we integrate (8) in time to obtain

$$\begin{aligned}&\frac{1}{2}\Vert u(\cdot , \cdot , T)\Vert ^2_{L^2(\Omega )} + \int _0^T \left( \tfrac{\mu }{2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} - \tfrac{\delta }{2} \Vert u\Vert ^2_{L^2(\Omega )} - \tfrac{\beta }{2} \Vert u\Vert ^2_{H^1(\Omega )} \right) \, dt \\&\quad \lesssim \frac{1}{2} \Vert u(\cdot , \cdot , 0)\Vert ^2_{L^2(\Omega )} \\&\qquad + \int _0^T \left( \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )}\right. \\&\qquad \left. + \tfrac{1}{2\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )} \right) \right) \, dt \\&= \int _0^T \left( \tfrac{\mu }{2} \Vert \nabla w\Vert ^2_{L^2(\Omega )} + \tfrac{1}{2\delta } \Vert w_t\Vert ^2_{L^2(\Omega )} \right. \\&\left. \qquad +\tfrac{1}{2\beta } \left( \Vert g^{\scriptscriptstyle N}\Vert ^2_{L^2(\partial \Omega ^N)} + \Vert g^{\scriptscriptstyle R}\Vert ^2_{L^2(\partial \Omega ^R)} + \alpha ^2 \Vert w\Vert ^2_{H^1(\Omega )} \right) \right) \, dt, \end{aligned}$$

where we have used \(u|_{t=0} \equiv 0\) in the last step. Since \(\int _0^T \Vert u(\cdot , \cdot , t)\Vert ^2 \, dt \le constant\), we observe from this inequality that \(\nabla u \in L^2(0,T; L^2(\Omega ))\), and thus we have \(u \in L^2(0,T;H^1(\Omega ))\).

4 Weak Formulation of the Heat Equation

Next, we derive the weak formulation of (3). Let \(H^1_{\partial \Omega ^D_0}(\Omega )\) denote the space of \(H^1\) functions vanishing at the Dirichlet boundary. Furthermore, let \(\phi \in H^1(0,T; H_{\partial \Omega ^D_0}^1(\Omega ))\) be a test function that satisfies \(\phi (\textbf{x}, T) = 0\), \(\textbf{x}\in \Omega \). Multiply (3a) by \(\phi \) and integrate over \(\Omega \).

$$\begin{aligned} \int _\Omega \phi u_t \, d\textbf{x}&= \int _\Omega \phi \nabla \cdot (\mu \nabla u) \, d\textbf{x}+ \int _\Omega \phi F \, d\textbf{x}. \end{aligned}$$
(10)

Integrating by parts and inserting the boundary conditions given in (3b)–(3d), give

$$\begin{aligned} \int _\Omega \phi u_t \, d\textbf{x}&= - \int _\Omega \nabla \phi \cdot \mu \nabla u \, d\textbf{x}+ \int _{\partial \Omega ^D} \phi (\mu \nabla u \cdot \varvec{n}) \, ds + \int _{\partial \Omega ^N} \phi g^{\scriptscriptstyle N} \, ds \\ {}&\quad + \int _{\partial \Omega ^R} \phi (g^{\scriptscriptstyle R} - \alpha u) \, ds + \int _\Omega \phi F \, d\textbf{x}, \\ {}&= - \int _\Omega \nabla \phi \cdot \mu \nabla u \, d\textbf{x}+\int _{\partial \Omega ^N} \phi g^{\scriptscriptstyle N} \, ds + \int _{\partial \Omega ^R} \phi (g^{\scriptscriptstyle R} - \alpha u) \, ds + \int _\Omega \phi F \, d\textbf{x}, \end{aligned}$$

where we have used \(\phi |_{\partial \Omega ^D} = 0\). Using \(\phi |_{t=T} = 0\), \(u|_{t=0} = 0\) and partially integrating the left-hand side in time further yields the weak form of (3):

$$\begin{aligned} \int _0^T \int _\Omega \phi _t u \, d\textbf{x}dt&= \int _0^T \int _\Omega \nabla \phi \cdot \mu \nabla u \, d\textbf{x}dt - \int _0^T \int _{\partial \Omega ^N} \phi g^{\scriptscriptstyle N} \, dsdt \nonumber \\ {}&\quad - \int _0^T \int _{\partial \Omega ^R} \phi (g^{\scriptscriptstyle R} - \alpha u) \, dsdt - \int _0^T \int _\Omega \phi F \, d\textbf{x}dt, \end{aligned}$$
(11)

where F given by (4) satisfies

$$\begin{aligned} \int _\Omega \phi F \, d\textbf{x}&= - \int _\Omega \nabla \phi \cdot \mu \nabla w \, d\textbf{x}- \int _{\partial \Omega ^R} \alpha \phi w \, ds - \int _{\Omega } \phi w_t \, d\textbf{x}. \end{aligned}$$
(12)

Remark 4.1

Since the forcing function is not the main focus of this work, we use \(\int _\Omega \phi F \, d\textbf{x}\) as short-hand notation for (12) and make comments about it where necessary.

Remark 4.2

From (12), we see that \(w \in H^1(\Omega )\) is sufficient to bound the two first integrals on the right-hand side. Furthermore, the regularity of \(w_t\) is determined by the regularity of the boundary data (see e.g. [17]). Thus, for \(\gamma (w_t) = g_t^{\scriptscriptstyle D}\) to be satisfied (where \(\gamma \) is the trace function), we must have \(w_t \in H^1(\Omega )\), and that is why we assumed that \(g^{\scriptscriptstyle D} \in H^1(0,T;H^{1/2}(\partial \Omega ^D))\) in Sect. 2.

Definition 4.3

A function u satisfying (11) is called a weak solution of the problem (3).

5 Spatial Discretisation

Let \(\bar{\Omega }_h\) be a discretisation of \(\bar{\Omega } = \Omega \cup \partial \Omega \) into non-overlapping triangles \(K_n\), \(n = 1, \ldots , N\) such that \(\bar{\Omega }_h = \cup _{n=1}^N K_n\), and such that there are no hanging nodes in \(\bar{\Omega }_h\). The grid functions are defined on the vertices of the triangles. Furthermore, subdivide \(\bar{\Omega }_h\) into a dual grid consisting of dual cells, \(V_i\), \( i = 1, \ldots , I\), such that \(\bar{\Omega }_h = \cup _{i=1}^I V_i\). The dual cells are polygons surrounding a vertex, i. A dual-volume boundary consists of straight lines drawn from the mid-point of an edge adjacent to grid point i to the centroid of the triangles adjacent to the grid point (see Fig. 1). (These are the dual volumes of the standard node-centred finite-volume method, see e.g. [22]). We introduce the notation

$$\begin{aligned} \bar{\Omega }_h^V&: \text { the set of indices for interior and boundary nodes}, \\ \bar{\Omega }_h^K&: \text { the set of triangles in} \bar{\Omega }_h, \\ \Omega _h^V&: \text { the set of indices for interior nodes}, \\ \partial \Omega _h^V&: \text { the set of indices for boundary nodes}, \\ \partial \Omega _h^N&: \text { the set of indices for boundary nodes on } \partial \Omega ^N, \\ \partial \Omega _h^R&: \text { the set of indices for boundary nodes on } \partial \Omega ^R. \end{aligned}$$
Fig. 1
figure 1

Example of a triangulation and the corresponding dual cells

Full size image

The discretisation of the problem (3) utilises the approximation of the Laplacian and gradient operator proposed in [7] for the interior scheme. For triangles having at least one edge along the Dirichlet boundary, the Dirichlet condition was incorporated weakly in the gradient operator in [7]. Here, we use the approximation for interior triangles for every triangle in the grid. The approximation is given by

$$\begin{aligned} \nabla _h u^{n} = - \frac{1}{2|K_n|} \left[ u_i \hat{\varvec{n}}^n_i + u_j \hat{\varvec{n}}^n_j + u_k \hat{\varvec{n}}^n_k \right] , \end{aligned}$$
(13)

where \(|K_n|\) is the area of triangle \(K_n\); i, j, k are the vertices of the triangle, and \(\hat{\varvec{n}}^n_{i,j,k}\) are the outward pointing normal vectors of the triangle, opposite of the particular node (see Fig. 2). The length of the normal vectors, \(\hat{\varvec{n}}^n_{i,j,k}\), is equal to the length of the adjacent edge.

Fig. 2
figure 2

Triangle depicting the components of the gradient approximation (13)

Full size image

Next, we introduce the following notation.

$$\begin{aligned} I_n&= \{\text {all vertices of triangle} n\}, \\ N_i&= \{\text {all triangles with vertex} i\}, \\ E_i&= \{\text {all boundary edges having vertex} i \text {as an endpoint}\}, \\ \end{aligned}$$

Then the approximation of the Laplacian on a dual volume is found by approximating (see [7])

$$\begin{aligned} \int _{V_i} \Delta u \, d\textbf{x}= \int _{\partial V_i \setminus \partial \Omega } \!\!\!\!\!\!\! \nabla u \cdot \varvec{n}\, ds + \int _{\partial V_i \cap \partial \Omega } \!\!\!\!\!\!\! \nabla u \cdot \varvec{n}\, ds, \end{aligned}$$
(14)

by

$$\begin{aligned} (\Delta _h u)_i = \frac{1}{V_i} \left[ \frac{1}{2} \sum _{n \in N_i} \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{e \in E_i} \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e) \right] . \end{aligned}$$
(15)

Here, \(\hat{\varvec{b}}(e)\) denotes the outward pointing normal vector at boundary edge e (see Fig. 3a). The superscript n(e) signifies the triangle that has the boundary edge e. The components of the approximation (15) is depicted in Fig. 3b.

Fig. 3
figure 3

a Example of a vertex i belonging to three triangles (\(K_1, K_4, K_7\)) where two of them (\(K_1, K_7\)) have an edge along the boundary, depicted with the corresponding boundary normals \(\hat{\varvec{b}}(e)\). b Example of a dual cell, \(V_i\), and the components of the Laplace approximation (15)

Full size image

The approximation of the Laplacian (15) with Dirichlet boundary conditions taken into account, was demonstrated to satisfy the Summation-by-Parts (SBP) property in Theorem 1 in [7]. Here, we state the analogous result without any boundary conditions.

Theorem 5.1

Let \(u^h\) and \(v^h\) be two grid functions defined on \(\bar{\Omega }_h\) such that \(u^h = (u_1, u_2, \ldots , u_I)\), and correspondingly for \(v^h\). Then the discrete approximation of the Laplacian operator (15) satisfies the SBP property

$$\begin{aligned} \sum _{i \in \bar{\Omega }_h^V} v_i V_i (\Delta _h u^h)_i&= - \sum _{n \in \bar{\Omega }_h^K} \nabla _h u^{n} \cdot \nabla _h v^{n} |K_n| + \frac{1}{2} \sum _{i \in \partial \Omega _h^V} v_i (\nabla _h u^{n_{i,1}(e)} \cdot \hat{\varvec{b}}_{i,1}(e) \\ {}&\quad + \nabla _h u^{n_{i,2}(e)} \cdot \hat{\varvec{b}}_{i,2}(e)), \end{aligned}$$

where the subscripts \(\{i,1\}\) and \(\{i,2\}\) indicate the two edges adjacent to the boundary node i.

Proof

Multiply Eq. (15) by \(v_i V_i\) and sum over all vertices in the grid.

$$\begin{aligned} \sum _{i \in \bar{\Omega }_h^V} v_i V_i (\Delta _h u^h)_i&= \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} v_i \sum _{n \in N_i} \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} v_i \sum _{e \in E_i} \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e), \nonumber \\ {}&= \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} \sum _{n \in N_i} v_i \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{i \in \partial \Omega _h^V} \sum _{e \in E_i} v_i \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e). \end{aligned}$$
(16)

In the second equality, we have used that the set \(E_i\) is empty for interior nodes.

For the first term in the above equation, we change the order of summation and move \(\nabla _h u^n\) outside the summation over the vertices of a triangle \(K_n\) in (16), to obtain

$$\begin{aligned} \begin{aligned} \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} \sum _{n \in K_i} v_i \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n}&= \frac{1}{2} \sum _{n \in \bar{\Omega }_h^K} \sum _{i \in I_n} v_i \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} = \frac{1}{2} \sum _{n \in \bar{\Omega }_h^K} \nabla _{h} u^{n} \cdot \sum _{i \in I_n} v_i \hat{\varvec{n}}_i^{n}. \end{aligned} \end{aligned}$$
(17)

For the boundary nodes, we have

$$\begin{aligned} \frac{1}{2} \sum _{i \in \partial \Omega _h^V} \sum _{e \in E_i} v_i \nabla _hu^{n(e)} \cdot \hat{\varvec{b}}(e)&= \sum _{i \in \partial \Omega _h^V} v_i (\nabla _h u^{n_{i,1}(e)} \cdot \hat{\varvec{b}}_{i,1}(e) + \nabla _h u^{n_{i,2}(e)} \cdot \hat{\varvec{b}}_{i,2}(e)) \end{aligned}$$
(18)

With (17) and (18), (16) can be written as

$$\begin{aligned} \sum _{i \in \bar{\Omega }_h^V} v_i V_i (\Delta _h u^h)_i&= \frac{1}{2} \sum _{n \in \bar{\Omega }_h^K} \nabla _{h} u^{n} \cdot \sum _{i \in I_n} v_i \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{i \in \partial \Omega _h^V} v_i (\nabla _h u^{n_{i,1}(e)} \cdot \hat{\varvec{b}}_{i,1}(e) \\ {}&\quad + \nabla _h u^{n_{i,2}(e)} \cdot \hat{\varvec{b}}_{i,2}(e)), \\ {}&= - \sum _{n \in \bar{\Omega }_h^K} \nabla _h u^{n} \cdot \nabla _h v^{n} |K_n| + \frac{1}{2} \sum _{i \in \partial \Omega _h^V} v_i (\nabla _h u^{n_{i,1}(e)} \cdot \hat{\varvec{b}}_{i,1}(e) \\ {}&\quad + \nabla _h u^{n_{i,2}(e)} \cdot \hat{\varvec{b}}_{i,2}(e)). \end{aligned}$$

In the last equality we have used the approximation of the gradient (13). \(\square \)

6 The Numerical Scheme and Discrete A Priori Estimates

To approximate the problem (1) we use (15) for the Laplacian approximation at the interior nodes. The Dirichlet condition is imposed strongly by injection (see e.g. [15, 16]). The Neumann and Robin conditions are imposed weakly in the same way as in [7]. That is, by replacing the last term of (14) with the boundary data, we approximate the Neumann and Robin boundaries by:

$$\begin{aligned} \int _{\partial V_i \cap \partial \Omega } \!\!\!\!\!\! \nabla u \cdot \varvec{n}\, ds \approx {\left\{ \begin{array}{ll} \frac{1}{2} \sum _{e \in E_i} g_i^{\scriptscriptstyle N} |\hat{\varvec{b}}(e)| &{} \text {if}\quad i\,\text {is a Neumann boundary node,} \\ \frac{1}{2}\sum _{e \in E_i} (g_i^{\scriptscriptstyle R} - \alpha u_i) |\hat{\varvec{b}}(e)| &{} \text {if}\quad i\,\text {is a Robin boundary node.} \end{array}\right. } \end{aligned}$$

Remark 6.1

Imposing the Dirichlet condition by injection means in practice that the Dirichlet nodes are overwritten by the exact boundary data after each time step. (Equivalently, no equation is solved at these nodes, since u is equal to the boundary data.)

Remark 6.2

A boundary node is either of Dirichlet, Neumann or Robin type. The entire dual-cell boundary coinciding with the physical boundary is subsequently approximated as the same type as the boundary node, see Fig. 4. This means that in the junction between two boundary types, part of the computational boundary may be approximated as something different than the actual physical boundary. However, this is an \(\mathcal {O}(h)\) error of the boundary integral which tends to zero with decreasing mesh sizes. Note that this is only necessary for the Dirichlet nodes where the boundary conditions are injected. For Neumann and Robin nodes, we could split the outer dual-cell boundary into a Neumann and a Robin part since these boundary conditions are set weakly. However, in the scheme (19) below, we use the first approach to reduce notation.

Fig. 4
figure 4

Example of a grid with corresponding dual cells with an intersection of a Neumann and Robin boundary. For boundary nodes, the whole dual cell boundary is approximated as the type of the boundary node

Full size image

The above choices lead to the following discrete approximation scheme of (1)

$$\begin{aligned} \begin{aligned} \frac{dv_i}{dt}&= \left( g_ i^{\scriptscriptstyle D}\right) _t,&i \in \partial \Omega _h^{D}, \\ \frac{dv_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h v^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2V_i} \sum _{e \in E_i} g_i^{\scriptscriptstyle N} |\hat{\varvec{b}}(e)|,&i \in \partial \Omega _h^{N}, \\ \frac{dv_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h v^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2V_i} \sum _{e \in E_i} (g_i^{\scriptscriptstyle R}- \alpha u_i) |\hat{\varvec{b}}(e)|,&i \in \partial \Omega _h^R, \\ \frac{dv_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h v^{n} \cdot \hat{\varvec{n}}_i^{n},&i \in \Omega _h, \\ v_i|_{t=0}&= f_i,&i \in \Omega _h. \end{aligned} \end{aligned}$$
(19)

Remark 6.3

For readers familiar with the simultaneous approximation term (SAT) (see e.g. the review papers [8, 30]) we remark that the schemes for the Neumann and Robin nodes are equivalent to

$$\begin{aligned} \frac{dv_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h v^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2V_i} \sum _{e \in E_i} \mu \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e) + \textsf{SAT}^{\scriptscriptstyle N}_i, \quad i \in \partial \Omega _h^{N}, \\ \frac{dv_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h v^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2V_i} \sum _{e \in E_i} \mu \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e) + \textsf{SAT}^{\scriptscriptstyle R}_i, \quad i \in \partial \Omega _h^R, \end{aligned}$$

where the SATs take the form

$$\begin{aligned} \begin{aligned} \textsf{SAT}^{\scriptscriptstyle N}_i&= - \frac{1}{2V_i} \sum _{e \in E_i} \left( \mu \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e) - g^{\scriptscriptstyle N}_i |\hat{\varvec{b}}(e)| \right) , \\ \textsf{SAT}^{\scriptscriptstyle R}_i&= - \frac{1}{2V_i} \sum _{e \in E_i} \left( \mu \nabla _h u^{n(e)} \cdot \hat{\varvec{b}}(e) - (g^{\scriptscriptstyle R}_i - \alpha u_i) |\hat{\varvec{b}}(e)| \right) . \end{aligned} \end{aligned}$$
(20)

Remark 6.4

To simplify following energy analysis, we have defined the Dirichlet nodes in (19) as \((v_i)_t = (g^{\scriptscriptstyle D}_i)_t\). We emphasise that when implementing the scheme, the Dirichlet nodes should take the form \(v_i = g^{\scriptscriptstyle D}_i\) in order to avoid discretisation errors from the time-stepping algorithm.

As for the continuous problem, we transform the scheme (19) into one that imposes homogeneous Dirichlet boundary conditions. That is, we construct a function w as defined in Sect. 2 and introduce \(u = v-w\) (see again [1, 17]). Inserting \(v = u + w \) into the scheme (19), we obtain

$$\begin{aligned} \frac{d u_i}{dt}&= 0,&i \in \partial \Omega ^D_h \end{aligned}$$
(21a)
$$\begin{aligned} \frac{d u_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^n \cdot \hat{\varvec{n}}^n_i + \frac{1}{2V_i} \sum _{e \in E_i} g^{\scriptscriptstyle N}_{i} |\hat{\varvec{b}}(e)| + F_i,&i \in \partial \Omega ^N_h \end{aligned}$$
(21b)
$$\begin{aligned} \frac{d u_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^n \cdot \hat{\varvec{n}}^n_i + \frac{1}{2V_i} \sum _{e \in E_i} (g^{\scriptscriptstyle R}_{i} - \alpha u_i) |\hat{\varvec{b}}(e)| + F_i,&i \in \partial \Omega ^R_h \end{aligned}$$
(21c)
$$\begin{aligned} \frac{d u_i}{dt}&= \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^n \cdot \hat{\varvec{n}}^n_i + F_i&i \in \Omega _h, \end{aligned}$$
(21d)
$$\begin{aligned} u_i|_{t=0}&= 0,&i \in \Omega _h, \end{aligned}$$
(21e)

where

$$\begin{aligned} F_i = \left\{ \begin{array}{ll} \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h w^n \cdot \hat{\varvec{n}}^n_i - \frac{dw_i}{dt}, &{}\quad i \in \partial \Omega ^N_h, \\ \\ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h w^n \cdot \hat{\varvec{n}}^n_i - \frac{1}{2V_i} \sum _{e \in E_i} \alpha w_i |\hat{\varvec{b}}(e)| - \frac{dw_i}{dt}, &{}\quad i \in \partial \Omega ^R_h, \\ \\ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h w^n \cdot \hat{\varvec{n}}^n_i - \frac{dw_i}{dt}, &{}\quad i \in \Omega _h. \end{array}\right. \end{aligned}$$
(22)

Remark 6.5

By the Picard–Lindelöf theorem (see e.g. [25]), the ordinary differential equation (21) has a solution if the scheme is stable.

To obtain a priori estimates for the approximate solution \(u^h = (u_1, u_2, \ldots , u_I)\), we use the discrete energy method (see e.g. [17] for more details on the energy method). That is, we multiply the scheme (21) in each node, i, by \(u_i V_i\) and sum over all grid points.

$$\begin{aligned} \sum _{i \in \bar{\Omega }_h^V} u_i V_i \frac{du_i}{dt}&= \underline{\sum _{i \in \Omega _h^V} \!\! u_i V_i \left[ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} \right] } \\ {}&\quad + \sum _{i \in \partial \Omega _h^N} u_i V_i \left[ \underline{\frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n}} + \frac{1}{2V_i} \sum _{e \in E_i} g^{\scriptscriptstyle N}_i |\hat{\varvec{b}}(e)| \right] \\ {}&\quad + \sum _{i \in \partial \Omega _h^R} u_i V_i \left[ \underline{\frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n}} + \frac{1}{2V_i} \sum _{e \in E_i} (g^{\scriptscriptstyle R}_i - \alpha u_i) |\hat{\varvec{b}}(e)| \right] \\&\quad + \sum _{i \in \bar{\Omega }_h^V} u_i V_i F_i. \end{aligned}$$

Since \( \bar{\Omega }_h = \Omega _h^V \cup \partial \Omega _h^D \cup \partial \Omega _h^N \cup \partial \Omega _h^R\), and all the sets are disjoint, and since the scheme for the Dirichlet nodes is zero, the underlined terms amount to summing over all nodes in the grid. That is, the above is equivalent to

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \sum _{i \in \bar{\Omega }_h^V} V_i u_i^2&= \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} u_i \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{i \in \partial \Omega _h^N} u_i \sum _{e \in E_i} g^{\scriptscriptstyle N}_i |\hat{\varvec{b}}(e)| \\ {}&\quad + \frac{1}{2} \sum _{i \in \partial \Omega _h^R} u_i \sum _{e \in E_i} (g^{\scriptscriptstyle R}_i - \alpha u_i ) |\hat{\varvec{b}}(e)| + \sum _{i \in \bar{\Omega }_h^V} u_i V_i F_i. \end{aligned}$$

Using Theorem 5.1, we obtain

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \sum _{i \in \bar{\Omega }_h} V_i u_i^2&= - \sum _{n \in \bar{\Omega }_h^K} \nabla _h u^{n} \cdot \mu \nabla _h u^{n} |K_n| + \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} u_i g_i^{\scriptscriptstyle N} (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \nonumber \\ {}&\quad + \sum _{i \in \partial \Omega ^R_h} \tfrac{1}{2} \left( u_ig^{\scriptscriptstyle R}_i - \alpha u_i^2 \right) (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) + \sum _{i \in \bar{\Omega }_h^V} u_iV_iF_i, \nonumber \\ {}&\le - \mu \sum _{n \in \bar{\Omega }_h^K} |\nabla _h u^{n}|^2 |K_n| + \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} u_i g_i^{\scriptscriptstyle N} (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \nonumber \\ {}&\quad + \sum _{i \in \partial \Omega ^R_h} \tfrac{1}{2} u_i(g^{\scriptscriptstyle R}_i (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \nonumber \\&\quad + \sum _{i \in \bar{\Omega }_h^V} u_iV_iF_i, \end{aligned}$$
(23)

where we in the last inequality have used that \(\sum _{i \in \partial \Omega ^R_h} -\tfrac{1}{2} \alpha u_i^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \le 0\) since \(\alpha \ge 0\). We can further manipulate the Neumann boundary terms as follows

$$\begin{aligned} \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} u_i g_i^{\scriptscriptstyle N} (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \le \sum _{i \in \partial \Omega ^N_h} |u_ig_i^{\scriptscriptstyle N}| (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|). \end{aligned}$$

Using Young’s inequality, we obtain

$$\begin{aligned} \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} u_i g_i^{\scriptscriptstyle N} (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|)&\le \tfrac{\beta }{2} \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} |u_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \\ {}&\quad + \tfrac{1}{2 \beta } \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} |g^{\scriptscriptstyle N}_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|). \end{aligned}$$

The Robin boundary terms can be manipulated the same way. Thus, (23) reads

$$\begin{aligned} \begin{aligned} \frac{1}{2} \frac{d}{dt} \sum _{i \in \bar{\Omega }_h^V} V_i u_i^2&\le - \mu \sum _{n \in \bar{\Omega }_h^K} |\nabla _h u^{n}|^2 |K_n| +\tfrac{\beta }{2} \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} |u_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \\ {}&\quad + \tfrac{1}{2 \beta } \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} |g^{\scriptscriptstyle N}_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \\ {}&\quad + \tfrac{\beta }{2} \sum _{i \in \partial \Omega ^R_h} \tfrac{1}{2} |u_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \\ {}&\quad + \tfrac{1}{2 \beta } \sum _{i \in \partial \Omega ^R_h} \tfrac{1}{2} |g^{\scriptscriptstyle R}_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) + \sum _{i \in \bar{\Omega }_h^V} u_iV_iF_i. \end{aligned} \end{aligned}$$
(24)

We introduce the following notation for the discrete equivalents of the \(L^2\)-norms:

$$\begin{aligned} \Vert u^h\Vert _{L^2_V(\Omega )}^2&= \sum _{i \in \bar{\Omega }_h^V} |u_i|^2 V_i, \end{aligned}$$
(25)
$$\begin{aligned} \Vert \nabla _h u^h\Vert _{L^2_K(\Omega )}^2&= \sum _{n \in \bar{\Omega }_h^K} |\nabla _h u^{n}|^2 |K_n|, \end{aligned}$$
(26)
$$\begin{aligned} \Vert u^h\Vert ^2_{L^2_B{(\partial \Omega _h)}}&= \sum _{i \in \partial \Omega _h^B} \tfrac{1}{2} |u_i|^2 (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|). \end{aligned}$$
(27)

Using the definitions (25)–(27), we can recast (24) as

$$\begin{aligned} \begin{aligned} \frac{d}{dt} \Vert u^h\Vert ^2_{L^2_V(\Omega )}&\le - \mu \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^N_h)} + \tfrac{1}{2 \beta } \Vert g^{{\scriptscriptstyle N}, h}\Vert ^2_{L^2_B(\partial \Omega ^N_h)} \\ {}&\quad + \tfrac{\beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} + \tfrac{1}{2 \beta } \Vert g^{{\scriptscriptstyle R}, h}\Vert ^2_{L^2_B(\partial \Omega ^R_h)} + \sum _{i \in \bar{\Omega }_h^V} u_iV_iF_i. \end{aligned} \end{aligned}$$

To obtain an estimate analogous to (8), we must consider the forcing term \(\sum _{i \in \bar{\Omega }_h^V} u_i V_i F_i\). Except for the time-derivative term in (22), \(F_i\) takes the same form as the right-hand side of the scheme (21). By using the SBP property from Theorem 5.1 and Young’s inequality, we obtain

$$\begin{aligned} \sum _{i \in \bar{\Omega }_h^V} u_i V_i F_i&\le \tfrac{\epsilon }{2} \mu \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\mu }{2 \epsilon } \Vert \nabla _h w^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} \\ {}&\quad + \tfrac{\alpha ^2}{2 \beta } \Vert w^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} + \tfrac{\delta }{2} \Vert u^h\Vert ^2_{L^2_V(\Omega _h)} + \tfrac{1}{2 \delta } \Vert w^h_t\Vert ^2_{L^2_V(\Omega )}. \end{aligned}$$

Thus, we have

$$\begin{aligned} \begin{aligned} \frac{d}{dt} \Vert u^h\Vert ^2_{L^2_V(\Omega )}&\le - \mu \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^N_h)} + \tfrac{1}{2 \beta } \Vert g^{{\scriptscriptstyle N}, h}\Vert ^2_{L^2_B(\partial \Omega ^N_h)} + \tfrac{\beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} \\ {}&\quad + \tfrac{1}{2 \beta } \Vert g^{{\scriptscriptstyle R}, h}\Vert ^2_{L^2_B(\partial \Omega ^R_h)} \\ {}&\quad + \tfrac{\epsilon }{2} \mu \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\mu }{2 \epsilon } \Vert \nabla _h w^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} \\ {}&\quad + \tfrac{\alpha ^2}{2 \beta } \Vert w^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} + \tfrac{\delta }{2} \Vert u^h\Vert ^2_{L^2_V(\Omega )} + \tfrac{1}{2 \delta } \Vert w^h_t\Vert ^2_{L^2_V(\Omega )}. \end{aligned} \end{aligned}$$
(28)

Similarly as for the continuous problem, if we choose \(\epsilon = 1\), we obtain \(- \mu \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{\epsilon }{2} \mu \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} = - \tfrac{\mu }{2} \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )}\) in (28). Furthermore, \(\beta \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^N_h)} + \frac{ \beta }{2} \Vert u^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)} \lesssim \beta \Vert u^h\Vert ^2_{L^2_B(\partial \Omega )}\). Hence, we have

$$\begin{aligned} \frac{d}{dt} \Vert u^h\Vert ^2_{L^2_V(\Omega )}&\lesssim - \tfrac{\mu }{2}\Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} +\tfrac{\mu }{2} \Vert \nabla _h w^h\Vert ^2_{L^2_K(\Omega )} \\ {}&\quad + \tfrac{\delta }{2} \Vert u^h\Vert ^2_{L^2_V(\Omega )} + \tfrac{1}{2\delta } \Vert w^h_t\Vert ^2_{L^2_V(\Omega )} + \beta \Vert u^h\Vert ^2_{L^2_B(\partial \Omega )} \\ {}&\quad \tfrac{1}{2\beta } \left( \Vert g^{{\scriptscriptstyle N}, h}\Vert ^2_{L^2_B(\partial \Omega ^N_h)} + \Vert g^{{\scriptscriptstyle R}, h}\Vert ^2_{L^2_B(\partial \Omega ^R_h)} + \Vert w^h\Vert ^2_{L^2_B(\partial \Omega ^R_h)}\right) . \end{aligned}$$

Finally using the trace theorem, we arrive at a similar estimate as in (8):

$$\begin{aligned} \frac{d}{dt} \Vert u^h\Vert ^2_{L^2_V(\Omega )}&+ \tfrac{\mu }{2} \Vert \nabla _h u^h\Vert ^2_{L^2_K(\Omega )} - \tfrac{\delta }{2} \Vert u^h\Vert ^2_{L^2_V(\Omega )} - \beta \Vert u^h\Vert ^2_{H^1_K(\Omega )} \\&\lesssim \tfrac{\mu }{2} \Vert \nabla _h w^h\Vert ^2_{L^2_K(\Omega )} + \tfrac{1}{2\delta } \Vert w^h_t\Vert ^2_{L^2_V(\Omega )} \\ {}&\quad + \tfrac{1}{2\beta } \left( \Vert g^{{\scriptscriptstyle N}, h}\Vert ^2_{L^2_B(\partial \Omega ^N_h)} + \Vert g^{{\scriptscriptstyle R}, h}\Vert ^2_{L^2_B(\partial \Omega ^R_h)} + \Vert w^h\Vert ^2_{H^1_K(\Omega )}\right) . \end{aligned}$$

Note that we have arrived at a semi-discrete equivalent of (8). Thus, by using Grönnwall’s inequality followed by integration in time, as done in Sect. 3, we obtain \(u^h \in L^\infty (0,T;L^2_V(\Omega ))\) and \(\nabla _h u^h \in L^2(0,T;L^2_K(\Omega ))\). We may extend the numerical solution, \(u^h\), to the entire domain by a linear interpolation on the triangles. Let \(u^h_c\) denote this continuous piecewise linear function. We have that \(\nabla _h u^h_c = \nabla u^h_c = \nabla _h u^h\). Hence \(\nabla u^h_c \in L^2(0,T;L^2_K(\Omega ))\) (and also, \(\nabla u^h_c \in L^2(0,T;L^2(\Omega ))\) since \(\nabla u^h_c\) is piecewise constant). Furthermore, the norm \(\Vert u^h_c\Vert ^2_{L^2(\Omega )}\) can be bounded by \(\Vert u^h\Vert ^2_{L^2_V(\Omega )}\). Thus, we have \(u^h_c \in L^2(0,T;H^1(\Omega ))\).

7 Weak Convergence to a Weak Solution

Let \(\phi \in H^1(0,T; C^\infty _{\partial \Omega ^D_0}(\bar{\Omega }) )\). Since \(\phi \) is smooth (in space), \(\phi |_{V_i}\), which is the restriction of \(\phi \) to a dual cell, can be written as \(\phi |_{V_i} = \phi (x_i,y_i,t) + h p_i = \phi _i + h p_i\), where h is a characteristic mesh size and \(p_i(x_i, y_i, t)\) is a function of size \(\mathcal {O}(1)\). The gradient approximation (13) is \(\nabla _h \phi |_{K_n} = \nabla \phi |_{K_n} + \mathcal {O}(h)\). This can easily be checked for equilateral triangles. Thereafter, one can prove the relation for a general triangle by transforming it to an equilateral one using a linear transformation.

We denote the right-hand side of the scheme (21) by \(L_hu^h\). To prove convergence to a weak solution, we test the numerical scheme (21) against \(\phi \). That is, we calculate

$$\begin{aligned} \begin{aligned} \int _\Omega \phi u^h_t \, d\textbf{x}&= \int _\Omega \phi (L_h u^h) \, d\textbf{x}, \\ {}&= \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} \phi |_{V_i} (L_h u^h)|_{V_i} \, d\textbf{x}. \end{aligned} \end{aligned}$$
(29)

We now use that \(\phi |_{V_i} = \phi _i + h p_i\) to obtain

$$\begin{aligned} \int _\Omega \phi u^h_t \, d\textbf{x}&= \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} (\phi _i + h p_i) (L_h u^h)_i \, d\textbf{x}\end{aligned}$$
(30)
$$\begin{aligned}&= \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} \phi _i (L_hu^h)_i \, d\textbf{x}+ \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} h p_i (L_h u^h)_i \, d\textbf{x}\nonumber \\ {}&= \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} \phi _i (L_h u^h)_i \, d\textbf{x}+ \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} h p_i (u_i)_t \, d\textbf{x}, \end{aligned}$$
(31)

where we have used \(u^h_t = L_h u^h\) in the last step. Thus

$$\begin{aligned} \int _\Omega (\phi - h p )u^h_t \, d\textbf{x}&= \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} \phi _i (L_h u^h)_i \, d\textbf{x}. \end{aligned}$$
(32)

Inserting the specific form of \(L_h u^h\) (that is, the right-hand side of the scheme (21)) yields

$$\begin{aligned} \begin{aligned} \int _\Omega (\phi -h p)u^h_t \, d\textbf{x}&= \sum _{i \in \Omega _h^V} \int _{V_i} \phi _i \left[ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} \right] \, d\textbf{x}\\ {}&\quad + \sum _{i \in \partial \Omega _h^N} \int _{V_i} \phi _i \left[ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2V_i} \sum _{e \in E_i} g_i^{\scriptscriptstyle N} |\hat{\varvec{b}}(e)| \right] \, d\textbf{x}\\ {}&\quad + \sum _{i \in \partial \Omega _h^R} \int _{V_i} \phi _i \left[ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2V_i} \sum _{e \in E_i} (g_i^{\scriptscriptstyle N} - \alpha u_i) |\hat{\varvec{b}}(e)| \right] \, d\textbf{x}\\ {}&\quad + \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} \phi _i F_i \, d\textbf{x}. \end{aligned} \end{aligned}$$
(33)

Since \(\phi _i\) is constant on each dual cell, \(V_i\) and the Laplacian approximation is a scalar constant, the right-hand side above can be integrated exactly, leading to

$$\begin{aligned} \begin{aligned} \int _\Omega (\phi - h p)u^h_t \, d\textbf{x}&= \underline{\sum _{i \in \Omega _h^V} \phi _i V_i \left[ \frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} \right] } \\ {}&\quad + \sum _{i \in \partial \Omega _h^N} \phi _i V_i \left[ \underline{\frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n}} + \frac{1}{2V_i} \sum _{e \in E_i} g_i^{\scriptscriptstyle N} |\hat{\varvec{b}}(e)|\right] \\ {}&\quad + \sum _{i \in \partial \Omega _h^R} \phi _i V_i \left[ \underline{\frac{1}{2V_i} \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n}} + \frac{1}{2V_i} \sum _{e \in E_i} (g_i^{\scriptscriptstyle N} - \alpha u_i) |\hat{\varvec{b}}(e)| \right] \\ {}&\quad + \sum _{i \in \bar{\Omega }_h^V} \phi _i V_i F_i. \end{aligned} \end{aligned}$$
(34)

As in the discrete analysis in Sect. 6, the underlined terms can be written as the sum over all grid points in \(\bar{\Omega }_h\) as follows

$$\begin{aligned} \begin{aligned} \int _\Omega (\phi - hp)u^h_t \, d\textbf{x}&= \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} \phi _i \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{i \in \partial \Omega _h^N} \phi _i \sum _{e \in E_i} g_i^{\scriptscriptstyle N} |\hat{\varvec{b}}(e)| \\ {}&\quad + \frac{1}{2} \sum _{i \in \partial \Omega _h^R} \phi _i \sum _{e \in E_i} (g_i^{\scriptscriptstyle R} - \alpha u_i) |\hat{\varvec{b}}(e)| + \sum _{i \in \bar{\Omega }_h^V} \phi _i V_i F_i. \end{aligned} \end{aligned}$$
(35)

Using the SBP properties from Theorem 5.1 yields

$$\begin{aligned} \begin{aligned} \int _\Omega (\phi -hp)u^h_t \, d\textbf{x}&= - \sum _{n \in \bar{\Omega }_h^K} \nabla _h \phi ^{n} \cdot \mu \nabla _h u^{n} |K_n| + \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} \phi _i g_i^{\scriptscriptstyle N} ( |\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)| )\\ {}&\quad + \sum _{i \in \partial \Omega ^R_h} \tfrac{1}{2} \phi _i (g_i^{\scriptscriptstyle R} - \alpha u_i) ( |\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)| )+ \sum _{i \in \bar{\Omega }_h^V} \phi _i V_i F_i. \end{aligned} \end{aligned}$$
(36)

Since \(\nabla _h \phi ^n\) and \(\nabla _h u^n\) are constant on each triangle K, we have that \(- \sum _{n \in \bar{\Omega }_h^K} \nabla _h \phi ^n \cdot \mu \nabla _h u^n |K_n| = - \sum _{n \in \bar{\Omega }_h^K} \int _{K_n} \nabla _h \phi ^h \cdot \mu \nabla _h u^h \, d\textbf{x}\). Thus, (36) can be written as

$$\begin{aligned} \int _\Omega (\phi - hp)u^h_t \, d\textbf{x}&= - \sum _{n \in \bar{\Omega }_h^K} \int _{K_n} \nabla _h \phi ^n \cdot \mu \nabla _h u^n \, d\textbf{x}\\ {}&\quad + \sum _{i \in \partial \Omega ^N_h} \tfrac{1}{2} \phi _i g_i^{\scriptscriptstyle N} (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,1}(e)|) \\ {}&\quad + \sum _{i \in \partial \Omega ^R_h} \tfrac{1}{2} \phi _i (g_i^{\scriptscriptstyle R} - \alpha u_i) (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,1}(e)|) \\ {}&\quad + \sum _{i \in \bar{\Omega }_h^V} \int _{V_i} \phi _i F_i \, d\textbf{x},\\ {}&= - \int _{\Omega } \nabla _h \phi ^h \cdot \mu \nabla _h u^h \, d\textbf{x}\\ {}&\quad + \int _\Omega \phi ^h F^h \, d\textbf{x}+ \int _{\partial \Omega ^N_h} \phi ^h g^{{\scriptscriptstyle N},h} \, ds + \int _{\partial \Omega ^R_h} \phi ^h (g^{{\scriptscriptstyle R}, h} - \alpha u^h) \, ds. \end{aligned}$$

Partial integration in time yields

$$\begin{aligned} \begin{aligned} \int _0^T \int _\Omega (\phi - h p)_t u^h \, d\textbf{x}dt&= \int _0^T \int _\Omega \nabla _h \phi ^h \cdot \mu \nabla _h u^h \, d\textbf{x}dt - \int _0^T \int _{\partial \Omega ^N_h} \phi ^h g^{{\scriptscriptstyle N}, h} \, ds dt \\ {}&\quad - \int _0^T \int _{\partial \Omega ^R_h} \phi ^h (g^{{\scriptscriptstyle R}, h} - \alpha u^h) \, ds dt - \int _0^T \int _\Omega \phi ^h F^h \, d\textbf{x}dt \\ {}&\quad - \int _{\Omega } h p u^h(T) \, d\textbf{x}, \end{aligned} \end{aligned}$$
(37)

where we have used \(u^h|_{t=0} = 0\) and \(\phi ^h|_{t=T} = 0\).

Remark 7.1

Here, \(\int _\Omega \phi ^h F^h \, d\textbf{x}\) is the short-hand for the semi-discrete form of (12). By using the SBP property from Theorem 5.1, it can be written as

$$\begin{aligned} \begin{aligned} \int _\Omega \phi ^h F^h \, d\textbf{x}&= \sum _{i \in \bar{\Omega }_h^V} \phi _i V_i F_i, \\ {}&= - \sum _{n \in \bar{\Omega }_h^K} \nabla _h \phi ^n \cdot \mu \nabla _h w^n |K_n| - \sum _{i \in \partial \Omega _h^R} \alpha \phi _i w_i (|\hat{\varvec{b}}_{i,1}(e)| + |\hat{\varvec{b}}_{i,2}(e)|) \\ {}&\quad \quad - \sum _{i \in \bar{\Omega }_h^V} \phi _i V_i\frac{d w_i}{dt}, \\ {}&= - \int _\Omega \nabla _h \phi ^h \cdot \mu \nabla _h w^h \, d\textbf{x}- \int _{\partial \Omega ^R_h} \alpha \phi ^h w^h \, ds - \int _\Omega \phi ^h w^h_t \, ds. \end{aligned} \end{aligned}$$
(38)

We keep the symbolic expression to reduce notation.

Since \(\phi |_{V_i} = \phi ^h|_{V_i} + hp_i\) and \(\nabla _h \phi |_{K_n} = \nabla \phi + \mathcal {O}(h)\), the weak formulation (37) becomes

$$\begin{aligned} \begin{aligned} \int _0^T \int _\Omega (\phi _t - hp_t) u^h \, d\textbf{x}dt&= \int _0^T \int _\Omega (\nabla \phi + \mathcal {O}(h)) \cdot \mu \nabla _h u^h \, d\textbf{x}dt \\ {}&\quad - \int _0^T \int _{\partial \Omega ^N_h} (\phi - hp) g^{{\scriptscriptstyle N}, h} \, ds dt\\ {}&\quad - \int _0^T \int _{\partial \Omega ^R_h} (\phi - hp) (g^{{\scriptscriptstyle R}, h} - \alpha u^h) \, ds dt\\ {}&\quad - \int _0^T \int _\Omega (\phi - hp) F^h \, d\textbf{x}dt \\ {}&\quad - \int _{\Omega } h p u^h(T) \, d\textbf{x}. \end{aligned} \end{aligned}$$
(39)

We utilise the following functional analysis theorem (see e.g. [10], and [5] for a proof).

Theorem 7.2

Let \(\Omega _T \subset \mathbb {R}^n\) be an open domain and let \(\{u_n\} \in L^2(\Omega _T)\) be a bounded sequence. Then there exists a subsequence, \(\{u_{n_i}\} \in L^2(\Omega _T)\) that converges weakly to \(\bar{u} \in L^2(\Omega _T)\). That is,

$$\begin{aligned} \int _{\Omega _T} \phi u_{n_i} \, d\textbf{x}\rightarrow \int _{\Omega _T} \phi \bar{u} \, d\textbf{x}\quad \text {as} n_i \rightarrow \infty , \text {for all} \phi \in L^2(\Omega _T). \end{aligned}$$

Here, we take \(\Omega _T = \Omega \times [0,T]\). Consider the \(\mathcal {O}\)(1) term on the left-hand side of (39). Using Theorem 7.2, we have that

$$\begin{aligned} \int _0^T \int _\Omega \phi _t u^h \, d\textbf{x}dt&\rightarrow \int _0^T \int _\Omega \phi _t \bar{u} \, d\textbf{x}dt. \end{aligned}$$

The other \(\mathcal {O}(1)\) terms can be treated in a similar way. Turning to the second term in (39), we have

$$\begin{aligned} \int _0^T \int _\Omega hp_tu^h \, d\textbf{x}dt \rightarrow 0, \end{aligned}$$

as \(h \rightarrow 0\), since \(u^h \in L^\infty (0,T;L^2(\Omega ))\). Using the available bounds, similar arguments imply that \(\mathcal {O}(h) \mu \nabla _h u^h\), \(hpg^{{\scriptscriptstyle N},h}\), \(hpg^{{\scriptscriptstyle R},h}\), \(hpu^h\) and \(hpF^h\) vanish.

Remark 7.3

Since all terms in F (see (38)) are known and bounded in \(L^2(\Omega _T)\) (see the assumptions in Sect. 2), the weak convergence of the symbolic expression (38) follows trivially.

In summary, letting \(h \rightarrow 0\), (39) becomes

$$\begin{aligned} \int _0^T \int _\Omega \phi _t \bar{u} \, d\textbf{x}dt&= \int _0^T \int _\Omega \nabla \phi \cdot \mu \overline{\nabla u} \, d\textbf{x}dt \nonumber \\ {}&\quad - \int _0^T \int _{\partial \Omega ^N} \!\!\!\! \phi \bar{g}^{\scriptscriptstyle N} \, ds dt- \int _0^T \int _{\partial \Omega } \phi (\bar{g}^{\scriptscriptstyle R} - \alpha \bar{u}) \, ds dt \nonumber \\&\quad - \int _0^T \int _\Omega \phi \bar{F} \, d\textbf{x}dt, \end{aligned}$$
(40)

which is satisfied for all \(\phi \in H^1(0,T;C^\infty _{\partial \Omega _0^D}(\bar{\Omega }))\).

Remark 7.4

Note that the boundary integrals over the computational boundaries converge to the integrals over the physical boundaries as \(h \rightarrow 0\). That is,

$$\begin{aligned}&\int _{\partial \Omega ^N_h} \phi ^h g^{{\scriptscriptstyle N},h} \, ds \rightarrow \int _{\partial \Omega ^N} \phi ^h g^{{\scriptscriptstyle N},h} \, ds \text { and } \\ {}&\int _{\partial \Omega ^R_h} \phi ^h (g^{{\scriptscriptstyle R},h} - \alpha u^h) \, ds \rightarrow \int _{\partial \Omega ^R} \phi ^h (g^{{\scriptscriptstyle R},h} - \alpha u^h) \, ds, \end{aligned}$$

as \(h \rightarrow 0\).

Remark 7.5

The term \(\int _\Omega \phi _t u^h \, d\textbf{x}\) in (39) satisfies

$$\begin{aligned} \int _\Omega \phi _t u^h \, d\textbf{x}= \int _\Omega \phi _t u^h_c \, d\textbf{x}+ \mathcal {O}(h), \end{aligned}$$

(this can be verified by using the specific form of \(u^h_c\) on each triangle). Since \(u^h_c \in H^1(\Omega )\), Theorem 7.2 gives

$$\begin{aligned} \int _\Omega \phi _t u^h_c \, d\textbf{x}\rightarrow \int _\Omega \phi _t \bar{u}^h_c \, d\textbf{x}, \end{aligned}$$

in \(H^1(\Omega )\). Thus, \(\overline{\nabla u} = \nabla \bar{u}\) in (40).

Theorem 7.6

Equation (40) holds for all \(\phi \in H^1(0,T;H^1_{\partial \Omega _0^D}(\Omega ))\).

Proof

Since the space \(H^1(0,T;C^\infty _{\partial \Omega ^D_0} (\bar{\Omega }))\) is dense in \(H^1(0,T;H^1_{\partial \Omega _0^D}(\Omega ))\) (see [1]), the equality (40) holds for all \(\phi \in H^1(0,T;H^1_{\partial \Omega _0^D} (\Omega ))\). \(\square \)

Hence, \(\bar{u}\) is a weak solution of the problem (3).

8 Strong Convergence to a Weak Solution

Next, we prove strong convergence to the weak solution.

Definition 8.1

(Strong convergence, [10]) A sequence \(\{u_n\}_{n=1}^\infty \subset X\) is said to converge to \(u \in X\), i.e., \(u_n \rightarrow u\), if \(\lim _{n \rightarrow \infty } \Vert u_n - u\Vert _X = 0\).

We also need the following definition.

Definition 8.2

[25, Definition 6.76] We say that a domain, \(\Omega \subset \mathbb {R}^d\), has the k-extension property if there exists a bounded linear mapping \(E: H^k(\Omega ) \rightarrow H^k(\mathbb {R}^d)\) such that \(Eu|_\Omega = u\) for every \(u \in H^k({\Omega }).\)

As we have assumed the spatial domain to be Lipschitz, the following result applies.

Theorem 8.3

(see e.g. [1] or [27]) Any Lipschitz domain has the k-extension property.

For a bounded domain, \(\Omega \), with the k-extension property, we have that \(H^1(\Omega )\) is compactly embedded in \(L^2(\Omega )\) (see e.g. [25]), which in turn is continuously embedded in \(H^{-1}(\Omega )\). To prove strong convergence, we need the Aubin–Lions Lemma:

Lemma 8.4

(Aubin–Lions, see e.g. [26]) Let X, B and Y be Banach spaces such that \(X \subset B \subset Y\), where the embedding, \(X \subset B\) is compact and \(B \subset Y\) is continuous. Let \(U = \{u \in L^p(0,T;X) \, | \, u_t \in L^q(0,T; Y)\}\), \(1 \le p,q < \infty \). Then U is compactly embedded in \(L^p(0,T;B)\).

A Banach space X is compactly embedded in another Banach space Y, if the following two conditions hold (see [10]):

  1. (i)

    \(\Vert u\Vert _Y \le \mathcal {C} \Vert u\Vert _X\), \((u \in X)\), for some constant \(\mathcal {C}\).

  2. (ii)

    each bounded sequence in X is precompact in Y, i.e., for a bounded sequence \(\{u_n\}_{n=1}^\infty \), there exists a subsequence, \(\{u_{n_i}\}_{n_i=1}^\infty \subseteq \{u_n\}_{n=1}^\infty \) that converges to a \(\bar{u}\) in Y.

Herein, we use \(X = H^1(\Omega )\), \(B = L^2(\Omega )\) and \(Y = H^{-1}(\Omega )\) in Lemma 8.4. Thus, since we have \(u^h_c \in L^2(0,T,H^1(\Omega ))\), it suffices to show that \((u^h_c)_t \in L^1(0,T;H^{-1}(\Omega ))\) to establish the strong convergence. That is, we need to show that the norm (see e.g. [25])

$$\begin{aligned} \Vert (u^h_c)_t\Vert _{L^1(0,T;H^{-1}(\Omega ))} = \int _0^T \sup _{\begin{array}{c} \phi \in H^1_0(\Omega ), \\ \Vert \phi \Vert _{H^1_0(\Omega )}=1 \end{array}} \int _\Omega (u^h_c)_t \phi \,d\textbf{x}dt, \end{aligned}$$
(41)

is bounded. To this end, we test the scheme (21) with a function \(\phi \in C_0^\infty (\bar{\Omega })\).

$$\begin{aligned} \int _\Omega \phi u^h_t \, dx&= \int _\Omega \phi (L_h u^h) \, d\textbf{x}= \sum _{i \in \bar{\Omega }_h} \int _{V_i} \phi |_{V_i} (L_h u^h) |_{V_i} \, d\textbf{x}. \end{aligned}$$

Note the resemblance to (29) (the only difference being the function \(\phi \) that is now vanishing on the whole boundary \(\partial \Omega \)). From derivations analogous to (30)–(35), we can recast the above equation to

$$\begin{aligned} \int _\Omega (\phi - hp) u^h_t \, d\textbf{x}&= \frac{1}{2} \sum _{i \in \bar{\Omega }_h^V} \phi _i \sum _{n \in N_i} \mu \nabla _h u^{n} \cdot \hat{\varvec{n}}_i^{n} + \frac{1}{2} \sum _{i \in \partial \Omega ^N_h} \phi _i \sum _{e \in E_i} g^{\scriptscriptstyle N}_i |\hat{\varvec{b}}(e)| \\ {}&\quad + \sum _{i \in \partial \Omega _h^R} \phi _i \sum _{e \in E_i} (g^{\scriptscriptstyle R}_i - \alpha u_i) |\hat{\varvec{b}}(e)| + \sum _{i \in \bar{\Omega }_h} \phi _i V_i F_i. \end{aligned}$$

By using \(\phi |_{\partial \Omega } = 0\) and the SBP property (see Theorem. 5.1) we have

$$\begin{aligned} \int _\Omega (\phi - hp) u^h_t \, d\textbf{x}&= - \sum _{n \in \bar{\Omega }_h^K} \int _{K_n} \nabla _h \phi ^{n} \cdot \mu \nabla _h u^{n} \, d\textbf{x}+ \sum _{i \in \bar{\Omega }_h^V} \phi _i V_i F_i \nonumber \\ {}&= - \int _\Omega \nabla _h \phi ^h \cdot \nabla _h u^h \, d\textbf{x}+ \int _\Omega \phi ^h F^h \, d\textbf{x}. \end{aligned}$$
(42)

Remark 8.5

Here, \(\int _{\Omega } \phi ^h F^h \, d\textbf{x}\) takes the same form as in Remark 7.1, except for the boundary term \(\int _{\partial \Omega ^R_h} \alpha \phi ^h w^h \, ds\) which is zero in (42) since \(\phi \) is vanishing on the entire boundary \(\partial \Omega \) in this case.

Inserting \(\phi = \phi ^h + hp\) and \(\nabla _h \phi = \nabla \phi + \mathcal {O}(h)\), we obtain

$$\begin{aligned} \int _\Omega (\phi - hp) u^h_t \, d\textbf{x}&= - \int _\Omega \left( \nabla \phi \cdot \mu \nabla _h u^h + \mathcal {O}(h) \cdot \mu \nabla _h u^h \right) \, d\textbf{x}+ \int _\Omega \left( \phi F^h - hp F^h \right) \, d\textbf{x}. \end{aligned}$$

Since \(\nabla _h u^h \in L^2(0,T;L^2_K(\bar{\Omega }_h))\) and all terms of \(F^h\) are properly bounded (see the assumptions in Sect. 2), letting \(h \rightarrow 0\) yields

$$\begin{aligned} \int _\Omega \phi u^h_t \, d\textbf{x}&= - \int _\Omega \nabla \phi \cdot \mu \overline{\nabla u} \, d\textbf{x}+ \int _\Omega \phi \bar{F} \, d\textbf{x}, \end{aligned}$$

as \(\lim _{h \rightarrow 0} (\phi - hp) = \phi \). By inserting the specific form of \(\int _\Omega \phi \bar{F} \, d\textbf{x}\) and using the Cauchy–Schwarz inequality, we obtain

$$\begin{aligned} \int _\Omega \phi u_t^h \, dx&\le \frac{1}{2} \left( \Vert \nabla \phi \Vert ^2_{L^2(\Omega )} + \mu \Vert \overline{\nabla u}\Vert ^2_{L^2(\Omega )} + \Vert \nabla \phi \Vert ^2_{L^2(\Omega )} + \Vert \overline{\nabla w}\Vert ^2_{L^2(\Omega )}\right. \nonumber \\ {}&\quad \left. + \Vert \phi \Vert ^2_{L^2(\Omega )} + \Vert \bar{w}_t\Vert ^2_{L^2(\Omega )}\right) . \end{aligned}$$
(43)

This holds for all \(\phi \in C^\infty _0(\bar{\Omega })\), and by density, it follows that the inequality holds for all \(\phi \in H^1_0(\Omega )\). Integration in time finally yields

$$\begin{aligned}&\int _0^T \sup _{\begin{array}{c} \phi \in H^1_0(\Omega ), \\ \Vert \phi \Vert _{H^1_0(\Omega )}=1 \end{array}} \int _\Omega \phi u_t^h \, d\textbf{x}dt \\&\le \int _0^T \sup _{\begin{array}{c} \phi \in H^1_0(\Omega ), \\ \Vert \phi \Vert _{H^1_0(\Omega )}=1 \end{array}} \frac{1}{2}\left( \Vert \nabla \phi \Vert ^2_{L^2(\Omega )} + \mu \Vert \overline{\nabla u}\Vert ^2_{L^2(\Omega )} \right. \nonumber \\ {}&\quad \left. + \Vert \nabla \phi \Vert ^2_{L^2(\Omega )} + \Vert \overline{\nabla w}\Vert ^2_{L^2(\Omega )} + \Vert \phi \Vert ^2_{L^2(\Omega )} + \Vert \bar{w}_t\Vert ^2_{L^2(\Omega )}\right) \, dt. \end{aligned}$$

Hence \(u_t^h \in L^1(0,T;H^{-1}(\Omega ))\), and since \((u^h_c)_t\) is \(u^h_t\) extended to the entire domain using a linear interpolant on the triangles, we also have \((u^h_c)_t \in L^1(0,T;H^{-1}(\Omega ))\). Thus, by Aubin–Lions’ lemma 8.4, the family of functions, \(U = \{u^h_c \in L^2(0,T; H^1(\Omega )) \, | \, (u^h_c)_t \in L^1(0,T; H^{-1}(\Omega ))\}\), is compactly embedded in \(L^2(0,T; L^2(\Omega ))\), meaning that \(u^h_c\) converges strongly to the weak solution.

9 Uniqueness of the Weak Solution

Assume that there are two weak solutions u, v to the problem (1) satisfying the boundary and initial data. Then \(w = u-v\) is also a weak solution with homogenous data (\(F = g^{\scriptscriptstyle D} = g^{\scriptscriptstyle N} = g^{\scriptscriptstyle R} \equiv 0\)). Take \(\phi = w\) in (10) to obtain

$$\begin{aligned} \int _\Omega ww_t \, d\textbf{x}&= \int _\Omega w (\nabla \cdot \mu \nabla w) \, d\textbf{x}. \end{aligned}$$

Integrating the right-hand side by parts, and using the fact that the boundary data is zero, we obtain

$$\begin{aligned} \frac{1}{2} \frac{d}{dt} \Vert w\Vert ^2_{L^2(\Omega )}&= -\mu \Vert \nabla w\Vert ^2_{L^2(\Omega )} - \alpha \Vert w\Vert ^2_{L^2(\partial \Omega ^R)} \le 0, \\ \Vert w(\cdot , \cdot , T)\Vert ^2_{L^2(\Omega )}&\le \Vert w(\cdot , \cdot , 0)\Vert ^2_{L^2(\Omega )} \equiv 0. \end{aligned}$$

Hence, \(\Vert w\Vert _{L^2(0,T;L^2(\Omega ))} = \Vert u - v\Vert _{L^2(0,T;L^2(\Omega ))} = 0\) and thus the weak solution is unique in \(L^2(0,T; L^2(\Omega ))\).

10 Numerical Simulations

We implement the scheme (1) and consider the manufactured solution used in [7]. That is, the exact solution is given by

$$\begin{aligned} u(x,y,t) = e^{-8\pi ^2t} \sin (2\pi x)\sin (2\pi y) + e^{-32\pi ^2t}\sin (4\pi x)\sin (4\pi y), \end{aligned}$$
(44)

which yields a zero forcing function. Furthermore, we let \(\mu = \alpha = 1\). We consider a square domain \(\Omega = [0,1] \times [0,1]\) containing a hole. The hole is located at \((x,y) = (0.5, 0.5)\), and has radius \(r = \frac{1}{8}\). We pose Dirichlet boundary conditions on the boundary of the hole, Neumann boundary conditions on \(y = 0\), \(y = 1\) and Robin boundary conditions on \(x= 0\), \(x = 1\). The boundary data is given by (44). \(t = 0.05\) is used as the final time. The scheme was run on grids containing 398, 1394, 5097, 19457 and 76166 nodes. A typical grid is depicted in Fig. 5a. All grids were generated using Gmsh (see [14]). The scheme was implemented using the Julia programming language (see [4]).

Figure 5b shows the convergence rate together with a reference line representing second-order convergence. We conclude that the scheme converges at approximately a rate of two.

Fig. 5
figure 5

a A typical mesh. b Convergence rate obtained for simulations using \(N = 398, 1394, 5097, 19457, 76166\) grid points

Full size image

11 Conclusion

Herein, we have considered a slightly modified local finite-volume approximation of the Laplacian operator proposed by Chandrashekar in [7] for discretising the heat equation in two spatial dimensions on general triangular grids. The equation was augmented with Dirichlet, Neumann and Robin boundary conditions. The Dirichlet boundary condition was imposed strongly by injection, while the Neumann and Robin conditions were imposed weakly. We demonstrated that this modification satisfies the SBP property proved in [7]. By using the energy method, a priori estimates for the numerical solution were derived. From these estimates, we were able to prove the weak convergence of the numerical solution to a weak solution of the heat equation. Thus, consistency, in a weak sense, of the Laplacian operator was established. Subsequently, we demonstrated that the numerical solution converges strongly to a weak solution by using Aubin–Lions’ lemma. Finally, the weak solution was shown to be unique. To the best of our knowledge, this is the first proof of convergence for a local finite-volume method for the Laplacian on general triangular grids. The theory presented here is straightforwardly applicable to three spatial dimensions, provided that the Laplacian approximation can be generalised to such domains.

A numerical simulation, which included Dirichlet, Neumann and Robin conditions was run on an unstructured triangulated grid containing a hole. By using the method of manufactured solutions, we demonstrated that the numerical solution converged with a second-order rate.