<figure><img src="image_1.png" /></figure>

HAL
open science

# A novel optimization technique for explicit finite-difference schemes with application to AeroAcoustics

Guilherme Cunha, Stéphane Redonnet

► To cite this version:

Guilherme Cunha, Stéphane Redonnet. A novel optimization technique for explicit finite-difference schemes with application to AeroAcoustics. International Journal for Numerical Methods in Fluids, 2015, 78 (4), pp.189-216. 10.1002/fld.4010. hal-01226262

HAL Id: hal-01226262

https://hal.science/hal-01226262

Submitted on 4 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A novel optimization technique for explicit finite-difference
schemes with application to AeroAcoustics

Guilherme Cunha\*,† and Stéphane Redonnet

ONERA, CFD and Aeroacoustics Department, Châtillon, France

The present paper addresses the optimization of finite-difference schemes when these are to be used for numerically approximating spatial derivatives in aeroacoustics evolution problems. With that view in mind, finite-difference operators are firstly detailed from a theoretical point of view. Secondly, time, the way such operators can be optimized in a spectral-like sense is recalled, before the main limitations of such an optimization are highlighted. This leads us to propose an alternative optimization approach of innovative character. Such a novel optimization technique consists of enhancing the scheme's formal accuracy through a minimization of its leading-order truncation error. This so-called *intrinsic* optimization procedure is first detailed, before it is thoroughly analyzed, from both a theoretical and a practical point of view. The second part of the paper focuses on two particular intrinsically optimized schemes, which are carefully assessed via a direct comparison against their standard and/or spectral-like optimized counterparts, such a comparative exercise being conducted utilizing several academic test cases of increasing complexity. There, it is shown how intrinsically optimized schemes indeed constitute an advantageous alternative to either the standard or the spectral-like optimized ones, being allotted with both (i) the better scalability of the former scheme with respect to grid convergence effects when the grid density increases and (ii) the higher accuracy of the latter scheme when the discretization level becomes marginal. Thanks to that, such intrinsically optimized schemes offer very good trade-offs in terms of (i) accuracy; (ii) robustness; and (iii) numerical efficiency (CPU cost).

KEY WORDS: aeroacoustics; finite difference; explicit; optimization; aerodynamics; spectral

# 1. INTRODUCTION

Highly accurate numerical schemes for Computational AeroAcoustics (CAA) have seen widespread development over the two last decades. In particular, the need for accurately simulating acoustic waves over long propagation distances and time durations implies that numerical schemes of minimal dissipation and dispersion properties are used. Indeed, unlike for CFD simulations, CAA calculations cannot stand dissipation nor dispersion errors, which may induce a complete degradation of the acoustic signals to be simulated [1].

Consequently, given their higher accuracy and, thus, ability to simulate wavelike signals on coarser grids [2], high-order Finite-Difference (FD) methods have emerged as a valuable alternative to more traditional techniques coming from the world of CFD, for example, the low-order Finite-Volume ones. Nevertheless, making use of such high-order schemes generally results in less efficient calculations because, for a given error tolerance, these schemes cost more than low order ones, in terms of operations per node [3].

<footnote>\*Correspondence to: Guilherme Cunha, Department of CFD and Aeroacoustics, Onera (French Aerospace Lab), BP 72 - Avenue de la Division Leclerc, Châtillon, 92322, France.</footnote>

<footnote>†E-mail: guilherme.cunha@onera.fr</footnote>
With the purpose of relaxing the latter constraint, one can seek to optimize high-order FD schemes so that both their error level remains low enough and their associated CPU cost is kept at a reasonable level [4]. In particular, regarding space derivatives, several authors proposed to optimize FD schemes following a spectral approach, that is, via an optimization procedure that is carried out in the wavenumber space. Indeed, this so-called spectral-like optimization seeks to enhance the scheme's ability to preserve the acoustic dispersion and dissipation properties by minimizing the error it induces on the spectral components of any given harmonic signal. Among other spectral-like optimized schemes, one can mention here the implicit compact [4–8], the Essentially Non-Oscillatory (ENO) [9], the dispersion-relation-preserving [10], and the low dispersion [11] schemes.

In order to achieve such a minimization of the spectral error, however, one or several of the conditions needed for the FD scheme to reach its maximal order-of-accuracy (and, thus, a low level of truncation error) are to be dropped. The point is that, as highlighted in Section 2 of the present paper, such a reduction in their accuracy order may limit the effective range of applicability of spectral-like optimized schemes.

With a view to overcoming the limitations weighting on the spectral-like optimization, an alternative optimization technique was developed, which constitutes the content of the present article. Relying on a minimization of the scheme's leading-order truncation error, this novel and so-called *intrinsic-optimization* approach is thoroughly detailed in the second part of the paper, being described from both a theoretical and practical point of view in Sections 3 and 4, respectively. Then, in Sections 5 and 6, two particular *intrinsically-optimized* schemes are carefully assessed, which is achieved through a direct comparison of the schemes proposed here with their standard and/or spectral-like optimized counterparts, on the basis of various academic test cases of increasing complexity. After that, in Section 7, the advantages offered by the intrinsic optimization approach are discussed, before general conclusions are drawn. Finally, an appendix is provided, which both details the intrinsic optimization procedure and documents a straightforward way to numerically derive intrinsically-optimized schemes, by solely adjusting the values of coefficients characterizing standard FD schemes.

## 2. SPECTRAL-LIKE OPTIMIZATION OF FD SCHEMES FOR SPATIAL DERIVATION

### 2.1. *Spectral-like optimization*

Discrete differential operators based on explicit FD schemes can be optimized so that their resulting integrated dispersion error is reduced as much as possible. The basic principle of such optimization is to replace one or several identities of (A.7), Appendix A, by dedicated relations to be specifically related to the absolute or relative dispersion errors, which norm, *once integrated over a given frequency range*, is then to be minimized. As a result, one recovers an FD scheme whose coefficients $b_{j,1}$ are slightly modified compared with their standard counterpart ones. Here, it can be recalled that numerous spectral-like FD schemes were proposed in the literature, depending on the application domain envisaged (e.g., acoustics and aerodynamics) or the optimization technique proposed by their respective authors [4–8, 10, 11]. Regarding more particularly an application within a CAA context, the most popular spectral-like optimized FD scheme is certainly the so-called dispersion-relation-preserving one, which was originally proposed by Tam & Webb [10], before it was further extended by Bogey & Bailly into the so-called low-dispersion scheme [11].

As previously stated, these particular spectral-like optimized schemes correspond in all points to their standard counterpart of same stencil width, except their sets of coefficients whose values were adjusted through the respective optimization process. As an example, regarding both the centered seven-point and 11-point stencil schemes, Table I provides the set of optimized coefficients (as obtained by Tam & Webb and Bogey & Bailly, respectively), along with those coming originally from their standard counterpart. In the following, the two particular standard schemes introduced here will be referred to as *FDs7p6o* and *FDs11p10o*, respectively, with the first digit standing for the number of stencil points, while the second stands for the accuracy order. Analogously to the previous notation, the two particular spectral-like optimized schemes introduced earlier will be referred to as *FDo7p4o* and *FDo11p4o*, respectively. Finally, concerning the associated sets of
<table><caption>Table I. Stencil coefficients associated with various centered finite-difference derivative schemes, of either spectral-like or standard nature.</caption><thead><tr><th></th><th>FDs7p6o</th><th>FDo7p4o</th><th>FDs11p10o</th><th>FDo11p4o</th></tr></thead><tbody><tr><td>a<sub>0</sub></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>a<sub>1</sub> = -a<sub>-1</sub></td><td>3/4</td><td>0.77088238051822552</td><td>5/6</td><td>0.872756993962</td></tr><tr><td>a<sub>2</sub> = -a<sub>-2</sub></td><td>-3/20</td><td>-0.166705904414580469</td><td>-5/21</td><td>-0.286511173973</td></tr><tr><td>a<sub>3</sub> = -a<sub>-3</sub></td><td>1/60</td><td>0.02084314277031176</td><td>5/84</td><td>0.090320001280</td></tr><tr><td>a<sub>4</sub> = -a<sub>-4</sub></td><td></td><td></td><td>-5/504</td><td>-0.020779405824</td></tr><tr><td>a<sub>5</sub> = -a<sub>-5</sub></td><td></td><td></td><td>1/1260</td><td>0.002484594688</td></tr></tbody></table>

coefficients for more consistency with the literature, in the following, the notation $a_j = b_{j,1}$ is adopted (Appendix A).

## 2.2. On the efficiency of spectral-like optimized schemes

As was shown in [12], compared with its standard counterpart of same stencil width, each one of the two optimized schemes ($FDo7p4o$ and $FDo11p4o$) exhibits a lower dispersion error – at least over the mid-frequency to high-frequency ranges. At lower frequencies, however, this result does not hold anymore, and both optimized schemes appear to be less accurate than the centered standard seven-point sixth-order scheme, $FDs7p6o$. The same observation can be made regarding the respective error induced on the group velocity by each scheme.

This behavior of spectral-like optimized schemes in the low-frequency range is a direct consequence of the *accuracy order loss inherited from the optimization strategy* (which weighs on all spectral-like optimization procedures, because they are based on a formal order-of-accuracy reduction [3]). In particular, as mentioned earlier, when submitted to a spectral-like optimization, a given FD scheme sees its formal order-of-accuracy decreasing, ending up with a lower value than the one characterizing its standard counterpart. As an example, whereas the centered standard 7-point and 11-point schemes are sixth- and tenth-order accurate, respectively, their spectral-like optimized versions (Tam & Webb and Bogey & Bailly schemes, respectively) are both characterized by a fourth order-of-accuracy.

Consequently, and although they are indubitably better suited than standard ones over the mid-frequency to high-frequency ranges, at lower frequencies, spectral-like optimized schemes are characterized by higher errors, in terms of dispersion and group velocity. Although small, once accumulated over long propagation time durations and distances, these effects may become more significant for spectral-like schemes than for standard ones (with, as a possible consequence, a less accurate propagation of the acoustic energy resulting from the increased error made on the group velocity). Therefore, before privileging a spectral-like optimized FD scheme over its standard counterpart, one should preferably check that the spectral errors made on low-frequency contents remain inferior to an acceptable level, with that level to be defined in accordance with the problem to be solved (e.g., propagation time duration and/or distance) [12].

In addition to that, and still compared with their standard counterpart, spectral-like optimized FD schemes exhibit a relative efficiency that is inversely proportional to the grid refinement level. This indeed also comes from the spectral optimization strategy, which seeks to allot the schemes with a *uniform nominal error level over the entire frequency range of the optimization interval*. As a result, although spectral-like optimized schemes effectively outperform their standard counterpart of same stencil width whenever the mesh grid discretization level is marginal, when the resolution increases this tendency is inverted. As a matter of fact, although spectral-like optimized schemes effectively ensure a nominal error level all over their optimization frequency range, the latter error level cannot be easily reduced further, as this can be done for any standard scheme through a simple grid convergence. Therefore, and contrarily to what happens with standard schemes, mesh refinements are relatively inefficient for lowering the scheme’s error level [12]. This may be problematic whenever the overall error level of the scheme used must be monitored down to a given value – such as when multiple length scale problems are solved with the help of adaptive mesh refinement [13–16] and multiple-size-mesh [17, 18] techniques.
As a conclusion, ideally, it would be desirable to make use of spectral-like optimized schemes preferably when the discretization level is marginal, and to employ standard schemes whenever the mesh resolution is higher.

