Algorithm xxxx: HiPPIS A High-Order Positivity-Preserving Mapping Software for Structured Meshes

Both the DBI and PPI methods rely on the Newton polynomial (doi:10.1137/0912034, ; 10.2307/2004888, ) representation to build interpolants that are positive or bounded by the data values. The ability to adaptively choose stencil points to construct the interpolation, as in ENO methods (HARTEN19973, ), is the key feature employed to develop the data-bounded and positivity-preserving interpolants.

Consider a 1D mesh defined as follows:

where $x_{i-J}<\cdots<x_{i}<x_{i+1}<\cdots<x_{i+L}$, and $\{u_{i-J},\cdots,u_{i+L}\}$ is the set of data values associated with the mesh points in Equation (1). The subscripts $J$, $L$, $i,\in\mathbb{N}_{0}=\mathbb{N}\cup\{0\}$, and $x_{k}$, $u_{k}\in\mathbb{R}$ for $i-J\leq k\leq i+L$. The DBI and PPI procedure starts by setting the initial stencil $\mathcal{V}_{0}$,

The stencil $\mathcal{V}_{0}$ in Equation (2) is expanded by successively appending a point to the right or left of $\mathcal{V}_{j}$ to form $\mathcal{V}_{j+1}$. Once the final stencil $\mathcal{V}_{n-1}$ is obtained, the interpolant of degree $n$ defined on $I_{i}=\{x_{i},x_{i+1}\}$ can be written as

where $\pi_{0,i}(x)=(x-x_{i}),\pi_{1,i}(x)=(x-x_{i})(x-x_{1}^{e}),\cdots$ are the Newton basis functions. $x_{j}^{e}$ is the point added to expand the stencil $\mathcal{V}_{j-2}$ to $\mathcal{V}_{j-1}$ and can be explicitly expressed as

The divided differences are recursively defined as follows:

The polynomial $U_{n}(x)$ can be compactly expressed as

$S_{n}(x)$ in Equation (4) is defined as

where $s$, $t_{j}$, and $d_{j}$ are expressed as follows:

$s$ and $d_{j}$ in Equations (5), (6), and (8) are defined such that $s\in[0,1]$ and $d_{j}\geq 0$. The positivity-preserving and data-bounded interpolants are obtained by imposing some bounds on $\bar{\lambda}_{j}$, defined as

For a given interval inside a mesh, the DBI and PPI polynomial interpolant is constructed by adaptively selecting points near the target interval to build an interpolation stencil and polynomial that together meet the requirements for positivity or data-boundedness. Requiring positivity alone can lead to large oscillations and extrema that degrade the approximation. Positivity alone does not restrict how much the interpolant is allowed to grow beyond the data values. The large oscillations can be removed with the PCHIP, MQSI, and DBI methods. However, in the case where a given interval $I_{i}$ has a hidden extremum, PCHIP, MQSI, and DBI will truncate the extremum. As in  (berzins2010nonlinear, ; sekora2009extremumpreserving, ), the interval $I_{i}$ has an extremum when two of the three divided differences $\sigma_{i-1}=U[x_{i-1},x_{i}]$, $\sigma_{i}=U[x_{i},x_{i+1}]$, and $\sigma_{i+1}=U[x_{i+1},x_{i+2}]$ of neighboring intervals are of opposite signs. The constrained PPI algorithm addresses these limitations by allowing the constructed interpolant to grow beyond the data values but not produce extrema that are too large.

The positive polynomial interpolant is constrained as follows:

The bounds $u_{min}$ and $u_{max}$ in Equation (10) are defined as

The parameters $\Delta_{min}$ and $\Delta_{max}$ in Equation (11) are positive, and the data-bounded interpolant is obtained for $\Delta_{min}$=$\Delta_{max}=0.0$. These parameters are chosen according to

and

