One- and multi-dimensional CWENOZ reconstructions for implementing boundary conditions without ghost cells

We address the issue of point value reconstructions from cell averages in the context of third order finite volume schemes, focusing in particular on the cells close to the boundaries of the domain. In fact, most techniques known in the literature rely on the creation of ghost cells outside the boundary and on some form of extrapolation from the inside that, taking into account the boundary conditions, fills the ghost cells with appropriate values, so that a standard reconstruction can be applied also in boundary cells. In (Naumann, Kolb, Semplice, 2018), motivated by the difficulty of choosing appropriate boundary conditions at the internal nodes of a network, a different technique was explored that avoids the use of ghost cells, but instead employs for the boundary cells a different stencil, biased towards the interior of the domain. In this paper, extending that approach, which does not make use of ghost cells, we propose a more accurate reconstruction for the one-dimensional case and a two-dimensional one for Cartesian grids. In several numerical tests we compare the novel reconstruction with the standard approach using ghost cells.


Introduction
Computing in an efficient way accurate albeit non-oscillatory solutions of conservation laws requires the employment of high-order accurate numerical schemes.
The issue of boundary treatment for hyperbolic conservation laws is usually tackled by constructing ghost points or ghost cells outside the computational domain and by setting their values with appropriate techniques; after this, a high-order non-oscillatory reconstruction procedure can be applied also close to the boundary, despite its large stencil, thanks to the ghost values. This approach is of course delicate, especially with finite-difference discretizations on non-conforming meshes. In this context a very successful technique is the Inverse Lax-Wendroff approach introduced in [36], rendered more computationally efficient in [37], and further studied and extended for example in [22,24,25]; a quite up-to-date review may be found in [35]. A modified procedure enhancing its accuracy and stability has been proposed in [39]. Other approaches, still based on an inverse Lax-Wendroff procedure but more tailored to coupling conditions on networks can be found in [7,11]. A different approach, entirely based on WENO extrapolation is studied in [2,3].
In [26] a different approach was considered. There, in a one-dimensional finite volume context, ghost values are entirely avoided and the point value reconstruction at the boundary is performed with a CWENO type non-oscillatory reconstruction that makes use of only interior cell averages. The reconstruction stencil for the last cell at the boundary is not symmetric, but extends only towards the interior of the computational domain. Then, the boundary flux is determined from the reconstructed value and the boundary conditions. Achieving non-oscillatory properties when a discontinuity is close to the boundary requires the inclusion of very low degree (down to a constant one, in fact) polynomials in the CWENO procedure; this, in turn, calls for infinitesimal linear weights in order not to degrade the accuracy on smooth solutions. This type of CWENO reconstructions have been studied in general in [30].
In this paper, we first enhance the accuracy of the boundary treatment of [26] by employing an Adaptive Order CWENO-Z reconstruction from [30] and furthermore extend it to two space dimensions. In particular, in §2 we describe the new one-dimensional reconstruction that avoids ghost cells and, in §3, compare it with the one of [26] on numerical tests. The novel two-dimensional no-ghost reconstruction is then presented in §4 and the corresponding numerical results that compares it with the ghosted approach of [14] are presented in §5. Some final remarks and conclusions are drawn in §6.

