Fourth-order conservative non-splitting semi-Lagrangian Hermite WENO schemes for kinetic and fluid simulations

For convenience, we require $\Delta x_{i}\equiv\Delta x,\quad\Delta y_{j}\equiv\Delta y\quad\forall i,j$. Based on the $P^{1}$ DG solution, $u^{n}(x,y)$, we recover a piecewise $P^{3}$ polynomial,

$$H^{n}(x,y)=H^{(i,j)}(x,y),~{}~{}(x,y)\in I_{i,j},~{}~{}\forall(i,j),$$

(2.7)

where $H^{(i,j)}(x,y)\in P^{3}(I_{i,j})$. We define a set of local orthogonal basis of the high order polynomial denoted as $\{P_{l}(x,y)\}$ with:

$$\begin{split}&P_{1}(x,y)=1,~{}~{}P_{2}(x,y)=\mu_{i}(x):=\frac{x-x_{i}}{\Delta x},~{}~{}P_{3}(x,y)=\nu_{j}(y):=\frac{y-y_{j}}{\Delta y},\\ &P_{4}(x,y)=\mu_{i}^{2}-\frac{1}{12},~{}~{}P_{5}(x,y)=\mu_{i}\nu_{j},~{}~{}P_{6}(x,y)=\nu_{j}^{2}-\frac{1}{12},\\ &P_{7}(x,y)=\mu_{i}^{3}-\frac{3}{20}\mu_{i},~{}~{}P_{8}(x,y)=\left(\mu_{i}^{2}-\frac{1}{12}\right)\nu_{j},~{}~{}P_{9}(x,y)=\mu_{i}\left(\nu_{j}^{2}-\frac{1}{12}\right),~{}~{}\\ &P_{10}(x,y)=\nu_{j}^{3}-\frac{3}{20}\nu_{j},~{}~{}P_{11}(x,y)=\left(\mu_{i}^{2}-\frac{1}{12}\right)\left(\nu_{j}^{2}-\frac{1}{12}\right).\end{split}$$

(2.8)

We also define that $\overline{u}_{5}^{n}:=\overline{u}_{i,j}^{n}$, $\overline{v}_{5}^{n}:=\overline{v}_{i,j}^{n}$, $\overline{w}_{5}^{n}:=\overline{w}_{i,j}^{n}$, $I_{5}:=I_{i,j}$ and other $\{\overline{u}_{s}^{n}\}$, $\{\overline{v}_{s}^{n}\}$, $\{\overline{w}_{s}^{n}\}$, $\{I_{s}\}$ represent corresponding moments and Eulerian cells based on the serial numbers in Figure 2.2. Then, the 2-D HWENO reconstruction methods over $I_{i,j}$ are summarized as follows.

Step 1. Reconstruct the first-order moments.

The first-order moments, $\{\overline{v}^{n}_{5},\overline{w}^{n}_{5}\}$, can be extremely large when $u(x,y,t^{n})$ is discontinuous in $I_{i,j}$ since $\overline{v}^{n}_{5}\sim\frac{\Delta x}{12}\frac{\partial u}{\partial x}|_{(x_{i},y_{j})}$ and $\overline{w}^{n}_{5}\sim\frac{\Delta y}{12}\frac{\partial u}{\partial y}|_{(x_{i},y_{j})}$. Hence, before recovering a cubic polynomial, we reconstruct the first-order moments with the following two goals: it will provide high-order approximations of first-order moments when $u(x,y,t^{n})$ is smooth in $I_{i,j}$; when $u(x,y,t^{n})$ is discontinuous in $I_{i,j}$, first-order moments will be reduced to a reasonable level.

The two first-order moments can be regarded as local indicators of the changing rates for the $x$- and $y$-dimensions. Each of them is highly independent of the other dimension. Hence, the reconstruction is performed in a dimension-by-dimension manner. Below, we take the $x$ direction as an example to illustrate the procedure for reconstructing its first-order moment.

Step 1.1. Compute approximations to the first-order moment and smoothness indicators from 1-D polynomial reconstructions.