### 3. INTRINSIC OPTIMIZATION TECHNIQUE: THEORETICAL CONSIDERATIONS

Based on all of the above, one could try optimizing a given FD scheme so that (i) it offers the advantages of spectral-like optimized schemes whenever the resolution becomes marginal, whereas (ii) keeping a higher formal order-of-accuracy, that is, offering the same qualities as its standard counterpart when the discretization level is higher. With a view to answering such a need, in the present section, a novel optimization procedure is proposed.

#### 3.1. Intrinsic optimization: Theoretical background

The optimization technique proposed here takes direct benefit from particular properties characterizing the Lagrange interpolation basis. In the current subsection, firstly, a summary of such properties is given.

##### 3.1.1. Lagrange interpolation and approximate derivative.
Let $I_n$ be a set of $n$ ordered integers and $\{x_j\}_{j \in I_n}$ a finite suite of real numbers. The Cauchy Interpolation Error theorem [19] states that, if a given function $f: \mathbb{R} \to \mathbb{R}$ can be continuously differentiated $n$ times, then it can be approximated with a $n^{th}$-order Lagrange polynomial (of $(n-1)^{th}$ degree) such that

$$f(x) = \sum_{j \in I_n} P_j(x)f(x_j) + \frac{\partial_x^n[f](\xi)}{n!} \prod_{j \in I_n}(x - x_j), x \in \mathbb{R}, \quad (1)$$

for some $\xi$ in the interval spanned by $x$ and the interpolation points $\{x_j\}_{j \in I_n}$. The aforementioned expression relies on Lagrange polynomials, which are defined as

$$P_j(x) = \prod_{p \in I_n, p \neq j} \frac{x - x_p}{x_j - x_p} \quad (2)$$

The first term of the right-hand side of (1) can be seen as an approximation of the function $f(\cdot)$ – such approximation being based on a continuous reconstruction of its discrete form. More precisely, let us consider the continuous interpolated function, $\mathcal{I}[f]: \mathbb{R} \to \mathbb{R}$, to be based on the discrete values $f(x_j)$, with $j \in I_n$:

$$\mathcal{I}[f](x) = \sum_{j \in I_n} P_j(x)f(x_j). \quad (3)$$

After differentiation, $\mathcal{I}[f](x)$ becomes

$$\partial_x[\mathcal{I}[f]](x) = \sum_{j \in I_n} \partial_x[P_j](x)f(x_j). \quad (4)$$

##### 3.1.2. Particular case of equally-spaced points.
Let us consider a set of $2N + 1$ discrete points $x_j$ to be equally-spaced and indexed by $I_n = D_{2N+1} = \{-N, -N+1, \dots, 0, \dots, N-1, N\}$. For doing so, let us choose the suite $\{x_j\}_{j \in D_{2N+1}}$ such as $x_j = j\Delta x$. Now, let us evaluate $\partial_x[\mathcal{I}[f]](x)$ at a given point of the referred stencil. For example, if $x = x_0$, (4) becomes

$$\partial_x[\mathcal{I}[f]](x_0) = \frac{1}{\Delta x} \sum_{j \in D_{2N+1}} a_j f(x_j), \quad (5)$$

where $a_j/\Delta x = \partial_x P_j(x_0)$, $\forall j \in D_{2N+1}$ [20]. Here, the $a_j$ terms correspond exactly to the well-known coefficients characterizing the centered standard $2N^{th}$-order accurate FD scheme, associated
with a spatial derivation of $f(.)$ in the direction $x$, as introduced at Appendix A. In particular, the scheme is antisymmetric ($a_{-j} = -a_j, \forall j \in D_{2N+1}$).

At this stage, it can be underlined that (4) constitutes a more general and straightforward way of deriving FD scheme coefficients. Firstly, it avoids solving the linear system resulting from (A.7). Secondly, it is not restricted to the sole particular case of a centered scheme and equally spaced points, encompassing all possible cases (noncentered schemes, nonuniform grids, etc.) [20].

From a more global perspective, it can also be stated that approximating the derivative of a function, $f(.)$, with the help of an FD scheme, is indeed strictly equivalent to derivating the Lagrange approximation, $\mathcal{I}[f](.)$, which can be built up from the function's discrete values:

$$ \partial_x[f](x) \cong \partial_x[\mathcal{I}[f]](x), \qquad (6) $$

with $x \in \{x_j\}_{j \in D_{2N+1}}$.

Differentiating (1) and projecting the latter Lagrange approximation onto the set of equally-spaced points $\{x_j\}_{j \in D_{2N+1}}$ results in a leading order truncation error that is directly proportional to

$$ E_{\text{std}} = \frac{N!N!}{(2N+1)!}, \quad N \in \mathbb{N}. \qquad (7) $$

As explained earlier, a pertinent way of optimizing an FD scheme would consist of minimizing its leading-order truncation error independently of the function to be differentiated. By observing (1), it can be inferred that such a goal could be reached by defining the interpolation (and thus the derivation) points $\{x_j\}_{j \in D_{2N+1}}$, such that the polynomial term in the error expression is minimized, even after differentiation. This is the focus of the next subsection.

## 3.2. Intrinsic optimization: General strategy

At this stage, it can be recalled that the Chebyshev Minimal Amplitude theorem [19] ensures that the unique polynomial to be the smallest maximum in $[-1, 1]$ is $T_n(\tilde{x})/2^{n-1}$, $T_n$ being the $n^{th}$ Chebyshev polynomial of the first kind.

$$ T_n(\tilde{x}) = \cos(n\cos^{-1}(\tilde{x})), n \in \mathbb{N} \text{ and } \tilde{x} \in [-1, 1]. \qquad (8) $$

It can therefore be concluded that, whenever (i) a suite of points $\{\tilde{x}_k\}_{k \in I_n}$ corresponds to the roots of the $n^{th}$ Chebyshev polynomial and (ii) the interval spanned by $\tilde{x}$ and $\{\tilde{x}_k\}_{k \in I_n}$ is included in $[-1, 1]$, then

$$ \Pi_{j \in I_n}(\tilde{x} - \tilde{x}_j) = \frac{T_n(\tilde{x})}{2^{n-1}}, n \in \mathbb{N} \text{ and } \tilde{x} \in [-1, 1]. \qquad (9) $$

It can be noted here that the Chebyshev roots can be easily evaluated through the following generic formula:

$$ \tilde{x}_i = \cos \left( \frac{2i - 1}{2n} \pi \right), i = 1, \dots, \text{card}(I_n). \qquad (10) $$

Now, if one chooses a set of $(2N+1)$-points, $\{\tilde{x}_j\}_{j \in D_{2N+1}}$, as the Chebyshev roots that are included in $[-1, 1]$, (1) becomes

$$ f(\tilde{x}) = \sum_{j \in D_{2N+1}} P_j(\tilde{x}) f(\tilde{x}_j) + \frac{\partial^2 T_n(\tilde{x})}{2^2 N (2N+1)!} T_{2N+1}(\tilde{x}), N \in \mathbb{N} \text{ and } \tilde{x} \in [-1, 1]. \qquad (11) $$

Because the derivative of $T_n(\tilde{x})$, $n \in \mathbb{N}$, is given by

$$ \partial_{\tilde{x}}[T_n](\tilde{x}) = nU_{n-1}(\tilde{x}), n \in \mathbb{N} \text{ and } \tilde{x} \in [-1, 1], \qquad (12) $$
where $U_n(\tilde{x})$ is the Chebyshev polynomial of second kind,

$$ U_n(\tilde{x}) = \frac{\sin((n+1)\cos^{-1}(\tilde{x}))}{\sin(\cos^{-1}(\tilde{x}))}, n \in \mathbb{N} \text{ and } \tilde{x} \in [-1, 1], \quad (13) $$

the differentiation of (11) at the central point of the Chebyshev stencil, $\tilde{x}_0 = 0$, leads to

$$ \partial_{\tilde{x}}[f](\tilde{x}_0) = \sum_{j \in D_{2N+1}} \partial_{\tilde{x}}[P_j](\tilde{x}_0) f(\tilde{x}_j) \pm \frac{\partial^{2N+1}_{\tilde{x}}[f](\tilde{\xi})}{2^{2N}(2N)!}, N \in \mathbb{N}. \quad (14) $$

It has to be noticed that, because the derivative is evaluated at point $\tilde{x}_0 = 0$, $T_m(0) = 0$ (resp. $U_m(0) = \pm 1$), for $m$ odd (resp. $m$ even). More importantly, as can be seen, the error term decreases as $2^{2N}$.

From (14), and analogously to the reasoning made in Section 3.1, one can thus conclude that the leading-order truncation error term is now directly proportional to $E_{cheb}$, defined by

$$ E_{cheb} = \frac{1}{2^{2N}(2N)!}, N \in \mathbb{N}. \quad (15) $$

The key point is that such an error is always going to be smaller than the one associated with a standard FD scheme of same stencil size, which is proportional to $E_{std}$; see (7).

The optimization proposed here thus consists of defining FD schemes as being based on specific points given by the Chebyshev nodes, rather than on equidistant ones (as is implicitly done for both standard and spectral-like schemes; see (A.1)). In all of the following, such an optimization approach will be referred to as the intrinsic optimization.

### 3.3. Intrinsic optimization: Practical implementation

From a practical point of view, numerically differentiating a discrete signal, $f(.)$, with the help of an FD scheme first requires to pass from the physical space (where the function is primarily defined) to the numerical one (where its discrete derivative can be properly approximated), following relations (5) and (6).

Here, we can recall that any FD scheme is implicitly associated with a particular set of discrete points whose spatial distribution determines the respective values of the scheme's coefficients $a_j$. For instance, as it is detailed in Appendix A, both standard and spectral-like optimized FD schemes are based on equally-spaced points of unitary spacing (see Section 3.1.2). On the other hand, it was mentioned earlier that intrinsically-optimized FD (IOFD) schemes rely on Chebyshev roots whose spatial distribution is even more specific (10).

However, and whatever its nature is, each particular set of discrete points has no reason to match the one that is naturally inherited from the discretization of the problem to be solved, with the exception of those trivial situations where both sets may collapse, for example, when a standard FD scheme is used in conjunction with a Cartesian mesh grid.

Nevertheless, in practice, one may pass from the discretization (or physical) space, $(x, y, z)$, to the solving (or numerical) one, $(\tilde{x}, \tilde{y}, \tilde{z})$, by making use of a time invariant coordinate transformation $\mathbf{J}: x, y, z \mapsto \tilde{x}, \tilde{y}, \tilde{z}$ (also referred to as Jacobian transformation), which maps the former onto the latter.