The novel CWENOZb reconstruction in one space dimension
We start recalling here the operators that define a generic Central WENO reconstruction, which will be useful later. Central WENO is a procedure to reconstruct point values of a function from its cell averages; it is different from the classical WENO by the fact that it performs a single non-linear weight computation per cell and outputs a polynomial globally defined in the cell, which is later evaluated at reconstruction points.
In defining a Central WENO, one starts selecting an optimal polynomial, denoted here by P opt , which should be chosen to have the maximal desired accuracy; the CWENO reconstruction polynomial, in fact, will be very close to this one when the cell averages in the stencil are a sampling of a smooth enough function. For the cases when a discontinuity is present in the stencil of P opt , a sufficient number of alternative polynomials, P 1 , . . . , P m , typically with lower degree and with a smaller stencil, should be made available to the blending procedure.
The CWENO operator than computes a nonlinear blending of all polynomials as follows. First a set of positive linear or optimal coefficients is chosen, with the only requirement that d 0 + d 1 + . . . + d m = 1. Then, the reconstruction polynomial is defined by The quantities ω i appearing above are called nonlinear weights; when ω i = d i for i = 0, . . . , n, then P rec = P opt and the reconstruction will have the maximal accuracy. When a discontinuity is present in the stencil, the nonlinear weight should deviate from their optimal values in order to avoid the occurrence of spurius oscillations in the numerical scheme. In practice, the nonlinear weights are computed with the help of oscillation indicators associated to each polynomial, that should be o(1) when the polynomial interpolates smooth data and O(1) when the polynomial interpolates discontinuous data. The construction is independent from the specific form of these indicators, which here we denote generically as OSC[P ]; typically, the Jiang-Shu indicators from [19] are employed. When the reconstruction is denoted by CWENO(P opt ; P 1 , . . . , P m ), the nonlinear coefficients are computed as in the original WENO construction, namely as where is a small parameter and p ≥ 1. For detailed results on the accuracy of such a reconstruction, see [12] and the references therein. Better accuracy on smooth data, especially on coarse grids, without sacrificing the non-oscillatory properties, can be obtained by computing the nonlinear weights as in the WENOZ construction of [15], namely as In this case, we denote the reconstruction as CWENOZ(P opt ; P 1 , . . . , P m ). Here above, τ is quantity that is supposed to be much smaller than the individual indicators when the data in the entire reconstruction stencil are smooth enough. For efficiency, this global smoothness indicator should be computed as a linear combination of the other oscillators. For results on the optimal choices for τ in a CWENO setting and the accuracy of the resulting reconstructions, see [14] and the references therein. The accuracy results of both CWENO and CWENOZ require that certain relations among the accuracy of all polynomials involved are satisfied; the precise conditions for optimal accuracy involve also the parameters and p [12,14], but the as a rule of thumb one should always have deg(P opt ) ≤ 2deg(P k ) for k = 1 . . . , m. If controlling spurious oscillations require the inclusion in the nonlinear combination of polynomials with degree smaller than 1 2 deg(P opt ), optimal accuracy can still be achieved if the linear weights of these additional polynomials of very low degree are infinitesimal, i.e. chosen as O(∆x r ) for some r > 0. In order to distinguish this situation and to easily spot the polynomials with infinitesimal linear weights, we adopt for this case the notations CWENO-AO(P opt ; P 1 , . . . , P m ; Q 1 , . . . , Q m ), when (2) is used for the nonlinear weights, and CWENOZ-AO(P opt ; P 1 , . . . , P m ; Q 1 , . . . , Q m ), when (3) is used instead. This was studied on a specific example in [26] for the CWENO case and in general for CWENOZ-AO in [30]. This latter contains a thorough study of sufficient conditions on r and on the other parameters that guarantee optimal convergence rates for a generic CWENOZ-AO reconstruction.