The positive parameters $\epsilon_{0}$ and $\epsilon_{1}$, used for intervals with and without extrema, respectively, are introduced to adjust $\Delta_{min}$ and $\Delta_{max}$. This work extends the bounds in (Ouermi2022, ) by introducing the parameter $\epsilon_{1}$ to allow for more flexibility on how to bound the interpolants in cases where an extremum is detected. The choice for the positive parameters $\epsilon_{0}$ and $\epsilon_{1}$ depends on the underlying function and the input data used for the approximation. As both $\epsilon_{0}$ and $\epsilon_{1}$ get smaller, the upper and lower bounds get tighter and the PPI method converges to the DBI method. The choices for $\epsilon_{0}$ and $\epsilon_{1}$ are further discussed in Section 3.2. In Equation (12), the interval $I_{i}$ has a local maximum if $\sigma_{i-1}\sigma_{i+1}<0$ and $\sigma_{i-1}<0$. Correspondingly, in Equation (13), the interval $I_{i}$ has a local minimum if $\sigma_{i-1}\sigma_{i+1}<0$ and $\sigma_{i-1}>0$. In both Equations (12) and (13), the type of extremum is ambiguous if $\sigma_{i-1}\sigma_{i+1}$, and $\sigma_{i-1}\sigma_{i}<0$.

Equation (10) is equivalent to bounding $S_{n}(x)$ as follows:

where the factors $m_{\ell}$ and $m_{r}$ in Equation (14) are expressed as

1. (1)

   : $u_{i+1}>u_{i}$

   (15)

   $$m_{\ell}=min\bigg{(}0,\frac{u_{min}-u_{i}}{u_{i+1}-u_{i}}\bigg{)},\textrm{ and }m_{r}=max\bigg{(}1,\frac{u_{max}-u_{i}}{u_{i+1}-u_{i}}\bigg{)}$$

2. (2)

   : $u_{i+1}<u_{i}$

   (16)

   $$m_{\ell}=min\bigg{(}0,\frac{u_{max}-u_{i}}{u_{i+1}-u_{i}}\bigg{)},\textrm{ and }m_{r}=max\bigg{(}1,\frac{u_{min}-u_{i}}{u_{i+1}-u_{i}}\bigg{)}.$$

The DBI method can be recovered from the PPI methods by setting $m_{\ell}=0$ and $m_{r}=1$. For $u_{i}=u_{i+1}$, $m_{\ell}$, $m_{r}$ and $U_{n}(x)$ as written in Equations (15), (16), and (4) are not defined. This limitation is addressed by re-writing $U_{n}(x)$ as

where $S_{n}(x)$ is expressed as follows:

The summation starts at $j=1$ because the linear term $\frac{u_{i+1}-u_{i}}{x_{i+1}-x_{i}}(x-x_{i})=0$. Let

$\bar{\lambda}_{j}$ in this context is defined as

For $u_{i}=u_{i+1}$, the parameters $m_{\ell}$ and $m_{r}$ in Equation (14) are then defined according to

1. (1)

   : $U[x_{1}^{l},\cdots,x_{1}^{r}]>0$

   (17)

   $$m_{\ell}=min\bigg{(}0,\frac{u_{min}-u_{i}}{w}\bigg{)},\textrm{ and }m_{r}=max\bigg{(}1,\frac{u_{max}-u_{i}}{w}\bigg{)}$$

2. (2)

   : $U[x_{1}^{l},\cdots,x_{1}^{r}]<0$

   (18)

   $$m_{\ell}=min\bigg{(}0,\frac{u_{max}-u_{i}}{w}\bigg{)},\textrm{ and }m_{r}=max\bigg{(}1,\frac{u_{min}-u_{i}}{w}\bigg{)}.$$

For $U[x_{i},x_{i+1}]=U[x_{1}^{l},\cdots,x_{1}^{r}]=0$, the data $u_{i-1}$, $u_{i}$, $u_{i+1}$, and $u_{i+2}$ have the same value ($u_{i-1}=u_{i}=u_{i+1}=u_{i+2}$). In this case, the algorithm approximates the function in the interval $I_{i}$ with a linear interpolant.