1. 1.

   Construct a quartic polynomial $p_{0}(x)$, and three quadratic polynomials $\{p_{k}(x)\}_{k=1}^{3}$ satisfying

   $$\begin{split}&\frac{1}{\Delta x}\int_{I^{x}_{i+l}}p_{0}(x)dx=\overline{u}_{i+l,j}^{n},~{}~{}~{}l=-1,0,1,\\ &\frac{1}{\Delta x}\int_{I^{x}_{i+l}}p_{0}(x)\left(\frac{x-x_{i+l}}{\Delta x}\right)dx=\overline{v}_{i+l,j}^{n},~{}~{}~{}l=-1,1,\\ \end{split}$$

   (2.9)

   and

   $$\begin{split}&\frac{1}{\Delta x}\int_{I^{x}_{i+l}}p_{1}(x)dx=\overline{u}_{i+l,j}^{n}~{}~{}l=-1,0,~{}~{}~{}\frac{1}{\Delta x}\int_{I^{x}_{i-1}}p_{1}(x)\left(\frac{x-x_{i-1}}{\Delta x}\right)dx=\overline{v}_{i-1,j}^{n};\\ &\frac{1}{\Delta x}\int_{I^{x}_{i+l}}p_{2}(x)dx=\overline{u}_{i+l,j}^{n}~{}~{}l=-1,0,1;\\ &\frac{1}{\Delta x}\int_{I^{x}_{i+l}}p_{3}(x)dx=\overline{u}_{i+l,j}^{n}~{}~{}l=0,1,~{}~{}~{}\frac{1}{\Delta x}\int_{I^{x}_{i+1}}p_{1}(x)\left(\frac{x-x_{i+1}}{\Delta x}\right)dx=\overline{v}_{i+1,j}^{n}.\\ \end{split}$$

   (2.10)

2. 2.

   Compute the first-order moments of $\{p_{k}(x)\}_{k=0}^{3}$:

   $$\begin{split}&\widetilde{\overline{v}}_{i,j,0}^{n}:=\frac{1}{\Delta x}\int_{I^{x}_{i}}p_{0}(x)\left(\frac{x-x_{i}}{\Delta x}\right)dx=-\frac{5}{76}\overline{u}^{n}_{i-1,j}-\frac{11}{38}\overline{v}^{n}_{i-1,j}-\frac{11}{38}\overline{v}^{n}_{i+1,j}+\frac{5}{76}\overline{u}^{n}_{i+1,j},\\ &\widetilde{\overline{v}}_{i,j,1}^{n}:=\frac{1}{\Delta x}\int_{I^{x}_{i}}p_{1}(x)\left(\frac{x-x_{i}}{\Delta x}\right)dx=\frac{1}{6}\overline{u}^{n}_{i,j}-\frac{1}{6}\overline{u}^{n}_{i-1,j}-\overline{v}^{n}_{i-1,j},\\ &\widetilde{\overline{v}}_{i,j,2}^{n}:=\frac{1}{\Delta x}\int_{I^{x}_{i}}p_{2}(x)\left(\frac{x-x_{i}}{\Delta x}\right)dx=\frac{1}{24}\overline{u}^{n}_{i+1,j}-\frac{1}{24}\overline{u}^{n}_{i-1,j},\\ &\widetilde{\overline{v}}_{i,j,3}^{n}:=\frac{1}{\Delta x}\int_{I^{x}_{i}}p_{2}(x)\left(\frac{x-x_{i}}{\Delta x}\right)dx=\frac{1}{6}\overline{u}^{n}_{i+1,j}-\frac{1}{6}\overline{u}^{n}_{i,j}-\overline{v}^{n}_{i+1,j}.\\ \end{split}$$

   (2.11)

