High Order Boundary Extrapolation Technique for Finite Difference Methods on Complex Domains with Cartesian Meshes

Abstract

The application of suitable numerical boundary conditions for hyperbolic conservation laws on domains with complex geometry has become a problem with certain difficulty that has been tackled in different ways according to the nature of the numerical methods and mesh type. In this paper we present a technique for the extrapolation of information from the interior of the computational domain to ghost cells designed for structured Cartesian meshes (which, as opposed to non-structured meshes, cannot be adapted to the morphology of the domain boundary). This technique is based on the application of Lagrange interpolation with a filter for the detection of discontinuities that permits a data dependent extrapolation, with higher order at smooth regions and essentially non oscillatory properties near discontinuities.

Appendices

Appendix 1: Computations of Intersections and Normals for Meshing

We have implemented a method to automatically compute the interior cells to $\varOmega$, the ghost cells and the normal (to $\partial \varOmega$) vectors associated to them. We require a parametric definition of the curve described by $\partial \varOmega$.

The knowledge of the Lipschitz constants of this curve is enough for computing the intersections of the boundary with the mesh lines. Indeed, with the notation $\alpha =(\alpha _x,\alpha _y):[a,b]\rightarrow {\mathbb {R}}^2$ for the curve described by $\partial \varOmega$, if $L_x$ and $L_y$ are upper bounds for the Lipschitz constants of the respective components of $\alpha$, then the values $\varDelta _xs:=\frac{h_x}{L_x}, \varDelta _ys:=\frac{h_y}{L_y}$ satisfy

$$\begin{aligned} |\alpha _x(s+\varDelta _xs)-\alpha _x(s)|\le & {} h_x\quad \forall s\in [a,b-\varDelta _xs]\\ |\alpha _y(s+\varDelta _ys)-\alpha _y(s)|\le & {} h_y\quad \forall s\in [a,b-\varDelta _ys]. \end{aligned}$$

If we detect that a vertical or horizontal mesh line passes through two values of the parameter, we use Newton–Raphson’s method, combined with bisection to safely approximate the intersection.

A similar strategy is used for the computations of the normals. In this case, given a ghost cell center, $(x,y)$, we start from the parameter value of the intersection that has generated that cell, $s_0$, for this should not be far from the value of the parameter that corresponds to the point $N(x, y)\in \partial \varOmega$ that determines the normal vector to the boundary passing by $(x, y)$.

Our search criterion is based on the sign of the scalar product of the tangent vector, $(\alpha '_x(s), \alpha '_y(s))$ and the vector formed by $(\alpha _x(s), \alpha _y(s))$ and $(x,y)$. At any rate, we aim to find $s$ satisfying the equation

$$\begin{aligned} f(s):=\langle (\alpha '_x(s),\,\,\alpha '_y(s)),\,\,(\alpha _x(s)-x,\,\, \alpha _y(s)-y)\rangle =0. \end{aligned}$$

(19)

To accomplish this, we start from $s=s_0$ and evaluate this expression at points at each side of $s_0$ until a sign change in $f$ is found. At this point we can use Newton-Raphson’s method, together with bisection, to safely approximate the solution of (19) to the required precision.

Appendix 2: Linear Stability Analysis for Inflow Boundary Conditions

Since a GKS stability analysis (see [10]) is out of the scope of this paper, we perform a simpler stability analysis to illustrate the strong influence of the spacing of the closest node to the boundary in the stability of the global numerical method.

For the sake of the exposition, we consider a linear advection equation

$$\begin{aligned} u_x+u_t=0 \end{aligned}$$

defined on $\varOmega =(0,1)$ and inflow condition at $x=0$, $u(0, t)=g(t)$. Let $h_x>0$, $\varDelta t>0$, $k\in (0,1]$ and define our grid points as $x_j:=(k+j)h_x$, $0\le j\le E(\frac{1}{h_x})$ and $t_n:=n\varDelta t$.

For solving numerically the equation, we use a first order upwind scheme:

$$\begin{aligned} U_j^{n+1}=U_j^n-\frac{\varDelta t}{h_x}\left( U_j^n-U_{j-1}^n\right) \end{aligned}$$

Assume we want to use linear interpolation at the boundary. This involves the inflow condition itself and the node $U_0^n$ at each time step $n$. Naming $G^n:=g(t_n)$, if we want to extrapolate at $x_{-1}=(k-1)h_x$, where no information is available, at time step $n$ using these two points, it can be shown by performing a linear extrapolation that