The positivity-preserving result in Equation (10) is obtained by successively imposing bounds on the quadratic, cubic, and higher order terms in the expression of $S_{n}(x)$ defined in Equation (5). The reconstruction procedure begins by considering the linear and quadratic terms from $S_{n}(x)$ in Equation (5) and imposing the following bounds:

Equation (19) can be reorganized to obtain

Noting that $\frac{1}{s(s-1)}\leq-4$, $\frac{1}{s}\geq 1$, and $\frac{1}{s-1}\leq-1$, we obtain

The bounds from Equation (20) are extended to bound the cubic form by requiring that what multiplies $\bar{\lambda}_{1}$ must fit into the inequality in Equation (20). Thus, for the cubic case, Equation (20) becomes

When $t_{2}$ defined in Equation (7) is negative, $s-t_{2}$ has a maximum value at $s=1$ and a minimum value at $s=0$. $\bar{\lambda}_{2}$ is then bounded by

When $t_{2}$ is positive, $\frac{1}{1-t_{2}}$ is substituted by $\frac{1}{-t_{2}}$ and the inequalities $\leq$ with $\geq$ and vice versa are swapped.

This procedure is continued to quartic and higher order interpolants to produce the recursive expression for the bounds on $\bar{\lambda}_{j}$ for the PPI and DBI methods as follows:

$B_{1}^{-}$ and $B_{1}^{+}$ are defined as $-d_{1}$ and $d_{1}$ for the DBI method, whereas for the PPI method, they are defined as $(-4(m_{r}-1)-1)d_{1}$ and $(-4m_{\ell}+1)d_{1}$, respectively. We refer the reader to Theorems 1 and 2 in  (Ouermi2022, ) for more details on the mathematical foundation used to build the positivity-preserving software.

The algorithms provide the necessary elements to construct the data-bounded or positive interpolants. Rogerson and Meiburg  (rogerson1990numerical, ) showed that the ENO reconstruction can lead to a left- or right-biased stencil that causes stability issues when used to solve hyperbolic equations. Shu  (shu1990numerical, ) addressed this limitation by introducing a bias coefficient used to target a preferred stencil. As indicated in  (Ouermi2022, ), the left- and right-biased stencil can fail to recover hidden extrema. For a given interval $I_{i}$, the left- and right-biased stencil does not include the points $x_{i-1}$ or $x_{i+1}$, respectively. Algorithm I addresses these limitations by extending the algorithm in  (Ouermi2022, ) to introduce more options for the adaptive stencil selection process described below. In addition to the symmetry-based points selection in  (Ouermi2022, ), Algorithm I includes ENO-type and locality-based point selection processes.

At any given step $j$, the next point inserted into $\mathcal{V}_{j}$ can be to the left or right. Let $x_{p}$ and $x_{q}$ be the mesh points immediately to the left and right of $\mathcal{V}_{j}$, respectively.

We define $\bar{\lambda}_{j+1}^{-}$ and $\bar{\lambda}_{j+1}^{+}$ as follows:

The terms $\bar{\lambda}_{j+1}^{-}$ and $\bar{\lambda}_{j+1}^{+}$ in Equation 22 correspond to the case where the stencil inserted is to the left and right, respectively. Given $\mathcal{V}_{j}$, let $\mu_{j}^{l}$ be the number of points to the left of $x_{i}$ and $\mu_{j}^{r}$ the number of points to the right. Algorithm I extends the algorithm in  (Ouermi2022, ) by introducing a user-supplied parameter $st$ used to guide the procedure for stencil construction. In the cases where adding both $x_{p}$ (to the left) or $x_{q}$ (to the right) are valid, the algorithm makes the selection based on the three cases below:

• •

  If $st=1$ (default), the algorithm chooses the point with the smallest divided difference, as in the ENO stencil.

• •

  If $st=2$, the point to the left of the current stencil is selected if the number of points to the left of $x_{i}$ is smaller than the number of points to the right. Similarly, the point to the right is selected if the number of points to the right of $x_{i}$ is smaller than the number of points to the left. When both the number of points to the right and left are the same, the algorithm chooses the point with the smallest $\bar{\lambda}_{j+1}$.