3. 3.

   Compute the smoothness indicators $\{\beta_{k}\}_{k=0}^{3}$ of $\{p_{k}(x)\}_{k=0}^{3}$ [23, 21, 33]:

   $$\beta_{k}=\sum_{l=1}^{r}\frac{1}{\Delta x}\int_{I_{i}^{x}}\left(\Delta x^{l}\frac{\partial^{l}}{\partial x^{l}}p_{k}(x)\right)^{2}dx~{}~{}~{}k=0,1,2,3$$

   (2.12)

   with $r$ representing the degree of the corresponding polynomial. Here, the smoothness indicators $\{\beta_{k}\}_{k=1}^{3}$ can be explicitly expressed by

   $$\begin{split}&\beta_{1}=\left(12\widetilde{\overline{v}}_{i,j,1}^{n}\right)^{2}+156\left(\widetilde{\overline{v}}_{i,j,1}^{n}-\overline{v}^{n}_{i-1,j}\right)^{2},\\ &\beta_{2}=\left(12\widetilde{\overline{v}}_{i,j,2}^{n}\right)^{2}+\frac{13}{12}\left(\overline{u}^{n}_{i+1,j}-2\overline{u}^{n}_{i,j}+\overline{u}^{n}_{i-1,j}\right)^{2},\\ &\beta_{3}=\left(12\widetilde{\overline{v}}_{i,j,3}^{n}\right)^{2}+156\left(\overline{v}^{n}_{i+1,j}-\widetilde{\overline{v}}_{i,j,3}^{n}\right)^{2}.\end{split}$$

   (2.13)

   We refer to [40] for the explicit expression of $\beta_{0}$.

4. 4.

   Compute a full stencil global reference smoothness indicator [20]:

   $$\tau:=\left(\frac{|\beta_{0}-\beta_{1}|+|\beta_{0}-\beta_{3}|}{2}\right)^{2}.$$

   (2.14)

   By the Taylor expansion, we can find that $\tau=O(\Delta x^{6})$, if there is no discontinuity involved.

Step 1.2 Weight the collected first-order moments.

Below, we first introduce the weighting strategy for the HWENO-1.

1. 1.

   Compute the nonlinear weights by

   $$\omega_{k}=\frac{\overline{\omega}_{k}}{\sum_{l}{\overline{\omega}_{l}}}~{}~{}~{}\text{with}~{}~{}~{}\overline{\omega}_{k}=\gamma_{k}\left(1+\frac{\tau}{\beta_{k}+\epsilon}\right)~{}~{}~{}k=0,1,3,$$

   (2.15)

   where $\epsilon=10^{-40}$ is set to avoid the denominator being zero. The linear weights, $\{\gamma_{0},\gamma_{1},\gamma_{3}\}$, are chosen as $\{0.6,0.2,0.2\}$ in this paper. We refer to [43, 40] for more details on linear and nonlinear weights.

2. 2.

   Reconstruct the $x$-dimension first-order moment, $\overline{v}_{5}^{n}$, by:

   $$\widetilde{\overline{v}}_{5}^{n}=\frac{\omega_{0}}{\gamma_{0}}\left(\widetilde{\overline{v}}_{i,j,0}^{n}-\gamma_{1}\widetilde{\overline{v}}_{i,j,1}^{n}-\gamma_{3}\widetilde{\overline{v}}_{i,j,3}^{n}\right)+\omega_{1}\widetilde{\overline{v}}_{i,j,1}^{n}+\omega_{3}\widetilde{\overline{v}}_{i,j,3}^{n}.$$

   (2.16)

The HWENO-2 reconstructs the first-order moment $\overline{v}^{n}_{5}$ as follows.

1. 1.

   Separate the discontinuities from broad-band smooth fluctuations as illustrated in [16, 15] by first taking

   $$\eta_{k}=\left(1+\frac{\tau}{\beta_{k}+\epsilon}\right)^{6}~{}~{}~{}k=1,2,3,$$

   (2.17)

   where $\epsilon=10^{-40}$ is used to avoid the denominator being zero as in [2]. If these is no discontinuity involved, we can find that $\eta_{k}\approx 1$ for all $k$. If there is discontinuity involved for the global three-cells stencil, the $\eta$ value for a smooth small stencil can be greatly enlarge with a magnitude of $O(\Delta x^{-12})$. Then, we normalize $\{\eta_{k}\}_{k=1}^{3}$ with

   $$\kappa_{k}=\frac{\eta_{k}}{\sum_{l=1}^{3}\eta_{l}}~{}~{}~{}k=1,2,3.$$

   (2.18)