CWENO-boundary reconstruction of [26]
A third-order accurate reconstruction that does not make use of ghost cells has been introduced in [26]. The reconstruction coincides with the CWENO3 reconstruction of [21], with variable parameter as in [20,13,12]. In particular, for the j-th cell one considers the following polynomials: P (2) j , which is the optimal second degree polynomial interpolating u j−1 , u j , u j+1 , P (1) j,L and P (1) j,R , which are the linear polynomials interpolating u j−1 , u j and u j , u j+1 respectively. CWENO3 is a shorthand for CWENO(P j,R ). This reconstruction produces a second degree, uniformly third order accurate, polynomial defined in each cell, using the cell averages in a stencil of three cells; it can thus be computed on every cell in the domain except for the last one close to each boundary.
In the first cell of the domain, the reconstruction is replaced with an adaptiveorder reconstruction CWENO-AO(P   is the parabola interpolating the cell averages u 1 , u 2 , u 3 . The last cell is treated symmetrically. The inclusion of the constant polynomial P (0) is necessary to prevent oscillations when a discontinuity is present one cell away from the boundary and giving it an infinitesimal linear weight allows to guarantee the optimal order of convergence for the reconstruction procedure on smooth data. More precisely, in [26] it is shown that choosing the linear weights as d (0) = min(∆xm, 0.01) for the constant polynomial, d (1) = 0.25 for the linear one and consequently setting d 0 = 1 − d (1) − d (0) , guarantees the optimal accuracy on smooth data whenm ∈ [1, 2], provided = ∆x q with q ≥m.
In general a small yields good results on discontinuities, but keeping q = 1 seems desirable to avoid rounding problems in the computation of the nonlinear weights. The combinationm = 2 and q = 1, however does not fulfill the hypotheses of the convergence result of [26]; in practice, however, the reconstruction appears to give rise nevertheless to a third order accurate scheme but degraded accuracy can be observed at low grid resolutions. As an extreme example in this sense, let us consider the linear transport of a periodic initial datum in a periodic domain. Of course there would be no need to employ the no-ghost reconstruction in this case, since it would be trivial to fill in the ghost values (except maybe for considerations on parallel communication), but this example serves quite well to illustrate the situation on smooth data.
In Table 1 we report the 1-norm errors observed for the transport of u(x, 0) = sin(πx−sin(πx)/π) after one period (for full details on the numerical scheme, the reader is referred to the beginning of §3). It is evident that for the reconstruction of [26], third order error rates are observed only on very fine grids when d (0) ∼ ∆x; for d (0) ∼ ∆x 2 , the optimal rate predicted by the theory is observed in practice, but the errors are still larger than its ghosted CWENO3 counterpart. Figure 1: Illustration of the stencil for the 3-rd order reconstruction in the last cell. Blue: cell where the reconstruction is computed. Black: the polynomials involved in P rec . Gray: additional polynomial for τ (b3) .

The novel CWZb3 reconstruction
The loss of accuracy at low resolution can be traced back to the relative inability of the smoothness indicators alone to detect a smooth flow on coarse grids. The net effect is that, when the grid is coarse, the nonlinear weight of the constant polynomial in the first and last cells is larger than it would be strictly needed, degrading the accuracy of the reconstruction there; the errors are then transported into the domain by the flow.
This issue can be successfully counteracted, even on coarse grids, by the employment of Z-weights in the construction. In fact, we recall that the idea behind WENO-Z, see [15], is to replace the standard WENO nonlinear weight computation (2) with (3) where the global smoothness indicator τ is supposed to be τ = o(I k ) if the cell averages represent a locally smooth data in the stencil. The improved performances of WENO-Z over WENO, and of CWENOZ over CWENO reconstructions are in fact linked to the superior ability of detecting smooth transitions, already at low grid resolution, granted by the global smoothness indicator τ . Moreover, the ability of detecting a smooth flow even at low grid resolution depends on how small is τ on smooth data; thus the goal in the optimal design of τ is to choose the coefficients of the linear combination λ 0 OSC[P opt ] + n i=1 λ k OSC[P k ] = O(∆x s ) that maximize s when the data in the stencil of the reconstruction is a sampling of a smooth function [14].
Our proposal thus consists in defining the new CWZb3 reconstruction to coincide with CWENOZ3 = CWENOZ(P  [19] is where Ω j is the cell where the reconstruction is applied.
On smooth data, in the domain interior, we have that so that the combination is O(∆x 4 ); this very low τ biases very strongly the non-linear weights (3) towards the optimal ones whenever the flow is smooth. In [14] it is shown that this is the optimal choice and that it is not possible to obtain a combination of the indicators that is o(∆x 4 ) in the third order setup. We now need to specify a suitable τ 1 for the uncentered stencil of the first cell and τ N for the last one. Recall that the role of τ is to indicate whether the data are smooth in the stencil, which is composed by the first three cells adjacent to the boundary. As argued in [30], only the polynomials with degree at least one are useful in the construction of τ . Here we could use the oscillators of the parabola P (2) fitting the three cell averages u 1 , u 2 , u 3 and the linear polynomial P (1) 1 interpolating the first two, and thus we cannot exploit the symmetry to obtain a global smoothness indicator of size O(∆x 4 ).
Using τ of O(∆x 3 ) however could make the reconstruction in the boundary cell less performing than the one in the domain interior. In order to overcome this difficulty, one could employ, in the construction of τ , also the indicator of the linear polynomial P (1) interpolating the averages u 2 , u 3 . However, since the role of τ is to detect smooth flows in the global stencil, which is composed by the first three cells, that is the same reconstruction stencil employed by the second cell, a simpler solution (which also allows to save some computations) is to take instead for the first cell the same value of τ that was computed in the second cell; this is O(∆x 4 ) on smooth flows and yields a better reconstruction.
The novel reconstruction procedure that we propose is thus: • in all cells except the first and last one, compute the CWENOZ3 reconstruction polynomial with the optimal definition (4) of τ j , as in [14]; • in the first cell, apply CWENOZ-AO(P • in the last cell, apply CWENOZ-AO(P After the analysis of §3.1.1 of [30], it is expected that this reconstruction has the third order of accuracy for As discussed in [14], the choice of parameters within the allowed ranges can trade better accuracy on smooth flows (largerm or smaller q) with a smaller spurious oscillations on discontinuities (smallerm or larger q). In [14] it was found that a good overall choice for CWENOZ3 was p = 1 andm = 2 and we will adopt these values in all our numerical tests. Regarding the infinitesimal linear weight, the choices d (0) = ∆x and d (0) = ∆x 2 will be compared.

One-dimensional numerical tests
All tests in this section are conducted with a finite volume scheme constructed with the method of lines, the Local Lax-Friedrichs numerical flux, and the third order TVD-RK3. The CWENO3 and the CWENOZ3 reconstruction in the first/last cell make use of one ghost cell outside each boundary, which is filled according to the boundary conditions before computing the reconstruction. In the same cells, the CWb3 and the CWZb3 reconstructions, instead, do not make use of ghost cells but extend their stencil for one extra cells inwards with respect to their ghosted counterparts.
In both cases, the flux on the boundary face is computed by applying a consistent numerical flux (here the Local Lax Friedrichs) to an inner value determined by the reconstruction and an outer value determined by the boundary conditions. More precisely, for periodic boundary conditions, the outer value on the left is copied from the inner value at the right boundary and viceversa; for reflecting boundary conditions in gasdynamics, the outer value is the same as the inner one but has the opposite sign for the velocity; for Dirichlet boundary, the outer value is set to the exact value of the boundary function at time t n + c i ∆t for the i-th stage of the Runge-Kutta scheme.
The CFL number is set to 0.45 in all tests. The numerical tests have been performed with the open-source code claw1dArena, see [29].

Linear transport
Periodic solution We consider again the linear transport equation u t + u x = 0 in the domain [−1, 1] with periodic boundary conditions. We evolve for one period the initial data u 0 (x) = sin(πx − sin(πx)/π), which has a critical point of order 1 (see [16]), with the CWENOZ3 and CWZb3 reconstructions. Table 2 shows that, according to the results of [30], CWZb3 can reach the optimal convergence rate already with d (0) ∼ ∆x and that the errors obtained without using ghosts are very close to those of the ghosted reconstruction CWENOZ3. As already pointed out in [14], also here we observe that using Z-weights in CWENO yields lower errors compared to the companion reconstructions with Jiang-Shu weights (compare Tab. 1).  imposing u(−1, t) = 0.25 − 0.5 sin(π(1.0 + t) and free-flow conditions on the (outflow) boundary at x = 1. We start with u 0 (x) = 0.25 + 0.5 sin(πx) and compare the computed cell averages with the exact solution u(t, x) = u 0 (x − t).
The final time is set to 1. This test was proposed in [36]. The results reported in Tab. 3 show that the CWZb3 reconstruction yields third order error rates already on coarse grids and with d (0) ∼ ∆x. No advantage is seen for the choice d (0) ∼ ∆x 2 .
We point out that applying a reconstruction that makes use of ghost cells, like CWENO3 or CWENOZ3, would not be straightforward in this case. In [9] was observed that accuracy would be capped at second order if the ghost cell values for the i-th stage were to be set by reflecting the inner ones in the exact boundary data at time t n + c i ∆t, where c i denotes the abscissa of the i-th stage of the Runge-Kutta scheme. In the same paper, also a suitable modification of the boundary data that preserve the accuracy of the Runge-Kutta scheme is also proposed. On the other hand, we point out that with the CWb3 and CWZb3 reconstructions this issue of filling the ghost cells is not present and that the exact boundary data can be employed in the numerical flux computation without observing losses of accuracy.
Discontinous solution Next we consider the same setup of the previous test, but impose the boundary value thus introducing a jump in the exact solution at t = 1. This test was proposed in [36] and, as there, we compute the flow until t = 1.5.
The computed solutions are shown in Fig. 2, where we compare the solution computed with CWENOZ3 using ghosts and the no-ghost CWb3 and CWZb3. No difference can be seen in the corner point at x = 0.5, which is originated by a continous but not differentiable boundary data. On the other hand, the jump at x = −0.5 in the final solution is generated by the discontinuity in the boundary data. The numerical solution around this jump has slightly more pronounced oscillations when using CWZb3 and d (0) = ∆x 2 and a more smoothed profile when using CWb3; CWZb3 with d (0) = ∆x produces an almost idential solution to the one computed by the ghosted CWENOZ3 reconstruction.

Burgers' equation
For a nonlinear scalar test, we consider the Burgers' equation u t + (u 2 ) x = 0 with initial data u 0 (x) = 1 − sin(πx) with periodic boundary conditions, so that a shock forms, travels to the right and is located exactly on the boundary at t = 1.
In Fig. 3 we compare the solutions computed with 25 cells. One can see that CWZb3 computes a solution which is almost exactly superimposed on the CWENOZ3, despite the fact that using the correct periodic ghost values should

Euler gas dynamics
Incoming wave from the left In this test we consider a gas initially at rest, with ρ = 1, p = 1, u = 0 everywhere. Through a time-dependent Dirichlet boundary condition on the left, we introduce the following disturbance The boundary introduces a smooth wave travelling right. Wall boundary conditions are imposed on the right and the final time is set at t = 1.25, when the wave is being reflected back from the wall. In Fig. 4 we report the solutions at time t = 1.25 computed on 50 cells with the third order ghosted and ghost-free reconstructions, together with a reference solution computed on 10000 cells with a second order TVD scheme. Spurious oscillations coming from the Dirichlet boundary conditions on the left side are completely absent when using CWb3 or CWZb3 instead of CWENOZ3. Also, a slightly better resolution is observed near the top of the wave. Here again, we stress that the CWb3 and the CWZb3 solutions have been computed by entirely neglecting the boundary conditions in the reconstruction phase and passing the exact Dirichlet value at t n + c i ∆t to the numerical flux as outer data on the left boundary.
Sod's Riemann problem In this test we use the initial data of the Sod problem, but we impose wall boundary conditions on both sides.
In Fig. 5 we show, in the left panel, the solution at time t = 0.2, which is before the waves reach the wall; the expected solution is thus the usual one. All three solutions are very close to each other and only a slight extra diffusion can be noticed for the reconstruction that is using the Jiang-Shu nonlinear weights instead of the Z-weights.
Letting the flow evolve past t = 0.2, the shock impinges on the wall and bounces back, interacting with the right-moving contact around t = 0.41; this in turn generates a left-moving shock, a very slow contact and a quite weak right-moving shock; the right-moving shock then bounces back from the wall and interacts with the contact at around t = 0.56, giving rise again to another shock-contact-shock interaction pattern.
In the right panel of Fig. 5 we show the solution at t = 0.6, and, counting from the left, we see a rarefaction which is reflecting in the left wall, two leftmoving shocks, a very slow contact (with speed 0.001) and a very weak rightmoving shock (density jump below 0.01). It can be appreciated that CWZb3, without using ghost cells, computes almost the same solution as the ghosted CWENOZ3. As in other tests, CWb3 is more diffusive. The very weak shock is barely captured at this resolution and even the reference solution almost misses it. The oscillations in the plateaux between the two left-moving shocks and the  hump left of the contact could be controlled with local characteristic projections, which was not employed in this computation.
Finally, we consider the d-dimensional version of the same problem. Following [38], in spherical symmetry this amounts to adding to the Euler equations the source term S(ρ, u, p) = − d−1 x [ρu, ρu 2 , up] T . In particular we show in Fig. 6 the solution for d = 3 at t = 0.5 and at t = 0.65. In this test, the source term contribution is computed in each cell with a two-point gaussian quadrature, which is fed by the reconstructed values. We thus test the CWENO-based reconstructions' capability of easily computing reconstructed values inside the cells. For all waves, we observe again that CWZb3 and CWENOZ3 produce very similar solutions, with CWb3 being slightly more diffusive.

Two-dimensional scheme
In this section we consider a two-dimensional Cartesian grid, with cells of size ∆x. We denote the cells as Ω i,j , with the pair of integers (i, j) referring to their position in the grid. As in the one dimensional case, the solution is advanced in time with the third order TVD-SSP Runge-Kutta scheme; the numerical fluxes on each edge of a cell are obtained with the two-point gaussian quadrature, with point values computed with a two-point numerical flux fed with the reconstructed values on each side of the edge. On boundaries, only the inner point value is computed from the reconstruction, while the outer one is computed according to the boundary conditions. For example, on a solid wall boundary, the outer value is equal to the inner one, except for the normal velocity, which is given the opposite sign.
The reconstruction from cell averages to point values in two space dimensions is not obtained by dimensional splitting, but is computed by blending polynomials in two spatial variables with a CWENOZ or a CWENOZ-AO construction. The reconstruction operator is called only once per cell and the polynomial returned is then evaluated at the eight reconstruction points where the numerical fluxes have to be computed.
Let Ω i,j be the cell in which the reconstruction is being computed. In every cell, the reconstruction is computed by a CWENOZ operator with optimal polynomial P opt of degree 2 in two spatial variables (6 degrees of freedom) associated with a 3 × 3 stencil containing Ω i,j (see later for the definition of the polynomial associated to a stencil). The reconstruction stencils are depicted in Fig. 7. In all panels, the cell in which the reconstruction is being computed is hatched, while the stencil of the optimal polynomial of degree 2 is shaded.
The CWENOZ operator is fully specified after the low degree polynomials and the global smoothness indicator τ are also chosen. The stencils of the low degree polynomials are indicated by circles joined by solid or dashed lines in Fig. 7.
In the bulk of the computational domain, the reconstruction coincides with the two-dimensional CWENOZ3 described in [14]; it is defined as a nonlinear combination of second and first degree polynomials: CWENOZ(P opt ; P N E , P SE , P SW , P N W ).
The optimal polynomial is associated with the 3 × 3 stencil of cells centered at Ω i,j (left panel in Fig. 7). The four polynomials P N E , P SE , P SW and P N W are linear polynomials in two variables associated to the four stencils depicted with solid lines in the figure. For example, P N E is associate to the stencil composed by the cells Ω r,s for r ∈ {i, i + 1} and s ∈ {j, j + 1}. As in [14], we define the global smoothness indicator by where OSC[P ] is the multidimensional Jiang-Shu smoothness indicator, as defined in [17]. The nonlinear weights are computed by (3) starting from the linear Next we consider the case of a cell adjacent to a domain boundary. We focus in particular on the case of the bottom boundary, which is depicted in the central panel of Fig. 7. Here the reconstruction is CWENOZ(P opt ; P N E , P N W ,P E ,P W ), where P N E and P N W are defined as in the domain bulk. The stencil of P opt is biased towards the interior of the domain and is composed by the cells Ω r,s for r ∈ {i−1, i, i+1} and s ∈ {j, j+1, j+2}. The other two polynomials,P E andP W are degree 1 polynomials that depend only on the tangential variable, x in the example, and that are constant in the direction normal to the boundary. Their stencils are indicated with dashed lines in the figure. The global smoothness indicator τ for the cell in the example is copied from the cell Ω i,j+1 . The linear weights are similar to the bulk case, i.e. d 0 = 3 /4 and d N E = d N W = d E = d W = 1 /16. The case of the other boundaries is obtained from this one by symmetry.
Finally we describe the reconstruction in a domain corner, focusing on the case of the south-west one, which is represented in the right panel of Fig. 7.
Here, for stability purposes, we must include also a constant polynomial in the CWENOZ operator, denoted withP 0 ,in order to avoid spurious oscillations when a strong wave hits the corner.P 0 has of course the constant value coinciding with the cell average of the corner cell and its 1-cell stencil is represented by the filled circle in the picture. Following [30], we assign to the constant polynomial and to theP E ,P N polynomials an infinitesimal weight of d 0 = d N = d E = ∆x 2 and the reconstruction in the south-west corner cell is CWENOZ-AO(P opt ; N E ;P E ,P N ,P 0 ).
The stencil of P opt is again biased towards the interior of the domain and is composed by the cells Ω r,s for r ∈ {i, i + 1, i + 2} and s ∈ {j, j + 1, j + 2}.P E , similarly to the previous case, is a degree 1 polynomial that is constant in the y direction, whileP N is a degree 1 polynomial that is constant in the x direction. The global smoothness indicator τ for the cell in the example is copied from the cell Ω i+1,j+1 . The case of the other corners is obtained from this one by symmetry.
The polynomials associated to the stencils are computed as follows. Let S be a collection of neighbours of the cell Ω i,j that includes the cell itself and let Π ⊂ P d (x, y) be the subspace of the polynomials of degree d in two spatial variables where P S is sought. If the stencil S contains as many cells as the dimension of Π, the polynomial P S is the solution of the linear system composed by the equations P S r,s = u r,s for all (r, s) ∈ S, where the operator · r,s denotes the cell average of its argument over the cell Ω r,s . In the examples above, all polynomials with a tilde in their name are computed in this way.
When the cardinality of S is larger than dim(Π), we associate to S the solution of the following constrained least-squares problem: In the examples above, the polynomials P opt , P N E , P SE , P SW , P N W are computed in this way.
On Cartesian grids, the constrained least square problem can be easily turned into an unconstrained one by choosing a basis of Π consisting of a constant function and of polynomials orthogonal to the constant one. Explicit expressions for the coefficients of the polynomials in the domain interior can be found in [10].

Two-dimensional tests
The numerical scheme has been implemented with the help of the PETSc libraries [5,4] for grid management and parallel communications; the tests were run on a multi-core desktop machine equipped with an Intel Core i7-9700 processor and 64Gb of RAM. We show the results obtained with the Local Lax-Friedrichs numerical flux.
In all the tests we consider the two-dimensional Euler equations of gas dynamics: where ρ, u, v, p and E are the density, velocity in x and y direction, pressure and energy per unit mass. We consider the perfect gas equation of state E = p γ−1 + 1 2 ρ(u 2 + v 2 ) with γ = 1.4.

Convergence test
We compare the novel reconstruction with the one of [14] that makes use of ghost cells on the isentropic vortex test [31]. Of course there would be no need to use a ghost-less reconstruction with periodic boundary conditions, since it would be trivial to set up and fill in the ghost cells, but we conduct this as a stress-test to verify the order of accuracy of the novel reconstruction. The initial condition is characterized by a uniform ambient flow with constant temperature, density, velocity and pressure T ∞ = ρ ∞ = u ∞ = v ∞ = p ∞ = 1.0, onto which the following isentropic perturbations are added in velocity and temperature: where r = x 2 + y 2 and the strength of the vortex is set to β = 5.0. The computational domain is the square [−5, 5] 2 with periodic boundary conditions and the final time is set to t = 10 so that the final exact solution is the same as the initial state. We observe third order convergence rates in all variables (1-norm errors in density and energy are shown in Tab. 4). Compared with the CWENOZ3 scheme, the errors are no worse, and in some cases slightly better.

Two-dimensional Riemann problem
We have run a number of Riemann problems, in particular configurations B, G and K from [27], in order to compare the performances of the novel reconstruction on flows with waves (almost) orthogonal to the boundary. In our numerical experiments we have noticed that choosing correctly the linear weights for the low degree polynomials in the boundary cells is important to avoid spurious waves and features generated by an anomalous diffusion in the tangential direction; this latter shows up for example when choosing infinitesimal weights for the planes with two cells in the stencil in the middle panel of Fig. 7. Since it is on contact waves that spurious diffusion can accumulate over time, we report only a comparison of the solutions computed with the ghosted and the no-ghost reconstruction on configuration B of [27], which involves four contact discontinuities.
We  Fig. 8. In the plot the colors stand for pressure (rainbow colorbar) and we also showing contour lines of density (grayscale colorbar). We are focusing on contact waves as they are a good indicator of numerical diffusion, since on this kind of waves its effects accumulate over time. No difference is visible between the two computed solutions, indicating that the reconstruction that does not make use of ghost cells does not introduce significant differences with respect to the standard approach that makes use of ghosts. In particular, no wave deformation is visible close to the boundary, indicating that, with our choice of linear weights, no extra tangential diffusion is introduced in the boundary cells with respect to the cells

Radial Sod test
Next we run the cylindrical Sod shock tube problem in two space dimensions. The initial conditions for velocity component are u, v = 0 everywhere, while density and pressure are (ρ H , p H ) = (1, 1) for the central region, i.e. where r = x 2 + y 2 < 0.5, and (ρ L , p L ) = (0.125, 0.1) elsewhere. The computational domain is set to Ω = [0; 1] 2 with symmetry boundary conditions along x = 0 and y = 0 and wall boundary conditions on x = 1 and y = 1.
In Fig. 9 we compare the solutions at t = 0.2 computed on a grid of 400×400 cells, with and without using ghost cells. The solution is colored by pressure and 25 equispaced contour lines of the density field, from 0.04 to 1.0, are also shown (grayscale colorbar), so that the type of wave can be easily recognized. All solutions were computed in the first quadrant only, but the CWENOZ3 one is shown reflected to the left to ease the comparison.
Almost no difference can be appreciated between the two solutions and even the small artifacts that render the numerical solution non symmetrical appear identical in both schemes. Further, in Fig. 10, we plot the density computed without using ghost cells against the cell center distance from the origin. An almost perfect radial symmetry is observed, despite the fact that the boundary cells are reconstructed with a different algorithm than the bulk ones.
The solution of this problem after t = 0.2 sees the reflection of the cylindrical  shock wave with the outer walls and its later interaction with the expanding contact. The reflected curved shock interacts with itself exactly at the upperright corner at t 0.56; in Fig. 11, we show the solution at t = 0.6, just after this event. In this way we are testing the numerical schemes on reflecting a non planar shock wave on a wall. Even more importantly, we are stressing the reconstruction procedure in the corner cell, since the shock convergence happening there implies that for some timesteps there is no non-discontinuous stencil available to the reconstruction procedure for the corner cell. Here too, no appreciable difference is visible between the solutions computed with the two reconstruction schemes, showing that not using ghost cells in the reconstruction does not impair the numerical scheme.

Implosion problem
Next we consider the problem of a diamond-shaped converging shock proposed in [18]. It is computed on a quarter plane with symmetry boundary conditions, so that an oblique shock interacts with the boundary for a long time in the initial stages of the evolution and later the converging shock hits exactly the origin and is reflected back from there, testing again the non-oscillatory properties of the reconstruction in a corner cell when a strong wave impinges there. The test is set in the square domain [0, 0.3] 2 with reflective boundary conditions on all four sides: those at x = 0 and y = 0 represent symmetry lines, while the other two are physical solid walls. The initial condition has zero velocity everywhere and ρ = 1, p = 1 in the outer region (x + y > 0.15) and ρ = 0.125, p = 0.14 in the interior one. A useful reference for this test is [23] and the first author's website cited therein. We show the solution computed with a grid of 800 × 800 cells; the final time was set to t = 2.5, saving snapshots every 0.005 until t = 0.1 and every 0.1 afterwards. Fig. 12 shows both solutions in an early stage of the evolution, at t = 0.03. Here and in all subsequent figures, we have reflected to the left the solution computed with ghosts. In the early stages of the evolution the initial discontinuity gives rise to a shock (indicated with "S" in the left panel) and a contact ("C"), both moving towards the origin, and to a rarefaction ("R") that moves outwards. At the boundary, the shock is reflected and the reflected waves interact with the incoming contact ("s" and "c" in the figure). In the right panel, the gas velocity is represented with arrows; notice the fast wind directed towards the origin blowing along the coordinate axis.
Later the main shock and the reflected shocks converge in the origin, hit there head to head and are bounced back outwards. The snapshot reported in Fig. 13 is taken at t = 0.06, just after this event. Here it is important to observe that no spurious waves and no difference among the two schemes can be observed close to the origin, testifying that the reconstruction procedure in the corner cells is able to employ correctly the constant polynomial when the waves are very close to the corner.
The reflected contact, further deformed by the interaction with the expanding reflected shock, is being deformed by the wind blowing along the coordinate Figure 12: Implosion test at t = 0.03. The rainbow colorbar is for pressure, the grayscale one is for the density isolines. In the right panel, the arrows represent the velocity. In the left panel, the main shock (S), contact (C) and rarefaction (R) are indicated with capital letters, some secondary waves with small letters.  In Fig. 14 we show the solution at time t = 0.1. At this time the rarefaction is still moving outwards, while the rounded shock bounced back from the origin has overcome the incoming contact, which shows its physical instability coming from the deformations along the coordinate axis.
When the evolution is computed for long times, there is no consensus among the different schemes about the form of the bubbles near the origin and along the main diagonal [23]. In Fig. 15 we report the solutions computed for t = 2.5 (note the different pressure colorbar than in the previous ones). At this time, many waves reflections, refractions and interactions have taken place and the solution exhibit a quite complex pattern. The two schemes compute the main waves and the pressure field almost identically, but differ from each other in the bubbles, which appear to be physically unstable slip lines and therefore it is quite natural that different schemes can represent them in a different way.

Shock-bubble interaction
The last computation that we show is the shock-bubble interaction problem from [8]. Here a right-moving shock hits a standing bubble of gas at low pressure.  The shock, in its movement towards the right, sets in motion, compresses and deforms the bubble; it interacts with it and the resulting refracted shocks are bounced back towards y = 0 by the outer walls, giving rise to a very complex interaction pattern. The bubble is an unstable pattern and computing the solution with different schemes or different grid resolution will give it different final shapes. In Fig. 16 we show the solutions at the final time of the computation, t = 0.4; the CWENOZ3 solution is reflected along the symmetry axis in order to ease the comparison with the no-ghost solution computed with CWZb3. Here the main difference can be seen in the portion of the bubble that remains attached to the symmetry plane. All the other waves, including the central portion of the bubble, are almost identical in both computations.

Conclusions and perspectives
In this paper we have pursued further the approach of [26] for reconstructions in finite volume schemes without using ghost cells. While the main motivation there was the application in internal nodes of networks, where it is difficult to prescribe appropriate extrapolations, here we have focused on the accuracy of the reconstructions and on the extension to higher dimensions. All the proposed reconstructions in the boundary cells are based on reconstruction stencils that are not symmetric (but extended only inwards) and nevertheless the optimal accuracy on smooth flows can be achieved.
Regarding the first point, we have proposed the employment of Z-weights instead of the Jiang-Shu nonlinear weights, obtaining a reconstruction that has less numerical diffusion than the one of [26]. Perhaps more importantly, the novel reconstruction can reach the optimal third order of convergence in a more ample subset of the parameter space and in particular can employ a very small without needing the infinitesimal weight of the constant polynomial to be O(∆x 2 ).
In the second part of the paper, we have proposed a (non dimensionally-split) reconstruction for two-dimensional Cartesian grids. This has been compared with the CWENOZ approach of [14] that uses ghost cells, showing that the employment of ghost cells can be entirely avoided without affecting the quality of the computed solutions.
The approach without using ghosts appears thus to be quite promising since setting them is always a tricky point in numerical schemes; it seems interesting to pursue further this line of research towards higher accuracy or for moving boundaries like in piston problems.