• •

  If $st=3$, the algorithm chooses the point that is closest to the starting interval $I_{i}$. It is important to prioritize the closest points in cases where the intervals surrounding $I_{i}$ vary significantly in size. These variations are found in computational problems for which different resolutions are used for different parts of the domain.

Algorithm II describes the 1D DBI and PPI methods built using the mathematical framework in Section 2 and Algorithm I. Algorithm II further extends the constraints in (Ouermi2022, ) by introducing the user-supplied positive parameter $\epsilon_{1}$ that is used to impose upper and lower bounds on the interpolants according to Equations (12) and (13). The positive parameters $\epsilon_{0}$ and $\epsilon_{1}$ are used for intervals without and with an extremum, respectively. The user-supplied parameter $im$ is used to choose between the DBI and PPI methods. Algorithm III and Algorithm IV describe the extension from 1D to 2D and 3D, respectively. Both Algorithm III and IV are constructed by successively applying Algorithm II to each dimension. Given that the DBI and PPI methods are nonlinear, the order in which Algorithm II is used can lead to different approximation results. In this paper and in (Ouermi2022, ), the 1D DBI and PPI are first applied to $x$, then the $y$ dimension, and finally the $z$ dimensions, as indicated in Algorithm III and IV. In Algorithm III and IV, the input and intermediate data values are modified only by Algorithm I, which preserves data-boundedness or positivity. Therefore, the resulting solutions from the 2D and 3D extensions will preserve data-boundedness and positivity. Similar to Algorithm II, the choices of parameters $st$, $\epsilon_{0}$, and $\epsilon_{1}$ influence the quality of the approximation in Algorithm III and IV. In the 1D, 2D, and 3D cases, the choice for parameters $st$, $\epsilon_{0}$, and $\epsilon_{1}$ is dependent on the data. In the case of 2D and 3D, applying the 1D PPI method (Algorithm II) along the $x$ and/or $y$ dimensions may introduce oscillations that could be amplified when applying the 1D PPI in the subsequent dimensions. These oscillations can be significantly reduced with small values for $\epsilon_{0}$ and $\epsilon_{1}$. The parameters $\epsilon_{0}$ and $\epsilon_{1}$ should be chosen to be small enough such that hidden extrema are recovered without introducing new large oscillations. In Algorithm II and IV, $\mathcal{M}_{\square x\times\square y}$ and $\mathcal{M}_{\square x\times\square y\square z}$ represent 2D and 3D meshes obtained by taking the tensor product of the 1D mesh along the $x$ and $y$ dimensions for the 2D mesh, and along the $x$, $y$, and $z$ dimensions for the 3D mesh. The square $\square$ represents $n$ or $m$.

Algorithm I
Input: $\mu_{j}^{l}$, $\mu_{j}^{l}$, $x_{p}$, $x_{i}$, $x_{q}$, $x_{i+1}$, $U[x_{p},\cdots,x_{j}^{r}]$, $U[x_{j}^{l},\cdots,x_{q}]$ $\bar{\lambda}_{j+1}^{-}$, $\bar{\lambda}_{j+1}^{+}$, and $st$.

1. (1)

   if $st=1$

   • •

     if $|U[x_{p}\cdots,x_{j}^{r}]|<|U[x_{j}^{l},\cdots,x_{q}]|$, then insert a new stencil point to the left;

   • •

     else if $|U[x_{p}\cdots,x_{j}^{r}]|>|U[x_{j}^{l},\cdots,x_{q}]|$, then insert a new stencil point to the right;

   • •

     else insert a new stencil point to the right if $|\bar{\lambda}_{j+1}^{-}|\geq|\bar{\lambda}_{j+1}^{+}|$; otherwise, insert a new point to left;