2. 2.

   Reconstruct the $x$-dimension first-order moment $\overline{v}_{5}^{n}$ by:

   $$\widetilde{\overline{v}}_{5}^{n}=\widetilde{\overline{v}}_{i,j,0}^{n},~{}~{}~{}\text{if}~{}~{}~{}\min_{k}\kappa_{k}>C_{T},$$

   (2.19)

   otherwise

   $$\widetilde{\overline{v}}_{5}^{n}=\begin{cases}\begin{split}&\widetilde{\overline{v}}_{i,j,1}^{n},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if}~{}~{}~{}\kappa_{1}>\kappa_{3},\\ &\widetilde{\overline{v}}_{i,j,3}^{n},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if}~{}~{}~{}\kappa_{3}>\kappa_{1}.\end{split}\end{cases}$$

   (2.20)

   where $C_{T}$ is a parameter deciding whether a corresponding polynomial should be involved [16]. We empirically choose $C_{T}=10^{-3}$ as in [20]. For treating discontinuity, unlike [16, 20], we only choose the smoothest polynomial for reconstruction, since the weighted summation does not increase the order of accuracy for our case. When $\{\kappa_{k}\}$ matches the smoothness of the three-cells stencil, we directly use $\widetilde{\overline{v}}_{i,j,0}^{n}$ for optimal accuracy.

The first-order moment in $y$ direction can be reconstructed in a similar way; we the reconstructed first-order moment in this direction by $\widetilde{\overline{w}}_{5}^{n}$.

Step 2. Recover the reconstructed $H^{n}(x,y)$ on $I_{i,j}$.

Step 2.1. Collect information from different 2-D polynomials.

1. 1.

   Construct a polynomial $\widetilde{q}_{0}(x,y):=\sum_{l=1}^{11}a^{q_{0}}_{l}P_{l}(x,y)$ such that

   $$\begin{split}&\frac{1}{\Delta x\Delta y}\iint_{I_{s}}\widetilde{q}_{0}(x,y)dxdy=\overline{u}_{s}^{n},~{}~{}~{}s=1,2,\ldots,9,\\ &\frac{1}{\Delta x\Delta y}\iint_{I_{5}}\widetilde{q}_{0}(x,y)\left(\frac{x-x_{i}}{\Delta x}\right)dxdy=\widetilde{\overline{v}}_{5}^{n},\\ &\frac{1}{\Delta x\Delta y}\iint_{I_{5}}\widetilde{q}_{0}(x,y)\left(\frac{y-y_{j}}{\Delta y}\right)dxdy=\widetilde{\overline{w}}_{5}^{n}.\\ \end{split}$$

   (2.21)

   Let $q_{0}(x,y)=\sum_{l=1}^{10}a^{q_{0}}_{l}P_{l}(x,y)$, which is the orthogonal projection of $\widetilde{q}_{0}(x,y)$ to $P^{3}(I_{i,j})$. We provide the explicit expressions of $\{a^{q_{0}}_{l}\}_{l=1}^{10}$ in Appendix A.

2. 2.

   Construct four quadratic polynomial $\{q_{k}(x,y)\}_{k=1}^{4}:=\{\sum_{l=1}^{6}a^{q_{k}}_{l}P_{l}(x,y)\}_{k=1}^{4}$ satisfying

   $$\begin{split}&\frac{1}{\Delta x\Delta y}\iint_{I_{s}}q_{k}(x,y)dxdy=\overline{u}^{n}_{s},\\ &\frac{1}{\Delta x\Delta y}\iint_{I_{5}}q_{k}(x,y)\frac{x-x_{i}}{\Delta x}dxdy=\widetilde{\overline{v}}^{n}_{5},\\ &\frac{1}{\Delta x\Delta y}\iint_{I_{5}}q_{k}(x,y)\frac{y-y_{i}}{\Delta y}dxdy=\widetilde{\overline{w}}^{n}_{5},\end{split}$$

   (2.22)

   where

   $$\begin{split}&s=1,2,4,5,~{}~{}\text{for}~{}~{}k=1;~{}~{}s=2,3,5,6,~{}~{}\text{for}~{}~{}k=2;\\ &s=4,5,7,8,~{}~{}\text{for}~{}~{}k=3;~{}~{}s=5,6,8,9,~{}~{}\text{for}~{}~{}k=4.\end{split}$$

   (2.23)

   The explicit expressions of $\{a^{q_{k}}_{l}\}$ are given in Appendix A.