$$\begin{aligned} U_{-1}^n=\frac{1}{k}G^n+\frac{k-1}{k}U_0^n. \end{aligned}$$

Since we want to focus on a stability analysis around the left boundary, where inflow conditions are prescribed, the scheme at $x_0$ reads

$$\begin{aligned} U_0^{n+1}= & {} U_0^n-\frac{\varDelta t}{h_x}\left( U_0^n-U_{-1}^n\right) = U_0^n-\frac{\varDelta t}{h_x}\left( U_0^n-\frac{1}{k}G^n+\frac{1-k}{k}U_0^n\right) \\= & {} \left( 1-\frac{\varDelta t}{kh_x}\right) U_0^n+\frac{\varDelta t}{kh_x}G^n=U_0^n-\frac{\varDelta t}{kh_x}\left( U_0^n-G^n\right) . \end{aligned}$$

With this, following a standard Von Neumann stability analysis from this new scheme, a necessary condition for stability is

$$\begin{aligned} 0\le \frac{\varDelta t}{kh_x}\le 2, \end{aligned}$$

which is equivalent to

$$\begin{aligned} k\ge \frac{\varDelta t}{2h_x}=\frac{\text {CFL}}{2}. \end{aligned}$$

Appendix 3: Accuracy of Schemes Obtained from Extrapolation at Ghost Cells

Related to Sect. 2.1, let us denote

$$\begin{aligned} \varPhi (a_{-p-1},\ldots ,a_q)&= \hat{f}(a_{-p},\ldots ,a_{q})-\hat{f}(a_{-p-1},\dots ,a_{q-1})\\ U(h)&=(u(x-(p+1)h),\ldots ,u(x+qh)) \end{aligned}$$

for fixed $x$ and sufficiently smooth $u$, $\hat{f}$. Then (5) is equivalent to

$$\begin{aligned} F(h)=\varPhi (U(h))=h f(u(x))_x + {\mathcal {O}}(h^{r+1}), \end{aligned}$$

(20)

which is in turn equivalent to

$$\begin{aligned} F'(0)=f'(u(x))u'(x),\quad F^{(n)}(0)=0,\quad n=2,\ldots ,r. \end{aligned}$$

(21)

If $\hat{f}$ were a linear function, these relations would immediately yield a linear system of equations for its coefficients. In general, a more complicated formula equivalent to (21) can be obtained as follows: By induction on $n$ the following vectorial version of Faà di Bruno’s Formula [5] can be proved:

$$\begin{aligned} F^{(n)}(h)=\sum _{m=1}^{n}\sum _{p_1,\ldots ,p_m=1}^{n-m+1} \alpha _{p_1,\dots ,p_m} \sum _{j_1,\ldots ,j_m=-p-1}^{q} \frac{\partial ^{m}\varPhi (U(h))}{\partial a_{j_1}\ldots \partial a_{j_m}} U_{j_1}^{(p_1)}(h)\ldots U_{j_m}^{(p_m)}(h), \end{aligned}$$

(22)

for suitable and uniquely defined coefficients $\alpha _{p_1,\ldots ,p_m}$. From (23) and the fact that $U_j^{(\nu )}$ $=j^\nu u^{(\nu )}(x)$, for $j=-p-1,\ldots ,q$ and $\nu \ge 0$, we deduce that

$$\begin{aligned}&F^{(n)}(0)=\sum _{m=1}^{n}\sum _{p_l=1}^{n-m+1} \alpha _{p_1,\ldots ,p_m} \sum _{j_l=-p-1}^{q} \varPhi ^{(j_1,\dots ,j_m)}(x) j_1^{p_1}\ldots j_m^{p_m}u^{(p_1)}(x)\ldots u^{(p_m)}(x),\nonumber \\&\varPhi ^{(j_1,\ldots ,j_m)}(x)= \frac{\partial ^{m}\varPhi }{\partial a_{j_1}\ldots \partial a_{j_m}}\left( u(x),\ldots ,u(x)\right) . \end{aligned}$$

(23)

For $n=1$, (23) is satisfied with $\alpha _{p_1}=1$:

$$\begin{aligned} F'(0)= \Big (\sum _{j_1=-p-1}^{q} j_1 \frac{\partial \varPhi }{\partial a_{j_1}}(u(x),\ldots ,u(x)) \Big ) u'(x). \end{aligned}$$

It follows from (21) that $F'(0)=f'(u(x))u'(x)$ for any sufficiently smooth $u$ if and only if