2. (2)

   if $st=2$

   • •

     if $\mu_{j}^{l}<\mu_{j}^{r}$, then insert a new stencil point to the left;

   • •

     else if $\mu_{j}^{l}>\mu_{j}^{r}$, then insert a new stencil point to the right;

   • •

     else insert a new stencil point to the right if $|\bar{\lambda}_{j+1}^{-}|\geq|\bar{\lambda}_{j+1}^{+}|$; otherwise, insert a new point to left;

3. (3)

   else $st=3$

   • •

     if $|x_{p}-x_{i}|<|x_{q}-x_{i+1}|$, then insert a new stencil point to the left;

   • •

     else if $|x_{p}-x_{i}|>|x_{q}-x_{i+1}|$, then insert a new stencil point to the right;

   • •

     else insert a new stencil point to the right if $|\bar{\lambda}_{j+1}^{-}|\geq|\bar{\lambda}_{j+1}^{+}|$; otherwise, insert a new point to left;

Algorithm II (1D)
Input: $\{x_{i}\}_{i=0}^{n}$, $\{u_{i}\}_{i=0}^{n}$, $\{\tilde{x}_{i}\}_{i=0}^{\tilde{n}}$, $d$, $st$ $\epsilon_{0}$, $im$, and $\epsilon_{1}$. Output: $\{\tilde{u}_{i}\}_{i=0}^{\tilde{n}}$.

1. (1)

   Select an interval $[x_{i},x_{i+1}]$. Let $\mathcal{V}_{0}=\{x_{i},x_{i+1}\}=\{x_{0}^{l},x_{0}^{r}\}$.

2. (2)

   If $\sigma_{i-1}\sigma_{i+1}<0$ or $\sigma_{i-1}\sigma_{i}<0$, then the interval $I_{i}$ has a hidden local extremum. For the boundary intervals, we assume that the divided differences to the left and right have the same sign.

3. (3)

   Compute $u_{min}$ and $u_{max}$ using Equations (11), (12), and (13).

4. (4)

   Compute $m_{r}$ and $m_{\ell}$ based on Equations (15) and (16) or Equations (17) and (18). For DBI, set $m_{r}=1$ and $m_{\ell}=0$.

5. (5)

   Given a stencil $\mathcal{V}_{j}$,

   • •

     if $B^{-}_{j+1}\leq\bar{\lambda}_{j+1}^{+}\leq B_{j+1}^{+}$ and $B^{-}_{j+1}\leq\bar{\lambda}_{j+1}^{-}\leq B_{j+1}^{+}$, choose the point to add based on Algorithm I

   • •

     else if $B^{-}_{j+1}\leq\bar{\lambda}_{j+1}^{-}\leq B_{j+1}^{+}$, then insert a new stencil point to the left;

   • •

     else if $B^{-}_{j+1}\leq\bar{\lambda}_{j+1}^{+}\leq B_{j+1}^{+}$, then insert a new stencil point to the right;

6. (6)

   This process (Steps $3$) iterates until the halting criterion that the ratio of divided differences lies outside the required bounds stated above or the stencil has $d+1$ points, with $d$ being the target degree for the interpolant.

7. (7)

   Evaluate the final interpolant $U^{l}(x)$ (for DBI) or $U^{p}(x)$ (for PPI) at the output points $\tilde{x}_{i}$ that are in $I_{i}$.

8. (8)

   Repeat Steps $1$–$7$ for each interval in the input 1D mesh.

Algorithm III (2D)
Input:$\{x_{i}\}_{i=0}^{nx}$,$\{y_{i}\}_{i=0}^{ny}$ $\{u_{i,j}\}_{i,j=0}^{nx,ny}$, $d$, $st$ $\epsilon_{0}$, $\{\tilde{x}_{i}\}_{i=0}^{mx}$, $\{\tilde{y}_{j}\}_{j=0}^{my}$, $im$, and $\epsilon_{1}$. Output: $\{\tilde{u}_{i,j}\}_{i,j=0}^{mx,my}$.