3. 3.

   Compute the smoothness indicator [23, 21, 33] $\{\hat{\beta}_{k}\}_{k=0}^{4}$ of $\{q_{k}(x,y)\}_{k=0}^{4}$:

   $$\begin{split}\hat{\beta}_{0}=&\frac{1}{\Delta x\Delta y}\sum_{l_{1}+l_{2}<=3}\iint_{I_{i,j}}\left(\Delta x^{l_{1}}\Delta y^{l_{2}}\frac{\partial^{|l_{1}+l_{2}|}}{\partial x^{l_{1}}\partial y^{l_{2}}}q_{0}(x,y)\right)^{2}dxdy;\\ \hat{\beta}_{k}=&\frac{1}{\Delta x\Delta y}\sum_{l_{1}+l_{2}<=2}\iint_{I_{i,j}}\left(\Delta x^{l_{1}}\Delta y^{l_{2}}\frac{\partial^{|l_{1}+l_{2}|}}{\partial x^{l_{1}}\partial y^{l_{2}}}q_{k}(x,y)\right)^{2}dxdy,~{}~{}k=1,2,3,4.\end{split}$$

   (2.24)

   The explicit expression of $\{\hat{\beta}_{k}\}_{k=0}^{4}$ can be given by

   $$\begin{split}\hat{\beta}_{0}=&\left(a^{q_{0}}_{2}+\frac{1}{10}a^{q_{0}}_{7}\right)^{2}+\left(a^{q_{0}}_{3}+\frac{1}{10}a^{q_{0}}_{10}\right)^{2}+\frac{13}{3}\left(a^{q_{0}}_{4}\right)^{2}+\frac{7}{6}\left(a^{q_{0}}_{5}\right)^{2}+\frac{13}{3}\left(a^{q_{0}}_{6}\right)^{2}\\ &+\frac{781}{20}\left(a^{q_{0}}_{7}\right)^{2}+\frac{47}{10}\left(a^{q_{0}}_{8}\right)^{2}+\frac{47}{10}\left(a^{q_{0}}_{9}\right)^{2}+\frac{781}{20}\left(a^{q_{0}}_{10}\right)^{2};\\ \hat{\beta}_{k}=&\left(a^{q_{k}}_{2}\right)^{2}+\left(a^{q_{k}}_{3}\right)^{2}+\frac{13}{3}\left(a^{q_{k}}_{4}\right)^{2}+\frac{7}{6}\left(a^{q_{k}}_{5}\right)^{2}+\frac{13}{3}\left(a^{q_{k}}_{6}\right)^{2}~{}~{}~{}\text{for}~{}~{}k=1,2,3,4.\end{split}$$

   (2.25)

4. 4.

   Compute a full stencil global reference smoothness indicator:

   $$\hat{\tau}=\left(\frac{|\hat{\beta}_{0}-\hat{\beta_{1}}|+|\hat{\beta}_{0}-\hat{\beta_{2}}|+|\hat{\beta}_{0}-\hat{\beta_{3}}|+|\hat{\beta}_{0}-\hat{\beta_{4}}|}{4}\right)^{2}.$$

   (2.26)

   Similarly, by the Taylor expansion, we can easily check that $\hat{\tau}=O(\Delta x^{6})$, if there is no discontinuity involved.

Step 2.2 Weight the collected 2-D polynomials.

The HWENO-1 weights the polynomials as follows.