$$ \mathbf{J} = \frac{\partial(\tilde{x}, \tilde{y}, \tilde{z})}{\partial(x, y, z)}. \quad (16) $$

In practice, the so-called intrinsic optimization consists thus of virtually evaluating the FD-scheme-based discrete derivative onto an orthogonal set of Chebyshev nodes $(\tilde{x}, \tilde{y}, \tilde{z})$, once the latter has been mapped with the discretization space, $(x, y, z)$, via a Jacobian transformation (see details in Appendix B). The differential operator then becomes
$$ \nabla_{(x,y,z)}[f](x_0, y_0, z_0) \cong \sum_{j \in D_{2N+1}} \mathbf{J}(\partial_{\tilde{x}}[\tilde{P}_j](\tilde{x}_0) f(\tilde{x}_j) \mathbf{e}_{\tilde{x}} + \partial_{\tilde{y}}[\tilde{P}_j](\tilde{y}_0) f(\tilde{y}_j) \mathbf{e}_{\tilde{y}} + \partial_{\tilde{z}}[\tilde{P}_j](\tilde{z}_0) f(\tilde{z}_j) \mathbf{e}_{\tilde{z}}). \quad (17) $$

Regarding more especially centered FD schemes, the intrinsic optimization procedure can be sketched as a succession of 4 steps, which are detailed in Appendix B, and summarized hereafter (only with respect to the direction $\tilde{x}$, for the sake of clarity).

1. Evaluation of the FD coefficients $a_j^c = \partial_{\tilde{x}}[\tilde{P}_j](x_0)$, $j \in D_{2N+1}$, that are associated with the discrete operator used in the Chebyshev-based numerical space, i.e. the one corresponding to the Chebyshev nodes, $\{\tilde{x}_j\}_{j \in D_{2N+1}}$, included in $[-1, 1]$. This is achieved in a straightforward manner, via analytical means.

2. Evaluation of the Jacobian transformation that maps the discretization points of the physical space $(x, y, z)$ onto the Chebyshev-based nodes of the numerical space, $(\tilde{x}, \tilde{y}, \tilde{z})$. This is achieved by manipulating the inverse Jacobian transformation $\mathbf{J}^{-1}$.

3. Evaluation of the values taken by the function $f(\cdot)$ on the Chebyshev nodes, $f(\tilde{x}_j)$. This step is required because the Jacobian transformation associated with the Chebyshev nodes is not bijective (see detail in Appendix B). Here, it is achieved by interpolating the function known values $f(x_j)$, such an interpolation step being conducted thanks to centered $n^{th}$-order Lagrange polynomial interpolations (which offer the advantage of being amplification-free [21]).

4. Evaluation of the approximate derivative on the Chebyshev nodes, which is achieved in a straightforward manner, following (17).

One can illustrate more easily the intrinsic optimization technique by considering a one-dimensional grid of equally-spaced points, and a virtual mesh constituted of their associated $\mathbf{J}^{-1}$-transformed Chebyshev points, as described in Figure 1 for the particular configuration $N = 4$ and $n = 6$. Here, differentiating a given discretized function in an intrinsically-optimized sense corresponds thus to (i) interpolating it with the help of a centered sixth-order Lagrange interpolation, before (ii) projecting the interpolated (continuous) result onto the transformed space of $2 \times 4 + 1 = 9$ Chebyshev points, and then (iii) virtually differentiating it, using the Chebyshev coefficients (see Table B.1, in Appendix B).

<figure><img src="image_3.png" /><figcaption>Figure 1. Noncoincidence of the discretization $(x)$ and of the $\mathbf{J}^{-1}$-transformed Chebyshev points, except at the origin $(O)$, where the derivative is to be calculated. The particular case of $N = 4$ (i.e. $2N + 1 = 9$) and $n = 6$ is illustrated, where the unknown values taken by the function at the Chebyshev nodes are to be obtained via an interpolation (the update of point $\mathbf{J}^{-1}\tilde{x}_2$ is illustrated here).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. The same as Figure 1 with, this time, the update of the point $\mathbf{J}^{-1}\tilde{x}_3$ illustrated. In the present case, a total of 11 donor points are required to properly evaluate the approximate derivative.</figcaption></figure>
In Figures 1 and 2, it can be seen how the interpolation step is handled, depending on the Chebyshev points location. In particular, in the present case, making use of a centered sixth-order Lagrange interpolation implies that the grid point $b + \Delta x$ (resp. $a - \Delta x$) value is also used for feeding the points $J^{-1}(\tilde{x}_3)$ and $J^{-1}(\tilde{x}_4)$ (resp. $J^{-1}(\tilde{x}_{-3})$ and $J^{-1}(\tilde{x}_{-4})$).

## 4. INTRINSIC OPTIMIZATION TECHNIQUE: SPECIFIC CONSIDERATIONS AND PRACTICAL ASPECTS

In view of what precedes, it can be concluded that any FD scheme to be implicitly associated with Chebyshev nodes shall constitute the ideally optimized differential operator, for the sole reason that its Chebyshev-based numerical (or solving) space is optimally designed.

However, as explained earlier and detailed in Appendix B, because the Jacobian transformation associated with the Chebyshev nodes is not bijective, the intrinsic optimization technique generally requires that an interpolation step is performed. One can thus wonder about how such an interpolation operation may impact the overall optimization procedure, whatever it is by questioning the effective quality (accuracy, robustness, etc.) of the resulting schemes [22, 23], or by limiting their practical use (numerical implementation, etc.). This will be examined in the next subsections.

### 4.1. Impact of the interpolation step on the accuracy of the intrinsically-optimized schemes

The impact of the interpolation stage on the overall accuracy of a given IOFD scheme is highlighted in Figure 3, with respect to the Lagrange interpolation basis used here. Firstly, Figure 3(a) compares the dispersion errors associated with both the centered standard 9-point 8-order scheme (FDs9p8o) and its intrinsically-optimized counterpart (FDc9p8o). The latter scheme is considered here as ideally optimized, that is, intrinsically-optimized without any interpolation applied. Its dispersion error can be accurately estimated with the help of a harmonic function, $f(x) = e^{i\alpha x}$, and an adapted version of (B.5), described in Appendix B:

$$ \partial_x[f](x_0) \approx J \sum_{j \in D_{2N+1}} a_j^c f(\mathbf{J}^{-1}(\tilde{x}_j)). \quad (18) $$

Given that, in this case, $f(.)$ is analytically known, no interpolation is needed.

As can be seen, in such an interpolation-free situation, the dispersion properties of the IOFD scheme are indubitably better than the ones offered by its standard counterpart. Nevertheless, Figure 3(b) shows that these good properties are partly degraded from the moment the interpolation stage is to be applied. In the present case, use was made of centered Lagrange interpolations of various stencil width ($n \in \{2, 4, 6, \dots\}$). Knowing that the stencil width drives the interpolation order-of-accuracy (which, in the present centered case, equals $n$), one can appreciate how improving the interpolation step helps in partially remedying the quality loss it induces.

<figure><img src="image_5.png" /><figcaption>Figure 3. Left side (a); dispersion error components for the centered 9-point stencil finite-difference scheme based on both equally-spaced (—FDs9p8o) and Chebyshev (•FDc9p8o) points. Right side (b); corresponding intrinsically-optimized schemes, for several Lagrange interpolation schemes of increased accuracy order (in bold dashed-lines, the 2$^{nd}$- and the 26$^{th}$-order accurate ones).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4. Dispersion error components for the FDs9p8o, FDs13p12o and the optimized scheme obtained via the use of a 8-knot cubic-spline (FDsp13p8k) instead of the Lagrange interpolation. —FDs9p8o, •FDs13p12o and — FDsp13p8k.</figcaption></figure>

Here, one can notice that, whatever the IOFD scheme is, its associated approximated wavenumber, $\tilde{\alpha}_1$, is purely real. This indicates that their intrinsic-optimization (and underlying interpolation step) does not question the amplification-free character of centered FD schemes, which remain solely characterized by dispersion effects (see Subsection A.3). This was indeed to be expected, since the intrinsic optimization procedure solely relies on amplification-free stages: FD approximations based on centered antisymmetric stencils, and Lagrange interpolations based on centered stencils composed of an even number of points.

At this point, it can be seen that the interpolation stage could even rely on an interpolation basis other than the centered Lagrange ones that were privileged here (e.g. splines, see Figure 4). Studying the impact of different interpolation bases on the overall efficiency of the intrinsic optimization process would constitute an interesting subject of research, and will be the focus of future works.

## 4.2. Equivalent intrinsically-optimized finite-difference schemes

In the present subsection, we discuss how one can easily and advantageously turn any standard FD scheme into a so-called equivalent intrinsically-optimized one, so that the latter can be numerically implemented in a straightforward manner, as well as compared against the standard ones on a same “accuracy order versus CPU cost” basis.

Indeed, applying as is the intrinsic optimization technique requires to first implement and then compute: (i) the interpolation (to occur from the physical towards the numerical spaces), and (ii) the derivation (to be achieved in the numerical space). Needless to say, this is likely to increase both the complexity of the effective implementation of IOFD schemes and the computational price induced by their use.

Therefore, in practice, it is much more advantageous to make use of the so-called equivalent IOFD schemes, which can be easily derived and implemented, and whose use does not imply any extra CPU cost. Indeed, any of these equivalent IOFD schemes corresponds in all points to its (standard or spectral-like optimized) FD counterpart of same stencil size, with the exception of its set of coefficient values. For deriving such an equivalent IOFD scheme, it is only required to consider a given couple of derivation and interpolation stencils (of $2N + 1$ and $n$ points, respectively), and to rearrange the right-hand side of (B.5), Appendix B, in the following equivalent form:

$$J \sum_{j \in D_{2N+1}} a_j^c \mathcal{I}[f](\mathbf{J}^{-1}\tilde{x}_j) = \frac{1}{\Delta x} \sum_{j \in D_{2N^{eq}+1}} a_j^{eq} f(x_j), \quad (19)$$

where $N^{eq} \in \mathbb{N}$ depends on both $(n, N)$, and verifies $N^{eq} \ge N$.
<table><caption>Table II. Two particular sets of (N, n), and their corresponding N<sup>eq</sup>-point equivalent intrinsically-optimized schemes.</caption><thead><tr><th></th><th>Chebyshev Stencil (2N + 1)</th><th>Interpolation Stencil (n)</th></tr></thead><tbody><tr><td>FDeq11p6o</td><td>9</td><td>6</td></tr><tr><td>FDeq15p8o</td><td>11</td><td>8</td></tr></tbody></table>

For instance, when considering both a 9-point (N = 4) derivation and a 6-point (n = 6) interpolation stencil, one ends up with a 11-point equivalent stencil (N<sup>eq</sup> = 5)†, as illustrated in Figures 1 and 2. Because of the fact that such a scheme naturally inherits the smaller order-of-accuracy between O(Δx<sup>2N</sup>) and O(Δx<sup>n</sup>), one can expect its effective accuracy to go along with the polynomial level of approximation used for the interpolation stage.