1. (1)

   Interpolate from the mesh $\mathcal{M}_{nx\times ny}$ to $\mathcal{M}_{mx\times ny}$ by applying the Algorithm II (1D) along the x dimension for each index $j$ along the y direction. More precisely, for each index $j$ ($0\leq j\leq nx$), Algorithm II (1D) uses $\{x_{i}\}_{i=0}^{nx}$, $\{u_{i,j}\}_{i=0}^{nx}$ to interpolate from $\{x_{i},y_{j}\}_{i=0}^{nx}$ to $\{\tilde{x}_{i},y_{j}\}_{i=0}^{mx}$. The interpolation solution is saved in $\{q_{i,j}\}_{i,j=0}^{mx,ny}$

2. (2)

   Use the solution$\{q_{i,j}\}_{i,j=0}^{mx,ny}$ from step (1) to interpolate from $\mathcal{M}_{mx\times ny}$ to the final output mesh $\mathcal{M}_{mx\times my}$ by applying the Algorithm II (1D) along the y dimension for each index $i$ along the x direction. For each index $i$ ($0<i<mx$), Algorithm II (1D) uses ($\{y_{j}\}_{j=0}^{ny}$, $\{q_{i,j}\}_{j=0}^{ny}$) to interpolate from $\{\tilde{x}_{i},y_{j}\}_{j=0}^{ny}$ to $\{\tilde{x}_{i},\tilde{y}_{j}\}_{j=0}^{my}$. The interpolation final results are saved in $\{\tilde{u}_{i,j}\}_{i,j}^{mx,my}$

Algorithm IV (3D)
Input:$\{x_{i}\}_{i=0}^{nx}$,$\{y_{i}\}_{i=0}^{ny}$, $\{z_{k}\}_{k=0}^{nz}$, $\{u_{i,j,k}\}_{i,j,k=0}^{nx,ny,nz}$, $d$, $st$ $\epsilon_{0}$, $\{\tilde{x}_{i}\}_{i=0}^{mx}$, $\{\tilde{y}_{j}\}_{j=0}^{my}$, $\{\tilde{z}_{k}\}_{k=0}^{mz}$, $im$, and $\epsilon_{1}$. Output: $\{\tilde{u}_{i,j,k}\}_{i,j=0}^{mx,my,mz}$.

1. (1)

   Interpolate from the mesh $\mathcal{M}_{nx\times ny\times nz}$ to $\mathcal{M}_{mx\times ny\times nz}$ by applying the Algorithm II (1D) along the x dimension for each pair $j,k$. More precisely, for each each pair $j,K$ ($0\leq j\leq nx$ and $0\leq k\leq nz$), Algorithm II (1D) uses $\{x_{i}\}_{i=0}^{nx}$, $\{u_{i,j,k}\}_{i=0}^{nx}$ to interpolate from $\{x_{i},y_{j},z_{k}\}_{i=0}^{nx}$ to $\{\tilde{x}_{i},y_{j},z_{k}\}_{i=0}^{mx}$. The interpolation solution is saved in $\{q_{i,j,k}\}_{i,j,k=0}^{mx,ny,nz}$

2. (2)

   Use the solution$\{q_{i,j,k}\}_{i,j,k=0}^{mx,ny,nz}$ from step (1) to interpolate from $\mathcal{M}_{mx\times ny\times nz}$ to $\mathcal{M}_{mx\times my\times nz}$ by applying the Algorithm II (1D) along the $y$ dimension. For each pair $i,k$ ($0\leq i\leq mx$ and $0\leq k\leq nz$), Algorithm II (1D) uses ($\{y_{j}\}_{j=0}^{ny}$, $\{q_{i,j,k}\}_{j=0}^{ny}$) to interpolate from $\{\tilde{x}_{i},y_{j},z_{k}\}_{j=0}^{ny}$ to $\{\tilde{x}_{i},\tilde{y}_{j},z_{k}\}_{j=0}^{my}$. The interpolation results are saved in $\{\tilde{g}_{i,j,k}\}_{i,j}^{mx,my,nz}$.