1. 1.

   Compute the nonlinear weights by

   $$\hat{\omega}_{k}=\frac{\widetilde{\omega}_{k}}{\sum_{l}{\widetilde{\omega}_{l}}}~{}~{}~{}\text{with}~{}~{}~{}\widetilde{\omega}_{k}=\gamma_{k}\left(1+\frac{\tau}{\beta_{k}+\epsilon}\right)~{}~{}~{}k=0,1,2,3,4,$$

   (2.27)

   where $\epsilon=10^{-40}$ is set to avoid the denominator being zero. The linear weights, $\{\gamma_{k}\}_{k=0}^{4}$ are chosen as $\{0.6,0.1,0.1,0.1,0.1\}$ in this paper.

2. 2.

   Construct $H^{(i,j)}(x,y):=\sum_{l=1}^{10}a_{l}P_{l}(x,y)$:

   $$\begin{split}H^{(i,j)}(x,y)=\frac{\hat{\omega}_{0}}{\hat{\gamma}_{0}}\left(q_{0}(x)-\sum_{k=1}^{4}\hat{\gamma}_{k}q_{k}(x,y)\right)+\sum_{k=1}^{4}\hat{\omega}_{k}q_{k}(x,y),~{}~{}~{}(x,y)\in I_{i,j}.\end{split}$$

   (2.28)

   Here, the coefficients $\{a_{l}\}$ can be explicitly given by

   $$\begin{split}&a_{1}=\overline{u}^{n}_{5};~{}~{}a_{2}=12\widetilde{\overline{v}}^{n}_{5};~{}~{}a_{3}=12\widetilde{\overline{w}}^{n}_{5};\\ &a_{l}=\frac{\hat{\omega}_{0}}{\hat{\gamma}_{0}}a^{q_{0}}_{l}+\sum_{k=1}^{4}\left(\hat{\omega}_{k}-\frac{\hat{\omega}_{0}}{\hat{\gamma}_{0}}\hat{\gamma}_{k}\right)a^{q_{k}}_{l}~{}~{}\text{for}~{}~{}l=4,5,6;\\ &a_{l}=\frac{\hat{\omega}_{0}}{\hat{\gamma}_{0}}a^{q_{0}}_{l}~{}~{}\text{for}~{}~{}l=7,8,\cdots,10.\end{split}$$

   (2.29)

The HWENO-2 weights the polynomials as follows.

1. 1.

   Compute the separation parameters $\{\hat{\eta}_{k}\}_{k=1}^{4}$ and the corresponding normalized parameters $\{\hat{\kappa}_{k}\}_{k=1}^{4}$ similarly:

   $$\hat{\kappa}_{k}=\frac{\hat{\eta}_{k}}{\sum_{l=1}^{3}\hat{\eta}_{l}}~{}~{}~{}\text{with $\hat{\eta}_{k}=\left(1+\frac{\hat{\tau}_{6}}{\hat{\beta}_{k}+\epsilon}\right)^{6}$ for $k=1,2,3,4,$}$$

   (2.30)

   where $\epsilon=10^{-40}$.

2. 2.

   Construct $H^{(i,j)}(x,y)$ by

   $$\begin{split}H^{(i,j)}(x,y)=q_{0}(x,y),~{}~{}~{}\text{if}~{}~{}~{}\min_{k}\hat{\kappa}_{k}>C_{T},\end{split}$$

   (2.31)

   otherwise

   $$\begin{split}H^{(i,j)}(x,y)=q_{K}(x,y)~{}~{}~{}\text{with $K$ being the index such that $\hat{\kappa}_{K}=\max_{k}\hat{\kappa}_{k}$}.\end{split}$$

   (2.32)

   Here, $C_{T}$ is chosen as $10^{-3}$.

In particular, if $\widetilde{\overline{v}}^{n}_{i,j}=\widetilde{\overline{v}}^{n}_{i,j,0}$, $\widetilde{\overline{w}}^{n}_{i,j}=\widetilde{\overline{w}}^{n}_{i,j,0}$ (with similar notation), and $H^{(i,j)}(x,y)=q_{0}(x,y)\quad\forall i,j,$ we call such a reconstruction linear reconstruction.

To evaluate the right-hand side of (2.5), it remains to provide the approximations of $\{I_{i,j}^{\star}\}$ and $w(x,y,t^{n})$. For a given index $(i,j)$, the cubic-curved quadrilateral upstream cell, $\widetilde{I}^{\star}_{i,j}$, and the cubic polynomial $\widetilde{w}(x,y)$ on $\widetilde{I}^{\star}_{i,j}$ are constructed by the following procedure.