However, increasing the interpolation accuracy (and thus n) requires enhancing the equivalent stencil width (i.e., N<sup>eq</sup>). Hence, one has to make sure that the numerical efficiency of the resulting scheme is not degraded when compared with its standard or spectral-like counterpart of identical stencil width (and thus of the same CPU price).

Therefore, it is preferable to choose the pair (N, n) ensuring that the resulting IOFD scheme has an optimal quality/price ratio. As an illustration, for deriving an equivalent intrinsically-optimized scheme that is more particularly suited for use within a CAA context, it is preferable to define the optimal pair (N, n) as that offering the best dispersion properties, while ensuring additional constraints such that

1. a minimal order-of-accuracy of 6 is respected for the error made on low-frequency contents remains minimal while a fast error level reduction is ensured when the mesh grid resolution increases; and

2. the resulting equivalent intrinsically-optimized scheme, FDeq(2N<sup>eq</sup>+1)p(n)o, is more accurate than its standard counterparts of same stencil width, FDs(2N<sup>eq</sup> + 1)p(2N<sup>eq</sup>)o, and of same formal accuracy order, FDs(n + 1)p(n)o, so that the low-frequency to high-frequency contents are accurately captured.

As an example, choosing the pair N = 4 and n = 6 leads to an equivalent intrinsically-optimized scheme (FDeq11p6o) that verifies the previous constraints, when compared with its standard counterparts of same stencil width (FDs11p10o) or order-of-accuracy (FDs7p6o). See sections 4.3 and 5 for more details.

Table II displays two particular couples (N, n) that were considered in the present work. It should be noted that only cases with n ≥ 6 were retained, so that the resulting order-of-accuracy of equivalent intrinsically-optimized schemes equals the interpolation one (n). Table III delivers the respective coefficient sets associated with each one of these two particular equivalent IOFD schemes, whose values were derived in a straightforward manner via the algorithm provided in Appendix C. Please note that such an algorithm may easily be implemented, allowing then to derive at no expense the set of coefficients associated with any IOFD scheme.

### 4.3. Spectral analysis of two equivalent intrinsically-optimized schemes

In the present Section, the two particular equivalent IOFD schemes introduced earlier are characterized from a spectral point of view, each one being then compared with its standard counterpart.

Figure 5 depicts the dispersion errors (A.13) induced by these two pairs of schemes when applied to a harmonic signal such as the one introduced in section A.3. As can be seen, despite exhibiting a lower order-of-accuracy, the two equivalent intrinsically-optimized schemes (FDeq11p6o and FDeq15p8o) offer better dispersion properties than their respective standard counterparts of same stencil width, FDs11p10o and FDs15p14o.

<footnote>†In all of the following, centered (2N<sup>eq</sup> + 1)-point n<sup>th</sup>-order equivalent intrinsically-optimized schemes are referred to as FDeq(2N<sup>eq</sup> + 1)p(n)o</footnote>
<table><caption>Table III. Stencil coefficients associated with the equivalent intrinsically-optimized finite-difference schemes FDeq11p6o and FDeq15p8o.</caption><thead><tr><th></th><th>FDeq11p6o</th><th>FDeq15p8o</th></tr></thead><tbody><tr><td>a<sub>0</sub><sup>eq</sup></td><td>0</td><td>0</td></tr><tr><td>a<sub>1</sub><sup>eq</sup> = -a<sub>-1</sub><sup>eq</sup></td><td>0.865816409136225</td><td>0.899944730463117</td></tr><tr><td>a<sub>2</sub><sup>eq</sup> = -a<sub>-2</sub><sup>eq</sup></td><td>-0.275690008254363</td><td>-0.326219480460150</td></tr><tr><td>a<sub>3</sub><sup>eq</sup> = -a<sub>-3</sub><sup>eq</sup></td><td>0.0811126692593992</td><td>0.124540388374523</td></tr><tr><td>a<sub>4</sub><sup>eq</sup> = -a<sub>-4</sub><sup>eq</sup></td><td>-0.0165005999177029</td><td>-0.0404962578862881</td></tr><tr><td>a<sub>5</sub><sup>eq</sup> = -a<sub>-5</sub><sup>eq</sup></td><td>0.00164559985302308</td><td>0.00985168108601892</td></tr><tr><td>a<sub>6</sub><sup>eq</sup> = -a<sub>-6</sub><sup>eq</sup></td><td></td><td>-0.00152738609091035</td></tr><tr><td>a<sub>7</sub><sup>eq</sup> = -a<sub>-7</sub><sup>eq</sup></td><td></td><td>0.000109143999162341</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 5. Dispersion error associated with various finite-difference schemes: —FDs11p10o; — FDeq11p6o; •FDs15p14o; and ··· FDeq15p8o.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. (a) Relative dispersion and (b) group velocity errors associated with various finite-difference schemes: —FDs7p6o; •FDs11p10o; — FDeq11p6o.</figcaption></figure>

With the purpose of better assessing the advantages offered by these two equivalent IOFD schemes, in the following, their effective gain is further analyzed.

### 4.3.1. Intrinsically-optimized FDeq11p6o scheme.
Firstly, the equivalent intrinsically-optimized FDeq11p6o scheme is compared against its standard counterparts of either (i) same order-of-accuracy, FDs7p6o, or (ii) same stencil width, that is, nominal CPU cost, FDs11p10o - such comparison being conducted in respect of either the relative dispersion (Figure 6(a)) or the group velocity (Figure 6(b)) errors.
<table><caption>Table IV. Minimal number of points per wavelength (PPW) and associated numerical efficiency (as $(2N^{(eq)} + 1) \times \text{PPW}$) for a global accuracy error inferior to 1%.</caption><thead><tr><th></th><th>PPW<sub>d</sub></th><th>(2N<sup>(eq)</sup> + 1)× PPW<sub>d</sub></th><th>PPW<sub>g</sub></th><th>(2N<sup>(eq)</sup> + 1)× PPW<sub>g</sub></th></tr></thead><tbody><tr><td>FDs7p6o</td><td>5.66</td><td>39.6</td><td>7.85</td><td>55.0</td></tr><tr><td>FDs11p10o</td><td>4.13</td><td>45.4</td><td>5.32</td><td>58.5</td></tr><tr><td>FDeq11p6o</td><td>3.51</td><td>38.6</td><td>4.18</td><td>46.1</td></tr></tbody></table>

PPW, points per wavelength.

<figure><img src="image_9.png" /><figcaption>Figure 7. (a) Dispersion and (b) group velocity errors associated with various finite-difference schemes: —FDs9p8o; •FDs15p14o; — FDeq15p8o.</figcaption></figure>

As can be seen, the maximum accuracy gain§ offered by the *FDeq11p6o* scheme corresponds to an accuracy limit of 1%, in terms of group velocity error. And, when compared with its standard counterparts, *FDs7p6o* and *FDs11p10o*, such *FDeq11p6o* scheme requires fewer points-per-wavelength (PPW) for ensuring a global accuracy of less than 1%, in terms of both dispersion (PPW<sub>d</sub>) and group velocity (PPW<sub>g</sub>) errors.

More precisely, Table IV summarizes the minimal PPW needed for all three schemes to ensure a global error of less than 1%. From that number, one can first conclude that, compared with its standard counterpart of same accuracy order (*FDs7p6o*), the *FDeq11p6o* scheme allows handling grids almost twice as coarse (7.85 PPW<sub>g</sub> against 4.18 PPW<sub>g</sub>, respectively), while ensuring a better numerical efficiency (55.0 against 46.1, respectively)¶.

This conclusion is important because it indicates that using the latter IOFD scheme instead of the (very popular) *FDs7p6o* one allows dealing with grids that are twice as coarse as in 1D, that is, lighter by a factor of 8 in 3D. This, knowing that the latter gain is likely to be multiplied by 2 (leading to a factor 16), from the moment the time step to be used in the simulation is solely driven by CFL criteria, and thus, by the mesh discretization level.

Now, when comparing it with its standard counterpart of same stencil width (*FDs11p10o*), it can be seen that the *FDeq11p6o* scheme allows dealing with grids up to 21% coarser (5.32 PPW<sub>g</sub> and 4.18 PPW<sub>g</sub>, respectively), while leading to a still higher relative numerical efficiency.

### 4.3.2. Intrinsically-optimized FDeq15p8o scheme.
As previously carried out for the equivalent intrinsically-optimized *FDeq11p6o*, the *FDeq15p8o* scheme is now compared with its standard counterparts of (i) same order-of-accuracy, *FDs9p8o*, and of (ii) same stencil width, *FDs15p14o*. Again, such comparison is conducted in respect of either the relative dispersion (Figure 7(a)) or the group velocity (Figure 7(b)) errors.

<footnote>§The accuracy gain between the considered IOFD scheme and its standard counterparts (of same stencil width and accuracy order), such gain being defined by the reduction of group velocity error.</footnote>

<footnote>¶The numerical efficiency is defined here as the product of PPW and the number of stencil points [11].</footnote>
<table><caption>Table V. Minimal number of points per wavelength (PPW) and associated numerical efficiency (as $(2N^{(eq)} + 1) \times \text{PPW}$), for a global accuracy error inferior to 0.14%.</caption><thead><tr><th></th><th>PPW<sub>d</sub></th><th>(2N<sup>(eq)</sup> + 1)× PPW<sub>d</sub></th><th>PPW<sub>g</sub></th><th>(2N<sup>eq</sup> + 1)× PPW<sub>g</sub></th></tr></thead><tbody><tr><td>FDs7p6o</td><td>8.08</td><td>56.6</td><td>11.22</td><td>78.5</td></tr><tr><td>FDs9p8o</td><td>6.16</td><td>55.4</td><td>8.16</td><td>73.4</td></tr><tr><td>FDs15p14o</td><td>4.26</td><td>63.9</td><td>5.28</td><td>79.2</td></tr><tr><td>FDeq15p8o</td><td>3.57</td><td>53.6</td><td>3.98</td><td>59.7</td></tr></tbody></table>

PPW, points per wavelength.

As revealed by a closer look at Figure 7(b), the maximum accuracy gain<fnref n="1" /> offered by the *FDeq15p8o* scheme is reached for an accuracy limit of 0.14%, in terms of group velocity error, whose associated minimal PPW is provided in Table V. Here again, when comparing it with its standard counterpart of same order-of-accuracy (*FDs9p8o*), it can be seen that such a *FDeq15p8o* scheme allows handling meshes of twice the coarseness (8.16 PPWg against 3.96 PPWg, respectively), while offering a better numerical efficiency (73.4 against 59.7, respectively).

On another hand, when compared with its standard counterpart of same stencil width (*FDs15p14o*), this *FDeq15p8o* scheme allows dealing with grids up to 25% coarser (5.28 PPWg against 3.96 PPWg, respectively), while offering a still higher numerical efficiency.