$$\begin{aligned} \sum _{j=-p-1}^{q} j \frac{\partial \varPhi }{\partial a_{j}}(u(x),\ldots ,u(x))=f'(u(x)), \end{aligned}$$

i.e.,

$$\begin{aligned} \sum _{j=-p-1}^{q} j \frac{\partial \varPhi }{\partial a_{j}}(u,\ldots ,u)=f'(u). \end{aligned}$$

Now, (20) may be obtained by any means but, as we have just seen, they are equivalent to (21), (23), which only depends on derivatives of $u$ at $x$, as long as $u$ is smooth enough at a convex open set containing the corresponding stencil.

Assume now that that $u(x-(p+1)h)$ is replaced by, $\widetilde{u}(x-(p+1)h)$, the evaluation at $x-(p+1)h$ of a polynomial that interpolates $u$ at another stencil $x+p'h,\ldots x+q'h$, with $p'\ge -p$. Let us denote $\overline{q}=\max (q, q')$ and

$$\begin{aligned} V(h)&=(\widetilde{u}(x-(p+1)h),u(x-ph), \ldots ,u(x+qh)) = W(u(x-ph),\ldots ,u(x+\overline{q}h)),\\ G(h)&=\varPhi (V(h)). \end{aligned}$$

If $u$ is smooth in a convex open set containing $[x-(p+1)h,\ldots ,x+\overline{q}h]$, then

$$\begin{aligned} u(x-(p+1)h)-\widetilde{u}(x-(p+1)h)={\mathcal {O}}(h^s), s\ge 2 \end{aligned}$$

and, therefore,

$$\begin{aligned} F(h)-G(h)=\frac{\partial {\varPhi }(\xi _h)}{\partial a_{-p-1}}\left( u(x-(p+1)h)-\widetilde{u}(x-(p+1)h)\right) = {\mathcal {O}}(h^s). \end{aligned}$$

It then follows from (20) that