1. 1.

   Tracing characteristics backward in time.

   We locate $4\times 4$ Gauss-Legendre-Lobatto (GLL) points on $I_{i,j}$. We determine the characteristic feet of these GLL points by solving (2.2) at $t=t^{n}$ (see Figure 2.3). In practice, we solve the ODEs (2.2) by a fourth-order Runge-Kutta (RK) method.

   \begin{overpic}[scale={0.07}]{Figures//v}\put(49.0,4.0){$i$} \put(4.0,44.5){$j$} \end{overpic}

   \begin{overpic}[scale={0.07}]{Figures//v_star} \put(49.0,4.0){$i$} \put(4.0,44.5){$j$} \end{overpic}

2. 2.

   Constructing the edges of the cubic-curved quadrilateral upstream cells.

   The edges of the cubic-curved quadrilateral upstream cells are constructed by a cubic interpolation procedure (see [42]). In particular, we prefer to use a parametric form to present any cubic-curved edge:

   $$\begin{cases}x(\xi)=x_{a}\xi^{3}+x_{b}\xi^{2}+x_{c}\xi+x_{d},\\ y(\xi)=y_{a}\xi^{3}+y_{b}\xi^{2}+y_{c}\xi+y_{d},~{}~{}~{}\xi\in[-1,1].\end{cases}$$

   (2.33)

3. 3.

   Constructing $\widetilde{w}(x,y)$.

   We denote the sixteen GLL points in Figure 2.3 (a) by $\{v_{k}\}$ and denote the corresponding characteristic feet by $\{v_{k}^{\star}\}$. By the adjoint problem (2.3), we have

   $$w(x(v_{k}^{\star}),y(v_{k}^{\star}))=W(x(v_{k}),y(v_{k}))~{}~{}~{}\text{for}~{}~{}k=1,2,\ldots,16.$$

   (2.34)

   Hence, by a standard least square procedure, we can find a cubic polynomial $\widetilde{w}(x,y)\in P^{3}$.

For numerically evaluating the right-hand side of (2.5), we integrate the piecewise polynomial $H^{n}(x,y)\widetilde{w}(x,y)$ over each $\widetilde{I}_{i,j}^{\star}$, which may cross different Eulerian background cells. Hence, $\widetilde{I}_{i,j}^{\star}$ is first clipped into curved polygons such that $\widetilde{I}_{i,j}^{\star}=\cup_{(p,q)}(\widetilde{I}^{\star}_{i,j}\cap I_{p,q})$ and the integrand is smooth in each polygon. For conciseness, the curved polygons $\left\{\widetilde{I}^{\star}_{i,j}\cap I_{p,q}\right\}$ are denoted by $\left\{\widetilde{I}_{i,j;p,q}^{\star}\right\}$. In particular, we are concerned about the information along the edges of $\left\{\widetilde{I}_{i,j;p,q}^{\star}\right\}$. The edges of $\left\{\widetilde{I}_{i,j;p,q}^{\star}\right\}$, overlapping $\partial\widetilde{I}_{i,j}^{\star}$, with counterclockwise direction with respect to $\widetilde{I}_{i,j}^{\star}$ are denoted by $\{\mathcal{L}_{i,j;p,q}^{k}\}$ (see Figure 2.4 (a)). We call $\{\mathcal{L}_{i,j;p,q}^{k}\}$ the outer integral segments of the upstream cell $\widetilde{I}_{i,j}^{\star}$. Similarly, the edges of $\left\{\widetilde{I}_{i,j;p,q}^{\star}\right\}$, overlapping mesh lines, with counterclockwise direction with respect to the corresponding curved polygons are denoted by $\{\mathcal{S}_{i,j;p,q}^{k}\}$ (see Figure 2.4 (b)). We call $\{\mathcal{S}_{i,j;p,q}^{k}\}$ the inner integral segments of $\widetilde{I}_{i,j}^{\star}$. For implementation, we refer to [42] for more details.