Finally, it can be seen here that, with no more than 3.98 PPWg, such an equivalent IOFD scheme of 15 points guarantees that the error made on the group velocity (on which depends the acoustic energy transport and, thus, the overall accuracy of the CAA stage) is less than 0.14%. Compared with the standard 7-point 6-order FD scheme (*FDs7p6o*), whose same error level is reached only for a 11.22 PPWg value, the present scheme, *FDeq15p8o*, thus offers a gain of roughly 3, in terms of grid density per direction, that is, a factor of 27 in 3D.

Again, this conclusion is of importance, because it indicates that using such an IOFD scheme, rather than the very popular *FDs7p6o* one, allows dealing with grids that are three times coarser, that is, lighter by a factor of 27 in 3D. Here again, this can be immediately translated into a reduction (by the same factor) of the memory needed. This, knowing that the latter gain is likely to be multiplied by 3 (leading to a factor 81), from the moment the time step used in the simulation is solely driven by CFL criteria, that is, by the grid discretization level. On another hand, the sole penalty to be possibly induced by such a scheme, *FDeq15p8o*, is the one to be logically expected from its larger stencil width (e.g., a possibly higher computational time because of parallel computing transmissions if the schemes are kept, as is the case for ensuring the interblock or intercore connections, etc.). Such an extra cost, however, is expected to be rather marginal, thanks to the very good numerical efficiency offered by the present *FDeq15p8o* scheme (59.7, as presented in Table V).

Finally, by cross-comparing Tables IV and V, it can be observed that, for nearly the same computational cost as the *FDs7p6o* scheme (59.7 against 55.0, respectively), the *FDeq15p8o* one leads to more accurate results, while requiring half the number of points per wavelength (7.85 PPWg against 3.98 PPWg, respectively).

## 5. ILLUSTRATION AND VALIDATION OF THE EQUIVALENT INTRINSICALLY-OPTIMIZED SCHEMES VIA AN APPLICATION TO 1D AND 2D ACADEMIC TEST CASES

In the present section, the two particular equivalent IOFD schemes introduced earlier are further assessed before they are compared with both their standard and their spectral-like optimized counterparts. Such assessment, which is conducted on the basis of two academic test cases of one and two dimensions, respectively, allows illustrating and validating further the advantages effectively offered by the intrinsic optimization technique.

<footnote n="1">¹Again, the accuracy gain between the considered IOFD scheme and its standard counterparts (of same stencil width and accuracy order), such gain being defined by the reduction in group velocity error.</footnote>
## 5.1. Advection of a 1D harmonic perturbation

With a view to illustrating further the capacities of IOFD schemes, we consider here the initial value problem given by the one-dimensional convective equation (20) whose solution consists of an initial disturbance that travels along the x-axis at a unitary dimensionless speed. The dimensionless variables are $\Delta x$ as the length scale and $\Delta x/c$ as the time scale, with c being the speed of sound.

$$ \partial_t u + \partial_x u = 0. \qquad (20) $$

The initial disturbance, which definition follows [11], is given (at $t = 0$) by

$$ u(x, t = 0) = Asin \left( \frac{2\pi}{a\Delta x} (x - x_0) \right) exp \left( -\ln(2) \left( \frac{(x - x_0)}{18\Delta x} \right)^2 \right), \qquad (21) $$

where A denotes the disturbance amplitude and $a\Delta x$ its wavelength. The latter, which equals the number of points-per-wavelength, PPW, is here successively taken as $a\Delta x = 4, 8$ or $16$.

The evolution problem associated with (20) is numerically handled via a six-stage Runge-Kutta time marching scheme. The latter is allotted a time step of $\Delta t = 0.10$, which provides a CFL number of 0.10. Such a low CFL value is here privileged so that the (unavoidable) approximation errors coming from the time marching scheme induce only negligible bias. Additionally, in order to keep the results free from any other side effects, the simulation is conducted without any filtering procedure; only the impact by the discrete differential operators are thus assessed.

Figure 8 plots the $L^2$-norm error (22) that is recorded for the various schemes, after the disturbance had been advected over $50a\Delta x$.

<figure><img src="image_10.png" /><figcaption>Figure 8. $L^2$-norm error made on the disturbance by various finite-difference schemes, for a propagation distance of $50a\Delta x$: (a) $+FDs7p6o; \circ FDs11p10o; *FDeq11p6o; \square FD011p4o$ and (b) $+FDs9p8o; \circ FDs15p14o; *FDeq15p8o; \square FD015p4o$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 9. (a) Fourier Transform (FT) of the initial disturbance and (b) its associated content, as numerically obtained by various finite-difference schemes after a numerical propagation over a distance of $50a\Delta x$, with 4 points per wavelength. $\circ$ Analytic solution; $\cdots FDs9p8o; - FDs15p14o; -FDeq15p8o$.</figcaption></figure>
$$E_{num} = \sqrt{\frac{\sum_{n}(u_{exact} - u_{num})^2}{\sum_{n}(u_{exact})^2}}. \qquad (22)$$

From Figures 8(a) and 8(b), one can draw several important conclusions.

1. Compared with the well-known standard seven-point scheme (FDs7p6o), and whatever the frequency range is (4, 8, or 16 PPW), the equivalent IOFD scheme FDeq15p8o delivers better results, with half the number of PPW. This confirms the conclusion previously made, according to which the latter scheme is more accurate than the former, although it induces nearly the same computational cost, and half the memory.

2. In addition to that, over the high-frequency range (4 PPW), an equivalent intrinsically-optimized scheme such as FDeq15p8o one outperforms its standard counterparts of same stencil width and accuracy order (FDs15p14o and FDs9p8o, respectively). Indeed, as illustrated in Figure 9 and quantified in Table VI, of all three schemes, FDeq15p8o is the sole one to accurately solve the signal with no more than 4 PPW.

3. Regarding now its spectral-like optimized counterpart of same stencil width (FDol11p4o), and as can be inferred from Figure 8(a), the equivalent IOFD scheme FDeq11p6o presents about the same accuracy over the mid-frequency to high-frequency ranges (4 and 8 PPW), while it turns out to be more accurate for signals of lower frequency (e.g. 16 PPW). In that last case, the spectral-like optimized FD scheme (FDol11p4o) presents a higher error level than the IOFD one (FDeq11p6o).

4. Finally, when the mesh refinement level increases, both spectral-like optimized schemes, FDol11p4o and FDol5p4o, do not see their error level decreasing substantially. On the contrary, both their standard and equivalent intrinsically-optimized counterparts of same stencil width adopt a much better behavior, because of their higher order-of-accuracy.

From what precedes, one can conclude that IOFD schemes constitute a good alternative to both standard and spectral-like optimized schemes, offering a very favorable tradeoff in terms of accuracy versus cost ratio.

In particular, the IOFD, FDeq11p6o, offers more advantages than its standard counterpart of same order (resp. stencil width) FDs7p6o (resp. FDs11p10o). Indeed, compared with the latter scheme, the former one ensures a higher accuracy on coarser grids, as well as a greater numerical efficiency. However, over the mid-frequency to high-frequency ranges, such a FDeq11p6o scheme behaves as well as its spectral-like optimized counterpart of same stencil width (FDol11p4o), outperforming it over the low-frequency range. This is, in particular, what makes such an IOFD scheme more suitable for application to problems that may involve various grid refinement levels. As a matter of fact, such a scheme (i) offers a high accuracy (that equals the spectral-like optimized one) when the discretization level is marginal, while also (ii) ensuring a fast error level reduction (that equals the standard one) when the mesh grid resolution increases.

Regarding the FDeq15p8o IOFD, it can be more advantageously used whenever still more accurate results are desired. In particular, such a scheme is more especially indicated whenever a high accuracy must be achieved for a broadband range of frequencies, especially if the acceptable error level is less than 1%.

<table><caption>Table VI. Convection problem: $L^2$-norm error (Equation (22)) made by various finite-difference schemes, for a propagation distance of $50a\Delta x$, with 4 points per wavelength.</caption><thead><tr><th></th><th>FDs9p8o</th><th>FDs15p14o</th><th>FDeq15p8o</th></tr></thead><tbody><tr><th>E<sub>num</sub></th><td>1.4941</td><td>0.9885</td><td>0.1118</td></tr></tbody></table>
## 5.2. Advection of a 2D Vortex

Let us now replace the initial disturbance by the vortex [2] defined in (23), and the previous one-dimensional homogeneous grid by a 2D Cartesian one with $\Delta x = \Delta y = 1$.