3. (3)

   Use the solution$\{q_{i,j,k}\}_{i,j,k=0}^{mx,ny,nz}$ from step (1) to interpolate from $\mathcal{M}_{mx\times ny\times nz}$ to $\mathcal{M}_{mx\times my\times nz}$ by applying the Algorithm II (1D) along the $y$ dimension. For each pair $i,j$ ($0\leq i\leq mx$ and $0\leq j\leq my$), Algorithm II (1D) uses ($\{z_{k}\}_{k=0}^{nz}$, $\{g_{i,j,k}\}_{k=0}^{nz}$) to interpolate from $\{\tilde{x}_{i},\tilde{y}_{j},z_{k}\}_{k=0}^{nz}$ to $\{\tilde{x}_{i},\tilde{y}_{j},\tilde{z}_{k}\}_{k=0}^{mz}$. The interpolation final results are saved in $\{\tilde{u}_{i,j,k}\}_{i,j}^{mx,my,mz}$.

The DBI and PPI software implementation is guided by the algorithms described above. HiPPIS is available at https://github.com/ouermijudicael/HiPPIS. The software can be organized into four major parts: (1) computation of divided differences, (2) calculations of upper and lower bounds for each interval, (3) a stencil construction procedure, and (4) 1D, 2D, and 3D DBI and PPI implementations.

The divided differences are essential to the DBI and PPI methods because they are used in the calculations of $\bar{\lambda}_{j}$ and the stencil selection process. The divided differences are computed using the standard recurrence form in Equation (2.1) and stored in a table of dimension $n\times(d+1)$ where $d$ is the maximum polynomial degree for each interpolant. Given that the maximum degree is $d$, it is sufficient to consider the $d+1$ divided differences for the stencil selection process and the construction of the final polynomial interpolant for each interval.

The bounds on each interpolant are obtained from Equation (11), (12), and (13) where the positive parameters $\epsilon_{0}$ and $\epsilon_{1}$ are user-supplied values used to adjust the bounds for the interval with and without extremum, respectively. The adjustment focuses on removing large oscillations as much as possible while still allowing high-degree polynomial interpolants that meet the positivity requirements.

The stencil selection process requires the computation of $B_{j}^{+}$ and $B_{j}^{-}$, which are both dependent on $d_{j}$, $t_{j}$, and $\bar{\lambda}_{j}$. The stencil $\mathcal{V}_{j}$ is constructed from $\mathcal{V}_{j-1}$ by appending a point to the left or right of $\mathcal{V}_{j-1}$. When appending to either the right or left meets the requirements for positivity, the software offers three possible options for choosing from both points that can be set by the user. The first and default option ($st=1$) chooses the stencil with the smallest divided difference, similar to the ENO-like approach. The second option ($st=2$) prioritizes the choice that makes the stencil more symmetric around $x_{i}$. The third option ($st=3$) chooses the point closest to the starting interval $I_{i}$, thus prioritizing locality.

The 1D DBI and PPI methods use (1), (2), and (3) as building blocks to construct the final approximation. Once the final stencil has been selected, the interpolant is built using a Newton polynomial representation and then evaluated at the corresponding output points. The Newton polynomial is used here because its coefficients/divided differences are available. The 2D and 3D implementations successively use the 1D version along each dimension to construct the final approximation on uniform and nonuniform structured meshes.

The interfaces for the 1D, 2D, and 3D DBI and PPI subroutines are designed to be similar to widely used interfaces for polynomial interpolation such as PCHIP, and can be incorporated into larger application codes. The interfaces require

• •

  the input mesh points and the data values associated with those points,

• •

  the maximum polynomial degree to be used for each interpolant,

• •

  the interpolation method to be used (DBI or PPI), and

• •

  the output mesh points.