$$\begin{aligned} G(h)=h f(u(x))_x + {\mathcal {O}}(h^{r'+1}), r'=\min (r, s-1), \end{aligned}$$

(24)

under the assumption of $u$ being smooth in a convex open set containing $[x-(p+1)h,\ldots ,x+\overline{q}h]$. As argued above, this is equivalent to

$$\begin{aligned} G'(0)=f(u(x))_x,\quad G^{(n)}(0)=0,\quad n=2,\ldots ,r', \end{aligned}$$

(25)

which translates into a functional relationship akin to that for $F$. But now $G$ depends only on the stencil $x-ph,\ldots ,x+\overline{q}h$ and therefore (25) on $u$ being smooth in a convex open set containing $[x-ph,\ldots ,x+\overline{q}h]$, i.e.,

$$\begin{aligned} \frac{G(h)}{h}=f(u(x))_x+{\mathcal {O}}(h^{r'}) \end{aligned}$$

if $u$ is smooth in a convex open set containing $[x-ph,\ldots ,x+\overline{q}h]$.

Therefore, we have concluded that $r$th order extrapolation applied to $r$th order schemes decreases the order of the scheme near the boundaries to $r-1$. It can be seen that this order is sharp. This is in apparent contradiction with the results obtained in Sect. 5.1, where fifth order extrapolation applied together with fifth order schemes yields fifth order global errors, both in 1-norm and $\infty$-norm.

We next justify that the global errors are of order $r$ in 1-norm if the local error is of order $r+1$ at all but a bounded number (with respect to $M$) of nodes where the order is a unit less.

Let us denote $v^n_j=u(x_j, t_n)$, $j=1,\ldots ,M$, for $M$ the total number of spatial nodes, and ${\mathcal {H}}={\mathcal {H}}_{\lambda ,h}$, $\lambda =\varDelta t/h$, the operator of the numerical scheme:

$$\begin{aligned} u^{n+1}={\mathcal {H}}u^n, \quad {\mathcal {H}}:{\mathbb {R}}^M\rightarrow {\mathbb {R}}^M. \end{aligned}$$

(26)

With this notation, the global error is $e^n=v^n-u^n$ and the local error can be written as $f_j^n=v_{j}^{n+1}-{\mathcal {H}}v^n_j,$ which yields:

$$\begin{aligned} v^{n+1}={\mathcal {H}}(v^n)+f^n. \end{aligned}$$

(27)

Subtracting (26) and (27) and using the mean value theorem for ${\mathcal {H}}$ we obtain a relationship between local and global error:

$$\begin{aligned} e^{n+1}= B_ne^n+f^n, \end{aligned}$$

(28)

where

$$\begin{aligned} B_n=\int _0^1{\mathcal {H}}'(u^n+s(v^n-u^n))ds, \end{aligned}$$

which, by induction on $n$, yields

$$\begin{aligned} e^n=\sum _{k=0}^{n-1}\left( \prod _{m=k+1}^{n-1}B_m\right) f^k. \end{aligned}$$

(29)

Now, we recall that the classical argument to obtain convergence from consistency and stability proceeds as follows: Given some fixed $T>0$, if for $h \in (0, h_0)$,

$$\begin{aligned} \Vert f^k\Vert \le Ah^{r+1}, \quad \forall k \text { such that } k\varDelta t \le T, \end{aligned}$$

(30)

and, for the induced matrix norm, $\Vert B_m\Vert \le 1+D h$, then it easily follows that $\Vert e^n\Vert \le C h^{r}$ for all $h\in (0, h_0)$ and all $n$ such that $n\varDelta t=n\lambda h \le T$.

This analysis is sufficient for the 1-norm if a bounded number of nodes $x_j$ satisfy $f_j^k={\mathcal {O}}(h^r)$, whereas the rest satisfy $f_j^k={\mathcal {O}}(h^{r+1})$, since then $\Vert f_j \Vert _1=h\sum _{j=1}^{M} |f_j^k| = {\mathcal {O}}(h^{r+1})$. But this analysis does not explain that $\Vert e^{n}\Vert _{\infty } = {\mathcal {O}}(h^r)$, since in this case we have $\Vert f^k\Vert _{\infty } ={\mathcal {O}}(h^{r})$, i.e., (30) does not hold, so this argument would yield $\Vert e^{n}\Vert _{\infty } = {\mathcal {O}}(h^{r-1})$. A finer analysis, based on the examination of the entries of the vectors involved in (29), is then needed to obtain $\Vert e^{n}\Vert _{\infty } = {\mathcal {O}}(h^r)$. The analysis of

$$\begin{aligned} e_p^n=\sum _{k=0}^{n-1}\sum _{q=1}^{M}\Big (\prod _{m=k+1}^{n-1}B_m\Big )_{pq}f_q^k, \end{aligned}$$

(31)

can be quite involved in general and it is out of the scope of this work. We will show a particular case in which we obtain global errors of order 1 (in $\infty$-norm) from local errors with some bounded number of entries of order 1 and the rest of order 2. We consider the upwind scheme

$$\begin{aligned} u_{i}^{n+1}=(1-\lambda ) u_{i}^{n} + \lambda u_{i-1}^{n},\quad i=1,\ldots , M \end{aligned}$$

for the linear advection equation $u_t+u_x=0$ with time dependent boundary condition at $x=0$, $u(0, t)=g(t)$. The value of $u_{0}^{n}\approx u(x_0, t_n)$ needs to be supplied according to this boundary condition. Prescribing $u_{0}^{n}=g(t_n)=u(0, t_n)$ gives a first order approximation of $u(x_0, t_n)=u(-h/2, t_n)$ and effectively yields a local error of order $r=1$ at $x_1$:

$$\begin{aligned} |f_1^{n}|=|u(h/2, t_{n}+\varDelta t)-(1-\varDelta t/h)u(h/2, t_n)-\varDelta t/h g(t_n)| \le A_1h, \end{aligned}$$

whereas the rest of the entries of the local error satisfy

$$\begin{aligned} |f_j^{n}| \le Ah^2, j=2,\ldots , M, \end{aligned}$$

i.e., are of order $r+1=2$, for suitable positive constants $A_1, A$.

For this linear case, all the matrices $B_m$ are $(1-\lambda )I_{M}+\lambda N_M$, where $N=N_M$ is the nilpotent matrix whose powers are given by $(N^j)_{pq}=\delta _{p-q,j}$, $j\ge 1$, $1\le p, q\le M$. We compute now, with the change of variables $k'=n-k-1 \rightarrow k$:

$$\begin{aligned} e_p^n&=\sum _{k=0}^{n-1}\sum _{q}^{M}((1-\lambda )I+\lambda N)^{n-k-1}_{pq}f_q^k\\&=\sum _{k=0}^{n-1}\sum _{q=1}^{M}\sum _{j=0}^{n-k-1} \left( {\begin{array}{c}n-k-1\\ j\end{array}}\right) (1-\lambda )^{n-k-1-j} \lambda ^{j}N^{j}_{pq}f_q^k\\&=\sum _{k=0}^{n-1}\sum _{q=1}^{M}\sum _{j=0}^{n-k-1} \left( {\begin{array}{c}n-k-1\\ j\end{array}}\right) (1-\lambda )^{n-k-1-j} \lambda ^{j}\delta _{p-q,j}f_q^k\\&=\sum _{k=0}^{n-1}\sum _{q=\max (1,p-n+k+1)}^{p}\left( {\begin{array}{c}n-k-1\\ p-q\end{array}}\right) (1-\lambda )^{n-k-1-p+q} \lambda ^{p-q}f_q^k\\&=\sum _{k=0}^{n-1}\sum _{q=\max (1,p-k)}^{p}\left( {\begin{array}{c}k\\ p-q\end{array}}\right) (1-\lambda )^{k-p+q} \lambda ^{p-q}f_q^{n-k-1}\\&=\sum _{q=1}^{p} \left( \lambda ^{p-q} \sum _{k=p-q}^{n-1}\left( {\begin{array}{c}k\\ p-q\end{array}}\right) (1-\lambda )^{k-p+q}\right) f_q^{n-k-1}. \end{aligned}$$

We can bound

$$\begin{aligned} \lambda ^{s} \sum _{k=s}^{n-1}\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s} \le \lambda ^{s} \sum _{k=s}^{\infty }\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s} \end{aligned}$$

for $s=p-q$ and use that $\lambda ^{s} \sum _{k=s}^{n-1}\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s} \le \lambda ^{s} \sum _{k=s}^{\infty }\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s}=\frac{1}{\lambda }$ (see next lemma) to conclude:

$$\begin{aligned} |e_p^n|\le \frac{1}{\lambda }\sum _{q=1}^{p} |f_q^{n-k-1}|\le \frac{1}{\lambda }\left( A_1 h+(p-1) Ah^{2}\right) \le \frac{1}{\lambda }\left( A_1 h+M Ah^{2} \right) =\frac{A_1+A}{\lambda }h, \end{aligned}$$

for any $p=1,\ldots ,M$, i.e., $\Vert e^n\Vert _{\infty }={\mathcal {O}}(h^r)$, $r=1$, as proposed.

Lemma 1

For any $\lambda \in (0,1)$ and $s \in \mathbb {N}$,

$$\begin{aligned} \sum _{k=s}^{\infty }\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s}= \lambda ^{-s-1}. \end{aligned}$$

Proof

Set $A_s=\sum _{k=s}^{\infty }\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s}$. From the identity $\left( {\begin{array}{c}k\\ s\end{array}}\right) =\left( {\begin{array}{c}k-1\\ s\end{array}}\right) +\left( {\begin{array}{c}k-1\\ s-1\end{array}}\right)$, for $0\le s < k$, and the change of variables $k'=k-1$ we get:

$$\begin{aligned} \sum _{k=s}^{\infty }\left( {\begin{array}{c}k\\ s\end{array}}\right) (1-\lambda )^{k-s}&= 1+\sum _{k=s+1}^{\infty }\left( \left( {\begin{array}{c}k-1\\ s\end{array}}\right) +\left( {\begin{array}{c}k-1\\ s-1\end{array}}\right) \right) (1-\lambda )^{k-s} \\&= 1+\sum _{k=s+1}^{\infty }\left( {\begin{array}{c}k-1\\ s\end{array}}\right) (1-\lambda )^{k-s}+\sum _{k=s+1}^{\infty }\left( {\begin{array}{c}k-1\\ s-1\end{array}}\right) (1-\lambda )^{k-s}\\&= 1+(1-\lambda )\sum _{k'=s}^{\infty }\left( {\begin{array}{c}k'\\ s\end{array}}\right) (1-\lambda )^{k'-s}+\sum _{k'=s}^{\infty }\left( {\begin{array}{c}k'\\ s-1\end{array}}\right) (1-\lambda )^{k'-s} \\&=(1-\lambda )\sum _{k'=s}^{\infty }\left( {\begin{array}{c}k'\\ s\end{array}}\right) (1-\lambda )^{k'-s}+\sum _{k'=s-1}^{\infty }\left( {\begin{array}{c}k'\\ s-1\end{array}}\right) (1-\lambda )^{k'-s}. \end{aligned}$$

We obtain that $A_{s}=(1-\lambda )A_s+A_{s-1}$ i.e., $A_{s}=\lambda ^{-1}A_{s-1}$, which immediately yields

$$\begin{aligned} A_{s}=\lambda ^{-s}A_0=\lambda ^{-s}\sum _{k=0}^{\infty }(1-\lambda )^k=\lambda ^{-s-1}. \end{aligned}$$

$\square$