$$ \left\{ \begin{array}{l} u = 1 - \frac{C(y-y_c)}{R^2} \exp \left( -\frac{r^2}{2} \right), \\ v = \frac{C(x-x_c)}{R^2} \exp \left( -\frac{r^2}{2} \right), \\ r^2 = \frac{(x-x_c)^2+(y-y_c)^2}{R^2}. \end{array} \right. \qquad (23) $$

In the aforementioned equation, $u, v$, and $R$ denote the two Cartesian velocity components and the vortex core radius, respectively. The Mach number is given by $M_\infty = 0.1$, and the parameter $C = 0.02 M_\infty R$.

After the vortex has been numerically advected over a distance equaling 100 times its core radius, the $L^2$-norm errors (22) associated with the various schemes are determined. Those errors are here evaluated for both the $u$-velocity and the $v$-velocity components at the position $y = y_c + R$ (Figure 10) and $y = y_c$ (Figure 11), respectively, corresponding thus to the maximal error made along a given grid line. From those figures, it can be clearly seen that IOFD schemes are extensively more accurate than spectral-like optimized ones, when it comes to solving functions whose support has predominantly low-frequency components that do not coincide exactly with the spectral-like optimization interval. In addition to speaking in favor of IOFD schemes, this observation confirms how important the effective order-of-accuracy can be, when it comes to numerically simulating phenomena that are not purely monochromatic.

<figure><img src="image_12.png" /><figcaption>Figure 10. Maximum $L^2$-norm error made on the $u$-velocity component at $y = y_c + R$ by various finite-difference schemes, after a convection distance of 100R. (a) +FDs7p6o; ◦FDs11p10o; *FDeq11p6o; □FDo11p4o and (b) +FDs9p8o; ◦FDs15p14o; *FDeq15p8o; □FDo15p4o.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Maximum $L^2$-norm error made on the $v$-velocity component at $y = y_c$ by various finite-difference schemes, after a convection distance of 100R. (a) +FDs7p6o; ◦FDs11p10o; *FDeq11p6o; □FDo11p4o and (b) +FDs9p8o; ◦FDs15p14o; *FDeq15p8o; □FDo15p4o.</figcaption></figure>
# 6. FURTHER ILLUSTRATION AND VALIDATION OF THE EQUIVALENT INTRINSICALLY-OPTIMIZED SCHEMES VIA AN APPLICATION TO 3D PROBLEM

With the purpose of illustrating and validating further the advantages offered by equivalent IOFD schemes in a more applied context, the present section addresses the problem of acoustic scattering by a 3D sphere within a quiescent medium.

## 6.1. Test case definition and generic calculation

We consider a sphere of unitary radius located at the origin $\mathbf{x} = (0, 0, 0)$ of a quiescent medium. Such sphere acts as a scattering agent onto a quadrupolar acoustic source that is located $\mathbf{x} = (-1.5, 0, 0)$ away from the sphere's center, and that pulsates with a dimensionless frequency $k = 2\pi$ (the spatial scale is taken as $R$, and the temporal one as $R/c$, with $R$ the unitary radius).

Please note that such a test case is directly inspired by the exercise validation proposed in the Benchmark on Airframe Noise Computations (BANC) initiative managed by NASA/LaRC. This problem, which corresponds to the BANC-II/Category 8/Test Case 1, was jointly designed by NASA, JAXA, and Onera to serve as a generic problem for assessing the various ways of far-field extrapolating a set of near-field unsteady data. Indeed, for this problem, there exists an analytical solution (which was provided here by Dr Ikeda, from JAXA) from which one can force any acoustic propagation tool. For doing so, the area encompassing both the sphere and the source is considered as constituting the near-field region, which is then acoustically extrapolated via a CAA calculation, thereby leading to an analytical-CAA weakly coupled calculation.

Such an analytical-CAA coupled calculation exercise had been achieved and validated in a previous effort [24–26]; for such a calculation, and although this could have been done (such as was carried out in another framework [27]), the sphere obstacle was not CAA-meshed, so that the propagation grid could be taken as an ideal one, that is, of a homogenous Cartesian type. The reason for doing so was to minimize at best all possible side effects (such as metric errors) that could have biased the analysis. Following the BANC guidelines, the near-field region (encompassing both the sphere and the quadrupole area) was considered as constituting the source area ($\Omega_{source} = [-2, 1.5] \times [-1.5, 1.5] \times [-1.5, 1.5]$). The latter was then CAA-extrapolated, the CAA perturbed field being implicitly forced at each time step via the so-called Non Reflective Interface (NRI) [24–26] technique, which was fed with the analytical dataset. The CAA stage calculation was run for six source periods ($T = 1$) with 20 iterations per period ($CFL = 1/2$). This CAA computation was achieved with the help of a centered standard 7-point FD derivative scheme (FDs7p6o), being performed over one core of a laptop (64 bits, Pentium 2.1 GHz, 4 Gb Ram), which required approximately 30 min of wall clock time. The resulting instantaneous perturbed field (Figure 12)

<figure><img src="image_14.png" /><figcaption>Figure 12. Analytical-CAA weakly coupled calculation. Instantaneous perturbed pressure at $t = 6T$.</figcaption></figure>
was very favorably compared with the analytical results (see [24–26]), thus validating both the Analytical-CAA extrapolation technique and the CAA numerical setup.

## 6.2. Alternative calculations based on equivalent IOFD schemes

With a view to assessing further the equivalent IOFD schemes, the previous generic Analytical-CAA calculation is repeated here, its computational setup being, however, modified in several aspects.

<figure><img src="image_15.png" /><figcaption>Figure 13. Comparison between different finite-difference schemes for the sphere test case with $k = 5\pi$. In solid lines, the analytical result and, in dashed ones, the results numerically obtained. (a), (b) FD\_eq\_11p6o; (c), (d) FD\_s11p10o; (e), (f) FD\_s7p6o. Inside the source zone, enclosed by the red square, only the analytical solution is defined.</figcaption></figure>
Firstly, the standard *FDs7p6o* scheme previously used is here replaced with one of the two equivalent IOFD ones introduced earlier (*FDeq11p6o* and *FDeq15p8o*). In addition, the filtering stage is handled with a centered standard 21-point (rather than the previous 11-point) filter, whose very high accuracy guarantees that the possible bias to be introduced by the filtering operation remains negligible. Finally, the acoustic source frequency previously used ($k = 2\pi$) is increased significantly, so that it reaches either $k = 5\pi$ or $k = 6\pi$. Regarding the computational setup used, here,

<figure><img src="image_16.png" /><figcaption>Figure 14. Comparison between different finite-difference schemes for the sphere test case with $k = 6\pi$. In solid lines, the analytical result and, in dashed ones, the results numerically obtained. (a), (b) *FDeq11p6o*; (c), (d) *FDs11p10o*; (e), (f) *Fs7p6o*. Inside the source zone, enclosed by the red square, only the analytical solution is defined.</figcaption></figure>
both signals are thus spatially discretized with 4 PPW and 3.3 PPW only. Please, note that these two
values were chosen so as to discriminate between the various FD schemes assessed here, knowing
that an FD scheme which cannot accurately capture a signal to be CAA-forced is likely to generate
spurious modes [23, 24, 28].

As was established earlier, with respect to the present computational setup, the first signal ($k = 5\pi$) is spatially discretized with 4 PPW only. For such a low PPW value, Table IV indicates that the equivalent IOFD scheme *FDeq11p6o* should provide more accurate propagation results than its standard counterpart of same stencil width or accuracy order (*FDs7p6o* and *FDs11p10o*, respectively). This is indeed confirmed by the comparison of numerical results delivered by the Analytical-CAA coupled calculations associated with each one of all three previous schemes. Indeed, as can be seen in Figure 13, which depicts the instantaneous perturbed fields obtained within the planes $z = 0$ and $x = 0$, the equivalent IOFD scheme, *FDeq11p6o*, considerably improves the propagation results, compared with its two standard counterparts.

Regarding now the second source signal ($k = 6\pi$), and as can be seen in Figure 14, the previous equivalent IOFD scheme, *FDeq11p6o*, keeps its superiority over its two standard counterparts (*FDs7p6o* and *FDs11p10o*). However, it can be seen that here, such a *FDeq11p6o* scheme has reached its accuracy limit (with only 3.3 PPW, see Section 4.3). When now replaced with its 15-point stencil equivalent (*FDeq15p8o*), the expectations of Table IV are verified; in such an extreme configuration, the latter scheme still delivers very accurate results, as can be seen in Figure 15.

Based on all the previous results, it can be concluded that the equivalent IOFD scheme *FDeq11p6o* outperforms any of its standard counterparts, allowing a much coarser grid to be used,

<figure><img src="image_17.png" /><figcaption>Figure 15. Comparison between $k = 5\pi$ (top) and $k = 6\pi$ (bottom) for the scheme *FDeq15p8o*. In solid lines, the analytical result and, in dashed ones, the results numerically obtained. Inside the source zone, enclosed by the red square, only the analytical solution is defined.</figcaption></figure>
while inducing no or only a small amount of extra cost. In the same way, the intrinsically-optimized scheme FDeq15p8o appears to be the one to be used whenever very highly accurate simulations and/or highly coarser grids are to be handled.

## 7. CONCLUSION AND PERSPECTIVES

The present work highlighted the main limitation associated with the spectral-like optimization approach, that is, the order-of-accuracy loss inherited from the optimization procedure. In order to overcome such a limitation, an innovative optimization technique of FD derivative schemes was proposed. The optimization procedure was developed thanks to a rigorous analysis of Lagrange polynomial functions as well as to a judicious use of the particular properties that these functions offer.

In contrast to the spectral-like optimization techniques that can be found in the literature, the present optimization is carried out in the real space, rather than in the spectral one. Additionally, whereas spectral-like optimizations seek to balance several aspects (scheme accuracy order, fidelity with respect to harmonic signals, etc.), the present optimization solely consists of minimizing the scheme's leading-order truncation error. In other words, the main advantage of the optimization technique proposed here lies in the fact that no particular constraint is imposed on the objective function of the optimization process, which makes it potentially more generic/universal than others (such as the spectral-like one, which makes explicit assumptions about frequency ranges to be privileged, etc.). Indeed, the sole penalty that weighs on the present optimization technique lies in its interpolation stage. However, the latter could potentially be improved further at no cost, for example, by replacing the Lagrangian basis that was privileged here with another one.

As a result, when compared with standard schemes of same stencil width or formal order-of-accuracy, the IOFD schemes proposed here enable to solve more accurately harmonic signals of high-frequency content, with a minimal discretization criteria of as little as 4 PPW. On another hand, when compared with spectral-like optimized schemes of same stencil width, the present IOFD schemes ensure a fast error level reduction when the mesh grid resolution increases.

As an example, with no more than 4 PPW, the so-called equivalent IOFD scheme of 11 (resp. 15) points guarantees that the error made on the group velocity (on which depends the acoustic energy transport and, thus, the overall accuracy of the CAA stage) is less than 1% (resp. 0.14%). Compared with the centered standard 7-point 6-order FD scheme, which same error level corresponds to 8 PPW (resp. 12 PPW), such an IOFD scheme offers a gain of 2 (resp. 3), in terms of grid density per direction, that is, a factor of 8 (resp. 27) in 3D.

In addition to that, thanks to a higher formal order-of-accuracy, IOFD schemes are particularly well-suited for handling not only harmonic signals but also signals of very different nature and/or order of magnitude (e.g., hydrodynamic occurrences). This feature makes these schemes better adapted to the solving of complex aeroacoustics problems as well as offering interesting potentialities for an application within other fields of research (such as the dynamic LES subgrid modeling [11]).

# APPENDIX A: DISCRETE APPROXIMATION OF DIFFERENTIAL OPERATORS VIA EXPLICIT FINITE-DIFFERENCE SCHEMES

## A.1. Theoretical background

Let us consider the exact $m^{th}$-order partial differential operator $\partial_x^m[\cdot]$. For a given 1D stencil of uniform discretization step $\Delta x$, we can approximate this exact differentiation with the help of a discrete operator, $\partial_x^m[\cdot]$, to be based on explicit FD schemes:

$$ \tilde{\partial}_x^m[\cdot] = \frac{1}{(\Delta x)^m} \sum_{j=-M}^{N} b_{j,m} \tau_{(j\Delta x)}[\cdot]. \quad (A.1) $$
In the aforementioned equation, $b_{j,m}$ indicate the scheme's characteristic coefficients, whose values are still to be determined. Additionally, $\tau$ is the translation operator that, for any continuous quantity $f: \mathbb{R} \to \mathbb{R}$ and a given point $x \in \mathbb{R}$, delivers

$$ \tau_h[f](x) = f(x + h), h \in \mathbb{R}. \qquad (A.2) $$

At this stage, one can notice that centered (resp. noncentered) schemes are such that $N = M$ (resp. $N \neq M$). One can also remark that, for $m = 1$, (A.1) is associated with the well-known FD first derivative in space (with $b_{j,1}$ coefficients usually indicated as $a_j$). On another hand, for $m = 0$, the same equation is this time associated with the FD filter in space (with $b_{j,0}$ coefficients usually indicated as $f_j$).

## A.2. Error operator

The formal error operator resulting from the aforementioned approximation of the exact $m^{th}$-order differential operator is defined as

$$ \epsilon_m[\cdot] = \partial_x^m[\cdot] - \tilde{\partial}_x^m[\cdot]. \qquad (A.3) $$

By inserting (A.1) in (A.3), and by considering the Taylor expansion of the translation operator

$$ \tau_{(j\Delta x)}[\cdot] = \sum_{p=0}^{+\infty} \frac{(j\Delta x)^p}{p!} \partial_x^p[\cdot], \qquad (A.4) $$

one gets

$$ \epsilon_m[\cdot] = \frac{1}{(\Delta x)^m} \sum_{p=0}^{+\infty} (\delta_{m,p} - \chi_p) \Delta x^p \partial_x^p[\cdot], \quad \text{with } \chi_p = \frac{1}{p!} \sum_{j=-M}^{N} b_{j,m} j^p. \qquad (A.5) $$

In the aforementioned equation, $\delta_{m,p}$ denotes the Kronecker symbol, which is of nonzero and unitary value if and only if $m = p$.

When applied to a given function, $f: \mathbb{R} \to \mathbb{R}$, the operator $\tilde{\partial}_x^m[\cdot]$ induces a specific error, which results from both the formal error operator and the function itself:

$$ \epsilon_m[f](x) = \frac{1}{(\Delta x)^m} \sum_{p=0}^{+\infty} (\delta_{m,p} - \chi_p) \Delta x^p \partial_x^p[f](x). \qquad (A.6) $$

Following (A.6), the discrete differential operator $\tilde{\partial}_x^m[\cdot]$ is said to be $P^{th}$-order accurate if $\epsilon_m[f](x) = O(\Delta x^P)$ when $\Delta x \to 0$, which occurs when the $P + m$ following relations are verified simultaneously:

$$ \chi_p = \delta_{m,p}, \text{ for } p = 0, 1, \dots, P + (m - 1). \qquad (A.7) $$

From the previous relations, one can infer that the maximal order-of-accuracy to be possibly reached by a given explicit FD scheme is $P = N + M$. In particular, for $m = 1$ (resp. $m = 0$), such a scheme will meet its maximal order-of-accuracy only if it is of both a centered and antisymmetric (resp. symmetric) nature, that is, $N = M$ and $b_{j,1} = -b_{-j,1}$ (resp. $b_{j,0} = b_{-j,0}$). The set of coefficients $b_{j,m}$ ensuring that the corresponding scheme is effectively $2N^{th}$-order accurate can then be determined through the inversion of the linear system associated with (A.7).

## A.3. Spectral errors

Let us consider a given function $f: \mathbb{R} \to \mathbb{R}$, whose Fourier transform exists and is defined by $\mathcal{F}[f](\alpha) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} f(x)e^{-i\alpha x} dx$, with $\alpha$ a given wavenumber. When considering the Fourier transform of (A.6), the Derivative and the Shifting Theorems [29] give