Listing 1 shows examples of how to use the 1D, 2D, and 3D interfaces for DBI and PPI in HiPPIS. The variables x, y, and z are 1D arrays used to define the input meshes, and xout, yout, and zout are used to define the output meshes. The variables v, v2D, and v3D correspond to the input data values associated with the input meshes. The parameters d and im (1, or 2) are used to indicate the target polynomial degree and the interpolation method to be used. For DBI and PPI, the parameter im is set to 1 and 2, respectively. The parameters $st$, $\epsilon_{0}$, and $\epsilon_{1}$ are optional parameters that are set to $3$, $0.01$, and $1$ by default, as explained below. The choice of the optional parameters depends on the underlying function and the input data.

In problems for which different resolutions are used for different parts of the computational domain, st=3 is a preferable choice. The algorithm prioritizes the closest points to the starting interval $I_{i}$ if st=3. This choice is particularly important in regions where the size of the intervals varies significantly. For cases when smoothness is the primary goal, st=1 is a suitable choice. For st=1, the Algorithm I prioritizes smoothness by choosing the points with the smallest divided differences during the stencil construction process. Both the st=1 and st=3 can lead to a left- or right-biased stencil. In these instances, st=2 can be used to remove the bias. For st=2, the algorithm prioritizes a symmetric stencil. The default value of st is set to 3 because the examples in this study indicate that st=3 leads to better approximations compared to st=1 or 2 and locality is often a highly desired property in many computational problems.

The positive parameters $\epsilon_{0}$ and $\epsilon_{1}$ are used to bound the interpolants for the intervals with and without extrema, respectively. The configurations in  (Ouermi2022, ) and  (arkivtajo, ) correspond to setting the parameters $\epsilon_{0}$ and $\epsilon_{1}$ to the default values of $0.01$ and $1$, respectively. The values of $\epsilon_{0}$ and $\epsilon_{1}$ are chosen such that the lower and upper bounds on each interpolant are relaxed enough to allow for a high-order polynomial that does not introduce undesirable oscillations. For profiles that are prone to oscillation such as the logistic functions, it is important to choose small values for $\epsilon_{0}$ and $\epsilon_{1}$. For $N\times N=17\times 17$, the approximation leads to large oscillations if $\epsilon_{0}$ and $\epsilon_{1}$ are greater than $10^{-4}$. For intervals without extrema, it is important to keep $\epsilon_{0}$ small to not introduce new extrema. For the intervals with extrema, $\epsilon_{1}$ needs to be large enough to allow for recovery of hidden extrema but small enough to not cause undesired large oscillations. This is very challenging given that the sizes of the peaks are not known a priori. The default values of $\epsilon_{1}=1$ are such that the interpolant maximum value is twice $max(u_{i},u_{I+1})$. This default value of one is sufficient for the modified Runge and TWP-ICE examples. However, in the case of BOMEX, smaller values of $\epsilon_{1}\leq 10^{-5}$ are required to remove undesired oscillations. In practice, it is prudent to start with a small value for $\epsilon_{0}$ and $\epsilon_{1}$ and increase them as needed if the approximation fails to recover hidden extrema or uses low-degree polynomial interpolants.

Figure 1 is a diagram of the different components of the main module of HiPPIS. The function divdiff(…) is used to calculate the divided differences needed for Algorithms I and II. Once the final stencil is constructed, the function newtonPolyVal(…) is used to build and evaluate the positive interpolant at the corresponding output points. The major part of the data-boundedness and positivity preservation including Algorithms I and II is in the function adaptiveInterpolation1D(…). This function is used for the 1D approximation or mapping problems and depends on the function divdiff(…) and newtonPolyVal(…). The functions adaptiveIterpolation2D(…) and adaptiveInterpolation3D(…) use adaptiveInterpolation1D(..) to construct the data-bounded or positive polynomial approximations on 2D and 3D structured tensor product meshes, respectively. The interfaces for the 1D, 2D, and 3D interpolations, in bold, require the parameter $im$, which is used to indicate the interpolation method chosen. For the DBI and PPI methods, the parameter $im$ is set using 1 and 2, respectively. HiPPIS does not allow for any other choices for the parameter $im$.