$$ \mathcal{F}[\epsilon_m[f]](\alpha) = \frac{1}{(\Delta x)^m} \sum_{p=0}^{+\infty} (\delta_{m,p} - \chi_p)(i\alpha\Delta x)^p \mathcal{F}[f](\alpha), \forall \alpha \in [-\infty, +\infty]. \qquad (A.8) $$
This expression can be further simplified by considering, for a given $j$, the expansion of the exponential function:

$$e^{ij\alpha\Delta x} = \sum_{p=0}^{+\infty} \frac{(ij\alpha\Delta x)^p}{p!}, \qquad (A.9)$$

and then, by rewriting the expression in (A.8) as

$$\mathcal{F}[\epsilon_m[f]](\alpha) = \frac{i^m}{(\Delta x)^m} ((\alpha\Delta x)^m - (\tilde{\alpha}_m\Delta x)^m) \mathcal{F}[f](\alpha), \forall \alpha \in [-\infty, +\infty], \qquad (A.10)$$

with

$$\tilde{\alpha}_m(\alpha\Delta x)\Delta x = \left[ -i^m \sum_{j=-M}^{N} b_{j,m} e^{i(j\alpha\Delta x)} \right]^{1/m}. \qquad (A.11)$$

In all of the following, we restrict ourselves to the study of FD schemes associated with the first derivative, that is, those discrete operators such that $m = 1$. We also voluntarily focus solely on those FD schemes that are of a nondissipative nature. Being associated with centered antisymmetric schemes, these operators are characterized by a purely real (i.e., attenuation free [3]) $\tilde{\alpha}_1$:

$$\tilde{\alpha}_1(\alpha\Delta x)\Delta x = 2 \sum_{j=1}^{N} b_{j,1} \sin(j\alpha\Delta x), \forall \alpha \in [-\infty, +\infty]. \qquad (A.12)$$

The latter quantity, which directly drives the spectral error (see (A.10)), then becomes a measure of the so-called dispersion error $e_d$, defined as

$$e_d(\alpha\Delta x) = \alpha\Delta x - \tilde{\alpha}_1\Delta x. \qquad (A.13)$$

Such $e_d$ quantity is particularly well suited for assessing the efficiency of a given FD scheme regarding an application of the latter to wavelike phenomena. Indeed, any aeroacoustic signal can be decomposed into a sum of harmonic components that, once projected onto a Cartesian axis (the x-axis, for example), take the form $f(x) = e^{ikx}$, with $k \in \mathbb{R}$ the wave number. For each one of these harmonic components, (A.10) gives

$$\mathcal{F}[\epsilon_1[f]](\alpha) = \frac{i}{\Delta x} e_d(\alpha\Delta x)\delta(\alpha - k), \forall \alpha \in [-\infty, +\infty], \qquad (A.14)$$

showing that, in this case, the spectral error (A.14) exactly matches the dispersion error (A.13), to a factor of $i/\Delta x$ excepted.

From what precedes, one can easily understand how important it is for FD schemes that are used within a CAA context to exhibit a low spectral error. Indeed, this constitutes a necessary condition for these schemes to induce only a minimal dispersion of harmonic-like elementary components of the acoustic signal to be simulated.

However, although necessary, the previous condition is not sufficient. Indeed, when it comes to accurately simulating wave phenomena, besides the dispersion error, another important measure of a scheme's efficiency is its ability to accurately render the group velocity, which drives the acoustic energy transport. The group velocity error (which notion is legitimate only when the governing equations are linear [3]) can be written as [30]

$$e_g(\alpha) = 1 - \frac{\partial \tilde{\alpha}}{\partial \alpha}. \qquad (A.15)$$
APPENDIX B: INTRINSIC OPTIMIZATION PROCEDURE FOR CENTERED FD SCHEMES

Regarding centered FD schemes, the intrinsic optimization procedure can be sketched as a succession of four steps, which are detailed hereafter (only with respect to the direction $\tilde{x}$, for the sake of clarity).

1. **Evaluation of the finite-difference coefficients associated with the Chebyshev nodes.** Regarding the Chebyshev points $\{\tilde{x}_j\}_{j \in D_{2N+1}}$ included in $[-1, 1]$, one first needs to identify the particular set of coefficients $a_j^c, j \in D_{2N+1}$, corresponding to the discrete operator in the Chebyshev transformed space. This can be carried out by analytically differentiating, for each $j$, the Lagrange polynomial with respect to the central point $\tilde{x}_0 = 0$:

$$a_j^c = \partial_{\tilde{x}} [\tilde{P}_j](0) = \prod_{p \in D_{2N+1}, p \neq j, p \neq 0} \frac{-\tilde{x}_p}{\tilde{x}_j - \tilde{x}_p}. \quad (B.1)$$

Table B.1 delivers the particular sets of coefficients obtained for the Chebyshev nodes associated with various stencils of different sizes. In the present paper, the $(2N + 1)$-point $(2N)^{th}$-order accurate intrinsically optimized schemes are referred to as FDC $(2N + 1)p(2N)\rho$.

2. **Evaluation of the Jacobian transformation** that maps the discretization points, $(x, y, z)$, onto the Chebyshev nodes, $(\tilde{x}, \tilde{y}, \tilde{z})$. The Jacobian transformation can be more easily derived from its inverse, which is defined as

$$\mathbf{J}^{-1}(\tilde{x}) = \frac{\partial(x, y, z)}{\partial(\tilde{x}, \tilde{y}, \tilde{z})}. \quad (B.2)$$

As an illustration, let us consider a set of equally spaced points corresponding to a Cartesian grid in the physical space, with $\Delta x = \Delta y = \Delta z = 1$ (Figure B.1). In this case, the inverse Jacobian transformation reduces to a diagonal matrix $\mathbf{J}^{-1} = \mathbf{J}^{-1} \mathbf{Id}$, with $\mathbf{Id}$ being the identity matrix. By making use of the coefficients listed at Table B.1, along with the FD scheme formula associated with each coordinate, $\mathbf{J}^{-1}$ can be numerically evaluated at the point $\tilde{x}_0 = 0$:

$$\mathbf{J}^{-1} = \sum_{j=-N}^{N} a_j^c j \Delta x. \quad (B.3)$$

Because such evaluation is carried out on the Chebyshev nodes, its truncation error is minimal. Finally, $\mathbf{J}$ is obtained through the relation $\mathbf{J}\mathbf{J}^{-1} = \mathbf{Id}$.

3. **Evaluation of the function values associated with the Chebyshev nodes.** Because the approximate derivative is to be virtually solved on the Chebyshev nodes, it is necessary to identify which values the function $f(\cdot)$ may take on those points. A straightforward way of doing that would consist of mapping the transformed Chebyshev nodes onto the discrete Cartesian

<table><caption>Table B.1. Stencil coefficients associated with centered finite-difference derivative schemes based on Chebyshev nodes.</caption><thead><tr><th></th><th>FDc7p6o</th><th>FDc9p8o</th><th>FDc11p10o</th><th>FDc13p12o</th></tr></thead><tbody><tr><td>a<sub>0</sub><sup>c</sup></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>a<sub>1</sub><sup>c</sup> = -a<sub>-1</sub><sup>c</sup></td><td>2.07652139657234</td><td>2.74747741945462</td><td>3.40568723888926</td><td>4.05715948563811</td></tr><tr><td>a<sub>2</sub><sup>c</sup> = -a<sub>-2</sub><sup>c</sup></td><td>-0.79747338882404</td><td>-1.19175359259421</td><td>-1.55603037297544</td><td>-1.90534081130307</td></tr><tr><td>a<sub>3</sub><sup>c</sup> = -a<sub>-3</sub><sup>c</sup></td><td>0.228243474390150</td><td>0.577350269189626</td><td>0.866504932543064</td><td>1.12876666008794</td></tr><tr><td>a<sub>4</sub><sup>c</sup> = -a<sub>-4</sub><sup>c</sup></td><td></td><td>-0.176326980708465</td><td>-0.456684697903556</td><td>-0.690250161969227</td></tr><tr><td>a<sub>5</sub><sup>c</sup> = -a<sub>-5</sub><sup>c</sup></td><td></td><td></td><td>0.143778293994974</td><td>0.379249975405874</td></tr><tr><td>a<sub>6</sub><sup>c</sup> = -a<sub>-6</sub><sup>c</sup></td><td></td><td></td><td></td><td>-0.121421983222261</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Figure B.1. Illustration of a 2D Cartesian grid defined as the physical domain (left side) and its corresponding numerical domain composed of Chebyshev nodes (right side) to be used with the centered five-point finite-difference scheme. The corresponding Jacobian mapping, J, is also represented. The filled circles stand for the Cartesian grid points that are mapped onto the Chebyshev nodes.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure B.2. 2D Cartesian grid of Figure B.1 and its stencil-corresponding Chebyshev nodes after the latter have been submitted to the inverse Jacobian transformation, J⁻¹.</figcaption></figure>

points, this being carried out by submitting the former to the inverse Jacobian transformation. However, such an operation is likely to be jeopardized by the fact that the Jacobian transformation may not be bijective for certain configurations. In that case, as can be seen in Figure B.2, the resulting mapping is inaccurate, which does not allow for properly allotting the function values f(\tilde{x}_j) with the f(x_j) ones, ∀j ≠ 0.

In order to be able to virtually differentiate the function f(.) on the Chebyshev points, it is thus necessary to calculate its values f(\tilde{x}_j), j ∈ D₂N+1 and j ≠ 0, in a different manner such as the one proposed here, which consists of interpolating the function's unknown values f(\tilde{x}_j) directly from the known ones, f(x_j). In the present work, such an interpolation step is achieved thanks to centered Lagrange interpolations, which offer the advantage of being amplification-free [21]. Following that, for any given j ∈ D₂N+1, the missing function value f(\tilde{x}_j) is approximated by a centered nᵗʰ-order Lagrange polynomial (with n not necessarily equal to 2N + 1):

$$f(\tilde{x}_j) \cong \mathcal{I}[f](\mathbf{J}^{-1}\tilde{x}_j) = \sum_{k \in I_n} P_k(\mathbf{J}^{-1}\tilde{x}_j) f(x_k), \quad (B.4)$$

where the suite {$x_k$}$_{k \in I_n}$ corresponds to the discrete equally spaced points associated with the Cartesian grid. Figure B.3 depicts the typical donor stencil that would be involved if the interpolation step is achieved with the help of a centered fourth-order Lagrange interpolator. Obviously, depending on the interpolation basis, additional donor grid points may be needed.
<figure><img src="image_20.png" /><figcaption>Figure B.3. Donor (Cartesian) and receiver (J⁻¹-transformed Chebyshev) grid points depicted in Figure B.2. In unfilled circles, the additional donor points needed for a centered fourth-order Lagrange interpolation.</figcaption></figure>

4. Evaluation of the approximate derivative on the Chebyshev nodes. Once all the previous steps are achieved, one can evaluate the approximate value of the spatial differential operator in a straightforward manner following (17).

$$ \partial_x[f](x_0) \approx J \sum_{j \in D_{2N+1}} a_j^c \mathcal{I}[f](\mathbf{J}^{-1}\tilde{x}_j). \quad (B.5) $$

## APPENDIX C: CALCULATION OF THE EQUIVALENT INTRINSICALLY-OPTIMIZED FINITE-DIFFERENCE SCHEMES

Rewriting the right-hand side of (B.5) in the equivalent form (19) implies the identification of a set of equivalent coefficients, which can be derived by making use of the following algorithm.

**input**: optimization parameters N, n and N<sup>eq</sup>; Chebyshev coefficients a<sub>j</sub><sup>c</sup>; Cartesian nodes x<sub>j</sub>; Chebyshev nodes x̃<sub>j</sub>

**output**: equivalent coefficients a<sub>j</sub><sup>eq</sup>

1. $a^{eq} := size(2N^{eq} + 1, 1)$

2. $\xi := size(2N + 1, 1)$

3. **for** j ∈ D<sub>2N+1</sub> **do**

4. **for** k ∈ I<sub>n</sub> **do**

5. $a_{\text{FLOOR}(\tilde{x}_j)-\text{FLOOR}(\tilde{x}_1)+1+k-\min(I_n)}^{\text{eq}} := J * a_j^c * P_k(\tilde{x}_j - \text{FLOOR}(\tilde{x}_j)) + a_{\text{FLOOR}(\tilde{x}_j)-\text{FLOOR}(\tilde{x}_1)+1+k-\min(I_n)}^{\text{eq}}$

6. **end for**

7. **end for**

In the aforementioned algorithm, the floor function FLOOR(z), z ∈ Z, rounds z to the nearest integer less than, or equal to, z. Additionally, min(I<sub>n</sub>) represents the smallest integer included in I<sub>n</sub>.

## APPENDIX D: EXTENSION TO THE SOLVING ON CURVILINEAR GRIDS

Sensu stricto, deriving finite-difference scheme coefficients a<sub>j</sub> through the inversion of the linear system associated with Equation (A.7) is valid only when discrete points are of a 1D equally-spaced nature. However, in real-life, it is generally needed to discretize and numerically solve problems with the help of multidimensional curvilinear grids (Figure D.1, left side, for a 2D example). Therefore, in order to avoid the issue of deriving specific sets of coefficients a<sub>j</sub> for each grid point and
<figure><img src="image_21.png" /><figcaption>Figure D.1. The 2D curvilinear discretization grid (left side) and its corresponding derivation mesh composed of Chebyshev nodes (right side), in the particular case of a five-point centered finite-difference scheme. The filled circles stand for the curvilinear grid points to be mapped onto the Chebyshev nodes via the Jacobian transformation, **J**.</figcaption></figure>

direction, it is common to map the physical curvilinear mesh onto a numerical grid, over which the regular finite-difference schemes can then be applied. In the present case, where the numerical grid is locally defined with the Chebyshev points, this technique still holds. One can then submit each space coordinate to a time-invariant curvilinear coordinate transformation **J**: $x, y, z \mapsto \tilde{x}, \tilde{y}, \tilde{z}$ (see (16)) that maps the discretized (i.e. physical) space, $(x, y, z)$ onto the derivation (i.e. numerical) one, $(\tilde{x}, \tilde{y}, \tilde{z})$.

## REFERENCES

1. Tam CKW. Computational aeroacoustics: An overview of computational challenges and applications. *Journal of Computational Fluid Dynamics* 2004; **18**(6):247-567.

2. Sherer SE, Scott JN. High-order compact finite-difference methods on general overset grids. *Journal of Computational Physics* 2005; **210**(2):459-496.

3. Colonius T, Lele SK. Computational aeroacoustics: Progress on nonlinear problems of sound generation. *Progress in Aerospace Sciences* 2004; **40**(6):345-416.

4. Lele SK. Compact finite difference schemes with spectral-like resolution. *Journal of Computational Physics* 1992; **103**(1):16-42.

5. Haras Z, Ta'asan S. Finite difference schemes for long-time integration. *Journal of Computational Physics* 1994; **114**(2):265-279.

6. Kim JW, Lee DJ. Optimized compact finite difference schemes with maximum resolution. *AIAA Journal* 1996; **34**(5):887-893.

7. Mahesh K. A family of high-order finite difference schemes with good spectral resolution. *Journal of Computational Physics* 1998; **145**(1):332-358.

8. Hixon R. Prefactored small-stencil compact schemes. *Journal of Computational Physics* 2000; **165**(2):552-541.

9. Lockard DP, Brentner KS, Atkins HL. High-accuracy algorithms for computational aeroacoustics. *AIAA Journal* 1995; **33**(2):246-251.

10. Tam C KW, Webb JC. Dispersion-relation-preserving finite difference schemes for computational acoustics. *Journal of Computational Physics* 1993; **107**(8):262-281.

11. Bogey C, Bailly C. A family of low dispersive and low dissipative explicit schemes for flow and noise computations. *Journal of Computational Physics* 2004; **194**(1):194-214.

12. Cunha G, Redonnet S. On the effective accuracy of spectral-like optimized finite-difference schemes for computational aeroacoustics. *Journal of Computational Physics* 2014; **263**(1):222-232.

13. Baeza A, Mulet P. Adaptive mesh refinement techniques for high-order shock capturing schemes for multi-dimensional hydrodynamic simulations. *International Journal for Numerical Methods in Fluids* 2006; **52**(4):455-471.

14. Mulet P, Baeza A. Highly accurate conservative finite difference schemes and adaptive mesh refinement techniques for hyperbolic systems of conservation laws. *Numerical Mathematics and Advanced Applications - Proceedings of ENUMATH 2005, the 6th European Conference on Numerical Mathematics and Advanced Applications Santiago de Compostela*, Spain, July 2005; 198-206.

15. Shen C, Qiu JM, Christlieb A. Adaptive mesh refinement based on high order finite difference weno scheme for multi-scale simulations original research article. *Journal of Computational Physics* 2011; **230**(10):3780-3802.
16. Yoon D, Kim H, Hwang W. Adaptative mesh refinement for weighted essentially non-oscillatory schemes. *Bulletin of the Korean Mathematical Society* 2008; **45**(4):781–795.

17. Tam C KW, Kurbatskii KA. Multi-size-mesh multi-time-step dispersion-relation-preserving scheme for multiple-scales aeroacoustics problems. *International Journal of Computational Fluid Dynamics* 2003; **17**(2):119–132.

18. Bogey C, Bailly C. High-order, low dispersive and low dissipative explicit schemes for multiple-scale and boundary problems. *Journal of Computational Physics* 2007; **24**(2):637–662.

19. Davis PJ. *Interpolation and Approximation*. Dover Publications, 1975.

20. Schmid P. *An Introduction to Computational Fluid Dynamics*. John Wiley and Sons Inc.: Mineola, New York, 2007.

21. Cunha G, Redonnet S. An innovative interpolation technique for aeroacoustic hybrid methods. *AIAA Paper 2011-2754*: Hoboken, New Jersey.

22. Tam CKW, Kurbatskii KA. An optimized extrapolation and interpolation method for computational aeroacoustics. *AIAA Paper 2001-0282*.

23. Cunha G, Redonnet S. On the signal degradation induced by the interpolation and the sampling rate reduction in aeroacoustics hybrid methods. *International Journal for Numerical Methods in Fluids* 2013; **71**(7):910–929. DOI: 10.1002/fld.3693.

24. Redonnet S, Cunha G. Towards a robust and accurate cfd-caa coupling procedure for hybrid methods in aeroacoustics. part 2 - on the application of the cfd-caa surface weak coupling methodology to realistic aircraft noise problems. *AIAA Paper 2012-2191*.

25. Redonnet S. On the numerical prediction of aerodynamic noise via a hybrid approach - part 1. *AIAA Paper 2010-3709*.

26. Redonnet S, Lockard D, Khorrami DP, Choudhari M. Cfd-caa coupled calculations of a tandem cylinder configuration to assess facility installation effects. *AIAA Paper 2011-2841*.

27. Redonnet S, Mincu D, Delattre G. Computational aeroacoustics of a realistic co-axial engine possibly equipped with acoustic liners. *AIAA Paper 2010-3717*.

28. Cunha G, Redonnet S. Towards a robust and accurate cfd/caa coupling procedure for hybrid methods in aeroacoustics. part 1: On the optimization of cfd/caa coupled calculations. *AIAA Paper 2012-2063*.

29. Boyd JP. *Chebyshev and Fourier Spectral Methods*. Springer-Verlag: Berlin Heidelberg, 1989.

30. Tam C KW, Webb JC, Dong Z. A study of the short wave components in computational acoustics. *Journal of Computational Acoustics* 1993; **1**(1):262–281.