Rémi Douvenot
email: remi.douvenot@recherche.enac.fr

Vincent Fabbro

Kevin Elis

Yvonick Hurtaud

Comparison of Metrics for Clutter Data Comparison

Keywords: refractivity from clutter, propagation, metrics, clutter

Some radar applications require to measure a distance between different sets of clutter data. For instance, refractivity from clutter (RFC) is a technique to retrieve the atmospheric conditions from the measured clutter in the absence of target. The clutter measurements are compared to simulations to infer the refractive index along the pathway. When a distance between two clutter data is required, the L2 norm between data in dB is usually chosen. This paper investigates other metrics and claims that the L2 norm is not the most appropriate.

I. INTRODUCTION

Refractivity from clutter (RFC) [START_REF] Rogers | Likelihood estimation of tropospheric duct parameters from horizontal propagation measurements[END_REF], [START_REF] Krolik | Tropospheric refractivity estimation using radar clutter from the sea surface[END_REF] is a technique to retrieve the low troposphere atmospheric conditions. The basic idea is to infer the propagation channel from the large-scale variations in the clutter measurements. These measurements are performed in a monostatic configuration, the radar illuminating the sea surface. Atmospheric ducts significantly modify the radar coverage. The aim of RFC is to infer the atmospheric ducts to predict the coverage of electromagnetic systems in the presence of the ducts.

To do so, inverse methods are required. In any inverse method, an objective function has to be defined. In the RFC literature, authors use the L2 norm as the objective function [START_REF] Karimian | Refractivity estimation from sea clutter: An invited review[END_REF]. So, the inverse problem is to retrieve the duct described by the vector parameter m that minimises

P dB o -P dB s (m) 2 2 , (1)
where P dB o is the observed clutter power with respect to the range expressed in dB, and P dB s is the power clutter with respect to the range expressed in dB and simulated in the presence of an atmosphere described by the parameters m.

• 2 is the L2 norm, defined for a vector X of length N by

X 2 = N n=1 X 2 n . (2)

II. PROPOSED METRICS

Recently, Douvenot et al. [START_REF] Douvenot | Parameter-based rules for the definition of detectable ducts for an RFC system[END_REF] have shown that taking values in dB in the objective functions overestimates the difference between two low power clutter data. Moreover, clutter data are functional data, and the L2 norm does not take this information into account. From this fact, two metrics are proposed in this paper, namely φ L1 and φ L2 .

The metrics φ L1 and φ L2 are derived from the L1 and L2 norm, respectively. Considering that clutter data are continuous, the covariance is introduced to define the metrics. φ L1 and φ L2 are defined as

φ L1 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 n 2 P o -P s (m) 2 1 cov(P o , P s ) if cov(P o , P s ) > 0, +∞ if cov(P o , P s ) ≤ 0, (3)
and

φ L2 = ⎧ ⎪ ⎨ ⎪ ⎩ 1 n 2 P o -P s (m) 2 2 cov(P o , P s ) if cov(P o , P s ) > 0, +∞ if cov(P o , P s ) ≤ 0. ( 4
)
where cov is the covariance function, n is the length of the vectors P o/s , and • 1 is the L1 norm, defined for a vector X of length N by

X 1 = N n=1 |X n | . (5)
Note that φ L1 and φ L2 are dimensionless. Compared to the usual L2 norm, the major change is the weighting by the covariance. Whereas L2 and L1 norms are point-to-point norms, the covariance gives an indication on the overall behaviour of the data [START_REF] Ferraty | NonParametric Functional Data Analysis. Theory and Practice, ser. Springer Series in Statistics[END_REF]. This is relevant since propagation factors are handled as nonparametric functional data. This improvement dismisses some major drawbacks and improves the metrics quality in general, as shown in section III.

For the metric φ L1 , the L1 norm is chosen. Indeed, L1 norm favours the "sparsity" of P o -P s (m) when the L2 norm favours a minimised mean error. Indeed, the simulated clutter data usually have minima more pronounced than the measured ones. This large local discrepancies affect less the L1 norm than the L2 norm.

Note that φ L1 and φ L2 are semimetrics because they do not satisfy the triangle inequality. However, for the sake of conciseness, we use the term metric to mention them. In inverse and optimisation methods, the triangle inequality is not necessary. Thus semimetrics are compatible with RFC requirements. 3) and ( 4). This study is based on 1998 Wallops Island data [START_REF] Rogers | Estimating evaporation duct heights from radar sea echo[END_REF] measured with the SPANDAR radar in the presence of surface-based ducts.

c 0 c 0 c 0 c 0 c 2 c 2 z b M-index Altitude z t z t δz b z b + z t

A. Principle of the comparison

A library of N propagation factors corresponding to N ducting conditions is generated by simulations. These simulations are performed by a parabolic equation method [START_REF] Levy | Parabolic Equation Methods for Electromagnetic Wave Propagation[END_REF]. The ducts are considered as trilinear with a simple distance variation. They are described by the parameters z b the duct height, c 2 the M-index slope in the inversion layer, and z t the inversion layer thickness. Finally, δz b is the inversion layer height variation with the distance, see figure 1 From these two ducts, two propagation factors are obtained with respect to altitude and distance, as plotted in figure 3. The propagation factor just above the sea level (figure 4) is extracted because in it proportional to the sea clutter received power, as detailed in section III-B. RFC consists in inferring the duct conditions from the clutter return with respect to distance.

RFC is then performed on two examples from Wallops measurement data. On these examples, the measured propagation factor and the simulated propagation factors that minimise the metrics are shown. The following norms are compared: L2 squared (2), φ L2 (4), L1 squared [START_REF] Ferraty | NonParametric Functional Data Analysis. Theory and Practice, ser. Springer Series in Statistics[END_REF], and φ L1 [START_REF] Karimian | Refractivity estimation from sea clutter: An invited review[END_REF].

Note that in this study, the aim is to compare the metrics by comparing the nearest neighbour with respect to each metric: the ducts responsible for the propagation factors are not considered.

B. The measured data

The data were measured on the Wallops Island, Virginia, on April 1998. They were collected by Space Range Radar (SPANDAR) in coastal environment. The radar operates at 2.84 GHz in vertical polarisation. The antenna is 30.78 m high, with an horizontal beamwidth of 0.4 • . Clutter to noise ratios (CNR) were measured with respect to range and azimuth [START_REF] Rogers | Estimating evaporation duct heights from radar sea echo[END_REF].

From the radar equation, the CNR measured with respect to the distance R on one azimuth is expressed as [START_REF] Douvenot | On the knowledge of radar coverage at sea using real time refractivity from clutter[END_REF]

CNR = P e G 2 λ 2 (4πR) 3 C 0 σ 0 r F 4 , (6)
where P e is the emitted power, G is the antenna gain, λ is the wavelength, F is the propagation factor, C 0 is a constant, and σ 0 r is the reduced normalised radar cross section considered as constant. It is defined such that σ 0 r θ 4 g = σ 0 , σ 0 being the normalised radar cross section, and θ g being the grazing angle. σ 0 r is considered as constant. Note that other interpretations lead to the same distance variation of the CNR [START_REF] Karimian | Refractivity estimation from sea clutter: An invited review[END_REF].

The propagation factor F 2 is deduced from the clutter return up to a constant multiplier. This constant is suppressed by normalising the signal at 10 km, obtaining a normalised propagation factor F 2 n (R). What is called the measured propagation factor in the following of the text is this normalised propagation factor.

Finally, the data are smoothed on 1 km in distance and 5 • in azimuth.

C. Results

The results presented in this section are not typical of all the data from 1998 Wallops Island campaign. Sometimes, the different metrics all lead to the same simulated propagation factor (when the library contains a simulation very similar to the measurement). The presented examples are chosen to highlight the behaviours specific to the different metrics.

A first result is presented in figure 5. It corresponds to a measurement carried out on April 2, 1998, at 17h20'47", on the azimuth 170 • . The propagation factors are plotted with respect to the distance R. The measured one is in black. In dotted red, continuous red, dotted blue, and continuous blue are the propagation factors from the library the closest to the measured one with respect to the L2 norm, φ L2 , the L1 norm, and φ L1 , respectively.

The L2 norm propagation factor tends to follow all the maxima and minima of the measured propagation factor. However, none of the data in the library being close to the measurement, the proposed best one is an unsatisfying compromise that misses a maximum at 40 km, and underestimates all the extreme values.

The L1 norm propagation factor follows the measurements on the first 30 km to the detriment of the rest of the fitting. Indeed, it favours sparseness (the objective function is almost equal to zero on these first 30 km).

When the covariance weighting is added to the L1 and L2 norms, the result is close to the measurement (the continuous blue and red curves coincide). The overall shape of the measured propagation factor is correctly described by the chosen simulation. The correlation emphasises the propagation factors that present the same minima/maxima than the measurement. Note that in this case, the metrics φ L1 and φ L2 give the same result.

A second result is presented in figure 6. It corresponds to a measurement carried out on April 2, 1998, at 17h10'25", on the azimuth 100 • . The L1 and L2 norms give poor results since their best solutions miss all the extrema of the measurement. This case emphasises on the need for a metric that takes into account the global behaviour of the data. The metrics φ L1 and φ L2 give results of similar shapes. The peak at 45 km has less influence with the L1 norm, that is why φ L1 result seems a bit better. From the presented results and the global observed behaviours, φ L1 and φ L2 metrics are much more appropriate than L1 and L2 norms for quantifying the similarity between two propagation factors obtained from clutter data. However, the debate between the choice of φ L1 or φ L2 should stay open without further studies.

IV. CONCLUSION

In this paper, new metrics including the covariance have been introduced. They are used to measure the distance between two propagation factors. On the presented results, the proposed metrics are more efficient than the usual L2 norm. However, further studies should be necessary to take a decision. Moreover, the signal after the radioelectric horizon (here around 26.4 km) is more important for this application since the effect of atmospheric ducts is prominent. It could be taken into account by modifying the proposed metrics as in [START_REF] Douvenot | Retrieve the evaporation duct height by least-squares support vector machine algorithm[END_REF].

This point is important for the RFC application because this latter requires inverse or optimisation algorithms that necessitate an objective function. The accuracy and the convergence speed of these algorithms depend on the objective function. Thus this latter should be carefully chosen. Note that the proposed functions are semimetrics that are compatible with the needs of inverse or optimisation algorithms.

Fig. 1 .

 1 Fig. 1. The trilinear duct, (left) at the antenna, (right) at the maximum range. δz b is the height variation of the duct.

  . The parameters are chosen within the following limits: 0 ≤ z b ≤ 150 m, 0 ≤ -c 2 ≤ 1.5 M-units.m -1 , 0 ≤ z t ≤ 70 m, and -20 ≤ δz b ≤ 20 m. N = 10000 ducts are generated with a Latin hypercube sampling for the library to statistically cover all the possible ducts [8]. From these latter, N propagation factors are generated with the parabolic equation software simulating the conditions of the 1998 Wallops Island measurement campaign. To highlight how the parameters and the clutter data are related, two samples are presented. The two ducts considered are the standard atmosphere (z b = 0 m, c 2 = 0 M-unit.m -1 , z t = 0 m, and δz b = 0 m) and a surface-based duct (SBD) of parameters z b = 10 m, c 2 = 0.25 M-unit.m -1 , z t = 40 m, and δz b = 0 m. They are represented with respect to altitude in figure 2. (They are invariant with the distance.)

Fig. 2 .

 2 Fig.2. The modified refractivity of the two samples corresponding to the standard atmosphere (Std) and to a surface-based duct (SBD).

AltitudeFig. 3 .

 3 Fig. 3. Propagation factors (in dB) obtained with respect to altitude and distance from the antenna for the two considered atmospheric conditions.

Fig. 4 .

 4 Fig. 4. Propagation factors with respect to distance just above the sea level corresponding to the standard atmosphere (Std) and to a surface-based duct (SBD).

Fig. 5 .

 5 Fig. 5. First example. The simulated propagation factors from the database the nearest to the measurement (black) with respect to the L2 norm (dotted red), φ L2 (continuous red), L1 norm (dotted blue), and φ L1 (continuous blue). The continuous blue and red curves coincide.

Fig. 6 .

 6 Fig.6. Second example. The simulated propagation factors from the database the closest to the one from measurement (black) with respect to the L2 norm (dotted red), φ L2 (continuous red), L1 norm (dotted blue), and φ L1 (continuous blue).

ACKNOWLEDGMENT

The presented work has been carried out in the framework of the French-German Technical Agreement 119 : "Modelling and measurements of sea clutter at centimetre and millimetre bands; application to atmospheric refractivity assessment" involving FHR, WTD71, DGA-MI, and ONERA. The authors would like to thank the French MoD for its financial support.


Pierre Cantin
email: pircantin@gmail.com

Alexandre Ern
email: alexandre.ern@enpc.fr

Ern

Vertex-Based Compatible Discrete Operator Schemes on Polyhedral Meshes for Advection-Diffusion Equations

Keywords: Polyhedral meshes, compatible discretization, advection, diffusion, Peclet robustness, divergencefree velocity AMS Subject Classification. 65N15, 65N12, 76Rxx

We devise and analyze vertex-based, Péclet-robust, lowest-order schemes for advection-diffusion equations that support polyhedral meshes. The schemes are formulated using Compatible Discrete Operators (CDO), namely primal and dual discrete differential operators, a discrete contraction operator for advection, and a discrete Hodge operator for diffusion. Moreover, discrete boundary operators are devised to weakly enforce Dirichlet boundary conditions. The analysis sheds new light on the theory of Friedrichs' operators at the purely algebraic level. Moreover, an extension of the stability analysis hinging on infsup conditions is presented to incorporate divergence-free velocity fields under some assumptions. Error bounds and convergence rates for smooth solutions are derived, and numerical results are presented on three-dimensional polyhedral meshes.

Introduction

The goal of this work is to approximate the scalar-valued function p : Ω → R solving the following advectiondiffusion problem:

-∇•(λ∇p) + β•∇p = s a.e. in Ω, (1.1a) p = p D a.e. on ∂Ω, (1.1b) where Ω is a bounded polyhedral connected subset of R 3 with boundary ∂Ω and outward unit normal n, λ a bounded, symmetric, uniformly positive-definite tensor-valued field in Ω, β a vector-valued field in W 1,∞ (Ω), s ∈ L 2 (Ω), and p D ∈ H t (∂Ω), t > 1. We use boldface fonts for vector-and tensor-valued quantities. In addition to the classical assumption on the sign of ∇•β, we also include in our analysis an extension to the case of divergence-free advection; see below. This extension is by no means straightforward and is rarely addressed in the literature. We also briefly discuss the (simpler) variants where the advection term is written in divergence form and where there is a zero-order reactive term. Of particular interest is the boundary conditions. Indeed, as for stabilized finite element and dG methods, weak enforcement of Dirichlet conditions yields better results for under-resolved outflow layers. To this purpose, we extend Nitsche's boundary penalty method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] to the CDO setting.

Let us put our work in perspective with existing schemes. For pure advection, we emphasize that the CDO scheme is essentially an upwind finite volume scheme on a dual mesh with vertex-based control volumes. Thus, the analysis uses similar techniques to those used for dG methods; see Johnson and Pitkäranta [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF], Brezzi et al. [START_REF] Brezzi | Discontinuous Galerkin methods for first-order hyperbolic problems[END_REF], and [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]. We also mention the following recent approaches to discretize pure advection equations in the setting of differential geometry. Using the notion of extrusion defined by

Bossavit in [START_REF] Bossavit | Extrusion, contraction: their discretization via Whitney forms[END_REF], Heumann and Hiptmair [START_REF] Heumann | Extrusion contraction upwind schemes for convection-diffusion problems[END_REF] and Mullen et al. [START_REF] Mullen | Discrete Lie advection of differential forms[END_REF] proposed a discretization of interior products respectively on triangular and Cartesian meshes. Stabilized Galerkin methods for differential forms are considered by Heumann and Hiptmair [START_REF] Heumann | Eulerian and semi-Lagrangian methods for convection-diffusion for differential forms[END_REF]. Palha et al. [START_REF] Palha | Mimetic spectral element advection[END_REF] proposed another approach using the wedge product as the adjoint operator of the interior product. Furthermore, for advection-diffusion, the present CDO scheme is, to our knowledge, the first polyhedral discretization that is only vertex-based and that is robust for dominant advection up to the limit of zero diffusion. Another framework for vertex-based polyhedral schemes for elliptic PDEs is that of Virtual Element Methods (VEM), see Beirão da Veiga et al. [START_REF] Beirão Da Veiga | Basic principles of virtual element methods[END_REF].

The difference is that we use explicit reconstruction functions (typically piecewise constant on subcells) and we treat dominant advection, but our schemes are only of lowest-order. A first alternative to vertex-based schemes are face-based schemes: an arbitrary-order, Péclet-robust, face-based scheme for advection-diffusion has been recently analyzed by Di Pietro et al. [START_REF] Di Pietro | A discontinuous-skeletal method for advection-diffusionreaction on general meshes[END_REF] (see also Beirão da Veiga et al. [START_REF] Beirão Da Veiga | A unified approach for handling convection terms in finite volumes and mimetic discretization methods for elliptic problems[END_REF] for the lowest-order and diffusion-dominated case). Another alternative is to use a cell-based dG method, but the treatment of diffusion requires introducing interior penalty parameters and using order k ≥ 1 in cells leading to an increase of DoFs.

The material is organized as follows. In Section 2, we introduce the main notation for the discrete setting. In Section 3, we devise and analyze CDO schemes for pure advection. In Section 4, we treat advection-diffusion. Both sections 3 and 4 focus on assumption (β1) for the velocity field for simplicity. In Section 5, we revisit the analysis in the case of a divergence-free velocity field under assumption (β2). In Section 6, we present numerical results on three-dimensional polyhedral meshes. Finally, in Section 7, we collect some proofs of technical results.

CDO schemes are formulated by considering a dual mesh M := { V, E, F, C} such that there is a one-toone pairing between primal vertices and dual cells, primal edges and dual faces, and so on (see Figure 1).

In particular, f (e) denotes the dual face associated with the primal edge e ∈ E, and c(v) the dual cell associated with the primal vertex v ∈ V. Dual mesh entities are oriented by the associated primal entity; for instance, n f (e) is the unit normal vector to f (e) oriented by t e . There are many possibilities to build a dual mesh. In this work, we assume that primal faces are planar and star-shaped w.r.t. their barycenter and that primal cells are star-shaped w.r.t. their barycenter, and we consider the fully barycentric dual mesh built using barycenters of all the primal mesh entities.

c(v) • v f (e) n f (e) e te n f f ẽ(f ) t ẽ(f ) t ẽ(f ) • ṽ(c 1 )
• ṽ(c 2 )

• ṽ(c 3 )

• ṽ(c 4 ) c(v) In what follows, we assume that the meshes M and M satisfy a regularity requirement stating that there exists a common simplicial sub-complex of M and M (i.e., any k-simplex, 1 ≤ k ≤ 3, in this sub-complex belongs to only one k-cell of M and of M) such that all the k-simplices are shape-regular in the usual sense of Ciarlet and any k-cell of M or M contains a uniformly bounded number of k-simplices. This mesh regularity assumption is only needed to analyze the schemes, but not for implementation. For any primal or dual mesh entity x, h x denotes the diameter of x; moreover, when deriving convergence rates for smooth solutions, we use h to denote the largest primal cell diameter. To alleviate the notation, we abbreviate A B the inequality A ≤ cB with positive constant c whose value can change at each occurrence as long as it is uniform with respect to the mesh-size and the model parameters.

• v f (e)
Since boundary conditions are enforced weakly in this work, we consider mesh entities at the boundary.

The trace of the primal mesh M at the boundary ∂Ω defines a cellular complex M ∂ := {V ∂ , E ∂ , F ∂ }, where V ∂ collects all the primal vertices lying at the boundary, and so on. Instead, the dual mesh has no entities lying at the boundary, so that we introduce an additional set of dual faces

F ∂ = { f ∂ (v) | v ∈ V ∂ } with f ∂ (v) := ∂c(v) ∩ ∂Ω; observe the one-to-one pairing between V ∂ and F ∂ .

Degrees of freedom

The degrees of freedom (DoFs) of discrete fields are attached to mesh entities according to their physical nature. For instance, the degrees of freedom of a discrete potential field (0-cochain) are attached to vertices, either primal or dual ones. In this work, we focus on vertex-based CDO schemes where these DoFs are attached to primal vertices. For a discrete potential q, we use the notation q ∈ V ≡ R #(V) , where V is the vector space composed of DoFs attached to primal vertices and #(V) denotes the cardinality of the set V.

We write q v for the value of q attached to the vertex v ∈ V. We also consider discrete circulation fields (1cochains) in E (attached to primal edges), discrete flux fields (2-cochains) in F (attached to dual faces), and discrete density fields (3-cochains) in C (attached to dual cells). Owing to the one-to-one pairing between primal and dual mesh entities, the vector spaces V and C are isomorphic, as well as the vector spaces E and F. This leads us to define the algebraic duality products q, s VC = v∈V q v s c(v) for all (q, s) ∈ V × C, and g, φ E F = e∈E g e φ f (e) for all (g, φ) ∈ E × F. To weakly enforce boundary conditions, we consider discrete fields at the boundary, and with obvious notation, we introduce the isomorphic vector spaces V ∂ and F ∂ , along with the algebraic duality product

q ∂ , φ ∂ (VF) ∂ = v∈V ∂ q ∂ v φ ∂ f ∂ (v)
. Furthermore, for all q ∈ V, we use the notation q ∂ ∈ V ∂ with q ∂ v = q v for all v ∈ V ∂ , i.e., q ∂ collects the DoFs of q attached to boundary vertices.

To measure the approximation error resulting from CDO schemes and to discretize the source terms and boundary conditions, we define DoFs for continuous fields. One possibility is to consider the classical de Rham maps (we also consider other choices below) for smooth enough fields. In what follows, we consider

the maps R V : S V (Ω) → V, R E : S E (Ω) → E, R F : S F (Ω) → F and R C : S C (Ω) → C such that (R V (p)) v = p(v) for all v ∈ V, (R E (g)) e = e t e •g de for all e ∈ E, (R F (φ)) f = f n f •φ d f for all f ∈ F, (R C (s)) c = c s dc for all c ∈ C.
Possible choices for the domains of the de Rham maps are

S V (Ω) = H s (Ω) with s > 3 2 , S E (Ω) = H s (Ω) with s > 1 or S E (Ω) = {g ∈ L p (Ω), ∇×g ∈ L t (Ω)} with p > 2 and t > 1 2 , and S F (Ω) = H s (Ω) with s > 1 2 or S F (Ω) = {φ ∈ L p (Ω), ∇•φ ∈ L 2 (Ω)} with p > 2 and S C = L 1 (Ω).
At the boundary, we use the maps

R V ∂ : S V ∂ (∂Ω) → V ∂ and R F ∂ : S F ∂ (∂Ω) → F ∂ such that (R V ∂ (p)) v = p(v) for all v ∈ V ∂ , (R F ∂ (φ)) f ∂ = f ∂ φ d f for all f ∂ ∈ F ∂ , with S V ∂ (∂Ω) = H t (∂Ω), t > 1, and S F ∂ (∂Ω) = L 1 (∂Ω).

Discrete differential operators

For all v ∈ V and all e ∈ E, we set ι v,e = 1 if v is the extremity of e toward which t e points, ι v,e = -1 if v is the other extremity of e, and ι v,e = 0 if v is not an extremity of e. The discrete gradient operator GRAD : V → E is defined such that (GRAD(q)) e = v∈V ι v,e q v for all q ∈ V; note that the algebraic representation of GRAD is a rectangular matrix with entries in {0, ±1}. We also define a discrete dual divergence operator DIV :

F → C such that ( DIV(φ)) c(v) = f (e)∈ F ι f (e),c(v) φ f (e) for all φ ∈ F, with ι f (e),c(v) = -ι v,e .
Observe that ι f (e),c(v) = 1 (resp., -1) if f (e) is a face of the dual cell c(v) such that n f (e) points outward (resp., inward) c(v), and ι f (e),c(v) = 0 if f (e) is not a face of c(v).

The following discrete adjunction property holds:

GRAD(q), φ E F = -q, DIV(φ) VC , ∀(q, φ) ∈ V × F. (2.1)
Other important properties are the two following commuting properties with the de Rham maps:

GRAD(R V (p)) = R E (∇p), ∀p ∈ S V (Ω), (2.2a) q, R C (∇•φ) VC = q, DIV(R F (φ)) VC + q ∂ , R F ∂ (n•φ) (VF) ∂ , ∀φ ∈ S F (Ω), ∀q ∈ V. ( 2

Restriction to primal cells and boundary faces

It

) f = f (n f •φ) d f for all f ∈ F c .
The local discrete gradient operator GRAD c : V c → E c is defined similarly to GRAD. We also define the following local norms:

|||q||| 2 2,Vc := h 3 c v∈Vc q 2 v , |||g||| 2 2,Ec := h c e∈Ec g 2 e , (2.3)
for all q ∈ V c and all g ∈ E c . The global counterparts of these norms are assembled cell-wise as |||q||| 2 2,V := c∈C |||q||| 2 2,Vc and |||g||| 2 2,E := c∈C |||g||| 2 2,Ec , for all q ∈ V and all g ∈ E. Let now f ∈ F ∂ be a primal boundary face. We define the local subset

V ∂ f = {v ∈ V ∂ | v ∈ ∂f }, collecting boundary vertices attached to f . For all v ∈ V ∂ f , we define f ∂ f (v) := f ∂ (v) ∩ f
f := { f ∂ f (v); v ∈ V ∂ f }. The vector space V ∂ f is composed of the DoFs of discrete boundary potential fields q ∂ attached to V ∂ f ; similarly for F ∂ f . The de Rham map R F ∂ f is such that (R F ∂ f (φ)) f ∂ = f ∂ φ d f ∂ for all f ∂ ∈ F ∂ f . • • • • e e Õ c fc(e) d(e Õ ) fl c 9 Def ⁄ ◊ e • v • w • f (e) • • • n f (e),c(v) n f (e),c(w)
• • v • • • • f f (v) f 9 Def ⁄ ◊ e • v • w • f (e) • • • n f (e),c(v) n f (e),c(w)
(β•n) ± = 1 2 (|β•n| ± β•n) ≥ 0. We introduce the graph space V β (Ω) = {q ∈ L 2 (Ω); β•∇q ∈ L 2 (
Ω)}; functions in the graph space have a trace in L 2 (|β•n|; ∂Ω) provided ∂Ω -and ∂Ω + are well separated (see [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF]). In this context, a well-posed weak formulation of problem (3.1) (see [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF][START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]) is as follows: Find p ∈ V β (Ω) such that

a β (p, q) = Ω qs dΩ + ∂Ω qφ D d∂Ω, ( 3.2)
for all q ∈ V β (Ω), with boundary flux φ D = (β•n) -p D , and with bilinear form

a β (p, q) = Ω q(β•∇p) dΩ + ∂Ω q(β•n) -p d∂Ω. ( 3.3)
Note that the boundary integrals vanish outside ∂Ω -.

CDO scheme

Vertex-based CDO schemes for pure advection are built using two discrete operators: a discrete contraction operator I E C β : E → C, which is the discrete counterpart of the map g → β•g, and a discrete boundary Hodge operator H ∂ α : V ∂ → F ∂ (indexed by a surface function α ∈ L ∞ (∂Ω)) which is the discrete counterpart of the map p → αp at the boundary. Using these operators, the following discrete problem can be formulated:

Find p ∈ V such that A β (p, q) = q, s VC + q ∂ , φ D (VF) ∂ , ∀q ∈ V, (3.4)
with bilinear form such that In the spirit of Friedrichs operators [START_REF] Ern | Discontinuous Galerkin methods for Friedrichs' systems. I. General theory[END_REF][START_REF] Ern | An intrinsic criterion for the bijectivity of Hilbert operators related to Friedrichs' systems[END_REF], we assume that there is a second discrete contraction operator I V F β , which is the discrete counterpart of the map p → βp, and such that the following two properties hold: (I1) [Discrete Leibniz rule] The bilinear map on V × V such that

A β (p, q) := q, I E C β (GRAD(p)) VC + q ∂ , H ∂ (β•n) -(p ∂ ) (VF) ∂ , ( 3
p, H VC -∇•β (q) VC := p, I E C β (GRAD(q)) VC -p, DIV(I V F β (q)) VC -p ∂ , H ∂ β•n (q ∂ ) (VF) ∂ , ( 3.6)
is symmetric and satisfies p, H VC -∇•β (q) VC ≥ (ess inf Ω -∇•β)|||q||| 2 2,V for all q ∈ V, so that H VC -∇•β is monotone under assumption (β1); see Remark 3.1 below for the boundary term in (3.6).

p oe V g oe E oe  F s oe  C p ˆoe V  F ˆGRAD Á DIV I E  C - H E F ⁄ H (-•n) ≠ RHS p oe V g oe E oe  F s oe  C p ˆoe V  F ˆGRAD Á DIV I V  F - H E F ⁄ H (-•n) ≠ RHS p oe V g oe E oe  F s oe  C p ˆoe V  F ˆGRAD Á DIV H E F ⁄ ' H ⁄/h N ⁄ RHS p oe V g oe E oe  F s oe  C p ˆoe V  F ˆGRAD Á DIV P V H E F ⁄ ' H ⁄/h N ⁄ RHS 2
V × V such that p, q upw,β := p, I E C β (GRAD(q)) VC + q, DIV(I V F β (p)) VC , ( 3.7)
defines a semi-inner product.

Concerning the discrete boundary Hodge operator, we assume the following:

(H ∂ ) For all α ∈ L ∞ (∂Ω), H ∂ α is self-adjoint,
α ≥ 0, q ∂ , H ∂ α (p ∂ ) (VF) ∂ defines a semi-inner product on V ∂ × V ∂ .
counterpart of Ω p(β•∇q)dΩ + Ω q∇•(βp)dΩ -∂Ω q(β•n)pd∂Ω = 0.
At the discrete level, this quantity can be non-zero owing to the use of stabilization. We also notice that the symmetry of the map p, q upw,β results from p, q upw,β -q, p upw,β = p, H VC -∇•β (q) VC -q, H VC -∇•β (p) VC = 0 where we have used the selfadjointness of H ∂ β•n and of H VC -∇•β . Finally, we observe that I E C β does not, in general, depend linearly on its argument β owing to the use of stabilization.

Remark 3.2 (Conservative advection).

A possible variant of (3.1) is to consider the conservative form of the advective derivative. The PDE becomes ∇•(βp) = s in Ω, and a Dirichlet boundary condition can still be enforced at the inflow boundary. Assumption (β1) is then modified as follows: There exists a real number

τ > 0 such that ∇•β ≥ τ -1 a.e.

in Ω. The discrete bilinear form then becomes

A β (p, q) = q, DIV(I V F β (p)) VC + q ∂ , H ∂ (β•n) + (p ∂ ) (VF) ∂ . (3.8)
The design of the discrete contraction and boundary Hodge operators still hinges on (I1)-(I2) and (H ∂ ).

Remark 3.3 (Reaction). Another possible variant is to include a zero-order reaction term in the PDE which becomes β•∇p + µp = s in Ω with Lipschitz reaction coefficient µ (the conservative form of the advective derivative can also be considered). Then, the reaction-related bilinear form A µ (p, q) = v∈V µ v p v q v is added to the discrete problem, where µ v denotes (for instance) the mean-value of µ in c(v).

Example: CDO scheme with upwinding

Let us give a concrete example for the CDO scheme (3.4). We introduce the notation

β e := (R F (β)) f (e) = f (e) β•n f (e) d f , ∀e ∈ E. (3.9)
We also set The reason to distinguish the properties β e Λ e ≥ 0 in (Λ1) and β e Λ e ≥ c Λ |β e | in (Λ2) is that the former is satisfied by the so-called centered scheme corresponding to Λ ve = 0 for all v ∈ V e , and the latter by an upwind scheme. Classical upwinding corresponds to the choice Λ ve = sign(ι f (e),c(v) β e ) (with sign function

E v := {e ∈ E | v ∈ e} for all v ∈ V,
sign(t) = -1 if t ∈ R <0 , sign(0) = 0, and sign(t) = 1 if t ∈ R >0 ), so that (Λ2) holds with c Λ = 1. With
this choice, the solution delivered by the CDO scheme coincides with that of the upwind FV scheme on the dual mesh.

The discrete contraction operator I E C

β : E → C is defined such that, for all g ∈ E,

(I E C β (g)) c(v) := e∈Ev g e 1 2 (1 -Λ ve )β e , ∀v ∈ V, (3.10)
while the companion operator I V F β : V → F is defined such that, for all q ∈ V,

(I V F β (q)) f (e) := v∈Ve q v 1 2 (1 + Λ ve )β e , ∀e ∈ E. (3.11)
Moreover, the discrete boundary Hodge operator

H ∂ α : V ∂ → F ∂ with α ∈ L ∞ (∂Ω) is defined such that, for all q ∂ ∈ V ∂ , (H ∂ α (q ∂ )) f ∂ (v) := q ∂ v f ∂ (v) α d f , ∀v ∈ V ∂ . (3.12)
Observe that H ∂ α is algebraically represented by a diagonal matrix.

Remark 3.4 (Upwinding design).

There are several possible variations in the geometric quantities considered for upwinding. Instead of considering the full dual face f (e) as in (3.9), one possibility is to consider the average of the normal advection velocity on the dual sub-faces fc (e), and to design the upwinding parameters based on the sign of these quantities. In general, the smaller the underlying geometric objects, the larger the dissipation introduced by upwinding. The advantage of considering the dual sub-faces fc (e) is that upwinding is then compatible with the assembly of the scheme on primal cells.

12). Assume (Λ1). Then, (I1)-(I2) hold with bilinear maps

p, H VC -∇•β (q) VC = v∈V p v q v c(v) -∇•β dc, (3.13a) p, q upw,β = e∈E [[p]] e [[q]] e Λ e β e , (3.13b)
so that H VC -∇•β is algebraically represented by a diagonal matrix, and where [[q]] e = v∈Ve ι f (e),c(v) q v is the jump of q across the dual face f (e) for all q ∈ V. Remark 3.5 (First-order diffusion). Since [[q]] e = -(GRAD(q)) e , the right-hand side of (3.13b) corresponds to the jump penalty term considered in dG methods; for k = 0, it can be interpreted as adding a first-order viscosity term.

Proof. Proof of (3.13a). Let v ∈ V and let p, q ∈ V. Using (3.10)-(3.11), we infer that

p, I E C β (GRAD(q)) VC = v∈V e∈Ev v ∈Ve p v q v ι v ,e 1 2 (1 -Λ ve )β e , p, DIV(I V F β (q)) VC = v∈V e∈Ev v ∈Ve p v q v ι f (e),c(v) 1 2 (1 + Λ v e )β e .
Using ι v,e = -ι f (e),c(v) , the definition of β e , and (Λ1) leads to

p, I E C β (GRAD(q)) VC -p, DIV(I V F β (q)) VC = v∈V p v q v e∈Ev -f (e) β•n f (e),c(v) d f . To conclude, we observe that if v ∈ V \ V ∂ , e∈Ev f (e) β•n f (e),c(v) d f = ∂c(v) β•n c(v) d∂c = c(v) ∇•β dc
owing to the divergence theorem, while for the boundary vertices, we use the definition (3.12) of the discrete boundary Hodge operator to infer that

v∈V ∂ p v q v e∈Ev f (e) β•n f (e),c(v) d f = v∈V ∂ p v q v c(v) ∇•β dc -p ∂ , H ∂ β•n (q ∂ ) (VF) ∂ .
Proof of (3.13b). Using (3.10)- (3.11), the definition of [[•]] e , and the adjunction property between GRAD and DIV, we infer that

p, I E C β (GRAD(q)) VC = v∈V e∈Ev p v [[q]] e 1 2 (Λ ve -1)β e , q, DIV(I V F β (p)) VC = e∈E v∈Ve [[q]] e p v 1 2 (1 + Λ ve )β e .
Exchanging the summations in the first line leads to

p, I E C β (GRAD(q)) VC + q, DIV(I V F β (p)) VC = e∈E v∈Ve [[q]] e p v Λ ve β e .
Since v∈Ve p v Λ ve = [[p]] e Λ e owing to (Λ1), we infer (3.13b).

Analysis: coercivity, consistency, and error bound

We define the following stability norm for all q ∈ V:

|||q||| 2 a,V := τ -1 |||q||| 2 2,V + |||q||| 2 upw,β + |||q||| 2 |β•n| , (3.14)
where τ > 0 results from assumption (β1), |||•||| 2,V is defined in Section 2.4, and we can define |||q||| 2 upw,β := q, q upw,β from assumption (I2), and

|||q||| 2 |β•n| := q ∂ , H ∂ |β•n| (q ∂ ) (VF) ∂ from assumption (H ∂ ).
|||q||| 2 a,V ≤ A β (q, q), ∀q ∈ V, (3.15
)

with = 1 2 . Consequently, (3.4) is well-posed. Proof. Let q ∈ V. Since (I1)-(I2) imply that q, I E C β (GRAD(q)) VC = 1 2 q, H VC -∇•β (q) VC + 1 2 q, q upw,β + 1 2 q ∂ , H ∂ β•n (q ∂ ) (VF) ∂ ,
we infer that the quantity A β (q, q) can be rewritten as

A β (q, q) = 1 2 q, H VC -∇•β (q) VC + 1 2 q, q upw,β + 1 2 q ∂ , H ∂ β•n (q ∂ ) (VF) ∂ + q ∂ , H ∂ (β•n) -(q ∂ ) (VF) ∂ .
Owing to (H ∂ ), the last two terms on the right-hand side can be recombined to yield

A β (q, q) = 1 2 q, H VC -∇•β (q) VC + 1 2 q, q upw,β + 1 2 |||q||| 2 |β•n| . (3.16)
Since (I1) and (β1) imply q, H VC -∇•β (q) VC ≥ τ -1 |||q||| 2 2,V , (3.15) holds and (3.4) is well-posed.

We now turn to the consistency of the CDO scheme (3.4) using commutators in the spirit of Bossavit [START_REF] Bossavit | Computational electromagnetism and geometry[END_REF],

Hiptmair [START_REF] Hiptmair | Discrete Hodge operators[END_REF], and [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF]. To write the consistency error, we consider the reduction map R V : L 1 (Ω) → V such that ( R V (p)) v equals the mean-value of p in the dual cell c(v), and the following three commutators:

I V F β , R (q) := R F (βq) -I V F β ( R V (q)), (3.17a)
H VC -∇•β , R (q) := R C ((-∇•β)q) -H VC -∇•β ( R V (q)), (3.17b)
H ∂ (β•n) + , R (q) := R F ∂ ((β•n) + q) -H ∂ (β•n) + (( R V (q)) ∂ ), (3.17c)
for all q ∈ H s (Ω), s > 1 2 , so that q is in the domain of the maps R F and R F ∂ .

Lemma 3.3 (Error bound).

Let p ∈ V β (Ω) be the unique solution of (3.1) and let p be the unique solution of (3.4). Assume p ∈ H s (Ω), s > 1 2 . Then, under the assumptions of Lemma 3.2, the following holds:

|||p -R V (p)||| a,V ≤ sup q∈V;|||q||| a,V =1 E β (p, q), (3.18)
with consistency error defined as follows:

E β (p, q) := q, H VC -∇•β , R (p) VC -GRAD(q), I V F β , R (p) E F + q ∂ , H ∂ (β•n) + , R (p) (VF) ∂ . (3.19)
Proof. Owing to Lemma 3.2, it suffices to show that A β (p-R V (p), q) = E β (p, q). In the context of Friedrichs' systems, the derivation of the error bound hinges on integration by parts. In the CDO framework, we use the continuous and discrete Leibniz formulas, as well as the properties of the discrete differential operators.

We observe that

A β (p, q) = q, R C (β•∇p) VC + q ∂ , R F ∂ ((β•n) -p) (VF) ∂ = q, R C (∇•(βp)) VC + q, R C (-(∇•β)p) VC + q ∂ , R F ∂ ((β•n) -p) (VF) ∂ = q, DIV(R F (βp)) VC + q, R C (-(∇•β)p) VC + q ∂ , R F ∂ ((β•n) + p) (VF) ∂ = -GRAD(q), R F (βp) E F + q, R C (-(∇•β)p) VC + q ∂ , R F ∂ ((β•n) + p) (VF) ∂ ,
where we have used the continuous Leibniz formula (recall that p is in the graph space), the discrete commuting property (2.2b), the fact that Moreover, setting p = R V (p), we observe that

β•n = (β•n) + -(β•n) -,
A β ( p, q) = q, I E C β (GRAD( p)) VC + q ∂ , H ∂ (β•n) -( p ∂ ) (VF) ∂ = q, DIV(I V F β ( p)) VC + q, H VC -∇•β ( p) VC + q ∂ , H ∂ (β•n) + ( p ∂ ) (VF) ∂ = -GRAD(q), I V F β ( p) E F + q, H VC -∇•β ( p) VC + q ∂ , H ∂ (β•n) + ( p ∂ ) (VF) ∂ ,
where we have used the discrete Leibniz formula, assumption (H ∂ ) (linearity) together with Let p be the unique solution of (3.1) and let p be the unique solution of (3.4). Assume that p ∈ H 1 (Ω). Then, the following holds:

β•n = (β•n) + - (β•n) -,
|||p -R V (p)||| a,V (|β| 1 2 + h 1 2 τ 1 2 |∇•β| )h 1 2 |p| H 1 (Ω) , (3.20)
with stability constant defined in Lemma 3.2, |β| := ||β|| L ∞ (Ω) , and

|∇•β| := ||∇•β|| L ∞ (Ω) .
Proof. We need to bound the three terms in the right-hand side of (3.19) for all q ∈ V such that |||q||| a,V = 1. A direct calculation shows that

q, H VC -∇•β , R (p) VC = v∈V q v c(v) (-∇•β) (p -( R V (p)) v ) dc.
The definition of R V together with the Poincaré inequality and the multiplicative trace inequality imply that

||q -( R V (q)) v || L 2 (c(v)) + h 1 2 c(v) ||q -( R V (q)) v || L 2 (∂c(v)) h c(v) |q| H 1 (c(v)) , (3.21)
for all q ∈ H 1 (Ω) and all v ∈ V. Hence, we infer that

| q, H VC -∇•β , R (p) VC | (τ -1 2 |||q||| 2,V )(h 1 2 τ 1 2 |∇•β| )h 1 2 |p| H 1 (Ω) .
Turning to the second term in (3.19), a direct calculation using (Λ1), the fact that p is single-valued on f (e), and recalling the definition of β e shows that

-GRAD(q), I V F β , R (p) E F = e∈E [[q]] e v∈Ve 1 2 (1 + Λ β ve ) f (e) β•n f (e) (p -( R V (p)) v ) d f , ( 3.22)
whence we infer using (3.21) and the fact that

|Λ ve | ≤ 1 that | GRAD(q), I V F β , R (p) E F |   e∈E [[q]] 2 e f (e) |β•n f (e) | d f   1 2 |β| 1 2 h 1 2 |p| H 1 (Ω) .
Owing to the triangle inequality, Lemma 3.1, and (Λ2), we infer that Using the assumption L β τ -1 , mesh regularity, and the definition of the discrete norm

  e∈E [[q]] 2 e f (e) |β•n f (e) | d f   1 2 ≤ c -1 2 Λ |||q||| upw,β +   e∈E [[q]] 2 e f (e) |β•n f (e) | d f -|β e |   1
|||•||| 2,V leads to   e∈E [[q]] 2 e f (e) |β•n f (e) | d f   1 2 |||q||| upw,β + τ -1 2 |||q||| 2,V .
Finally, a direct calculation shows that

q ∂ , H ∂ (β•n) + , R (p) (VF) ∂ = v∈V ∂ q v f ∂ (v) (β•n) + (p -( R V (p)) v ) d f ∂ , so that | q ∂ , H ∂ (β•n) + , R (p) (VF) ∂ | |||q||| |β•n| |β| 1 2 h 1 2 |p| H 1 (Ω)
. This completes the proof.

Remark 3.6 (Localization). The error estimate (3.20) can be localized to dual mesh cells.

Advection-diffusion

This section addresses the derivation and analysis of vertex-based CDO schemes for the advection-diffusion problem (1.1). The diffusion tensor λ takes symmetric, uniformly positive definite values. For simplicity, we assume that λ is constant in each primal cell c ∈ C with minimal and maximal eigenvalues λ ,c and λ ,c , respectively, and local anisotropy ratio ρ c = λ ,c /λ ,c ≥ 1. The analysis can be extended to locally Lipschitz diffusion tensors.

Preliminaries: boundary penalty for pure diffusion

In this section, we consider the pure-diffusion version of the model problem (1.1) with β ≡ 0: Formally, the weak formulation is as follows: For all q ∈ H 1 (Ω),

-∇•(λ∇p) = s a.e. in Ω, (4.1a
Ω ∇q•λ•∇p dΩ - ∂Ω q(n•λ•∇p) d∂Ω + ∂Ω ηqp d∂Ω = Ω qs dΩ + ∂Ω ηqp D d∂Ω, ( 4.2)
with some boundary penalty parameter η. It is also possible to consider a symmetric bilinear form on the left-hand side. Symmetry is an important property when invoking duality arguments for pure diffusion problems; it is also a relevant property when inverting the linear system. It is less important in the presence of advection.

CDO scheme

The vertex-based CDO scheme with weakly enforced boundary conditions is formulated in terms of a discrete

Hodge operator

H E F λ : E → F,
A λ (p, q) = q, s VC + q ∂ , φ D (VF) ∂ , ∀q ∈ V, (4.3)
with bilinear form such that

A λ (p, q) := GRAD(q), H E F λ GRAD(p) E F -q ∂ , N ∂ λ GRAD(p) (VF) ∂ + η 0 q ∂ , H ∂ λ/h (p ∂ ) (VF) ∂ , ( 4.4)
where η 0 > 0 is a real number to be chosen large enough (see below), s = R C (s), and

φ D = η 0 H ∂ λ/h (R V ∂ (p D )
). The bilinear form (4.3) extends that of [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF] where the Dirichlet boundary condition was strongly enforced.

The discrete Hodge operator

H E F λ is assembled cell-wise from local operators H (EF)c λ : E c → F c for all c ∈ C, so that g 1 , H E F λ (g 2 ) E F = c∈C g 1 , H (EF)c λ (g 2 ) (EF)c (4.5)
for all g 1 , g 2 ∈ E. Similarly, the discrete normal flux operator N ∂ λ is assembled face-wise from local operators

N ∂ f λ : E c → F ∂ f for all f ∈ F ∂ , where c = c(f ) is the primal cell containing the primal boundary face f , so that q ∂ , N ∂ λ (g) (VF) ∂ = f ∈F ∂ q ∂ , N ∂ f λ (g) (VF) ∂ f , ( 4.6)
for all q ∂ ∈ V ∂ and all g ∈ E. Note that this implies that N ∂ λ (g), for all g ∈ E, only depends on the components of g attached to an edge of a cell having a boundary face. The discrete boundary penalty operator H ∂ λ/h is such that, still with c = c(f ),

H ∂ λ/h (q ∂ ) f ∂ (v) := q ∂ v f ∈F ∂ v | f ∂ f (v)|λ ,c h -1 c , ( 4.7)
for all q ∂ ∈ V ∂ and all v ∈ V ∂ , where

F ∂ v := {f ∈ F ∂ | v ∈ f }.
Note that H ∂ λ/h is algebraically represented by a diagonal matrix.

Example

Let us give a concrete example of CDO scheme. We consider a reconstruction operator

L Ec : E c → L ∞ (c) for all c ∈ C. The discrete Hodge operator in each cell c ∈ C is defined such that g 1 , H (EF)c λ (g 2 ) (EF)c = c L Ec (g 1 )•λ•L Ec (g 2 ) dc, ( 4.8)
for all g 1 , g 2 ∈ E c , while the discrete normal flux operator in each boundary face f ∈ F ∂ is defined as follows (with c = c(f ) the primal cell containing f ):

N ∂ f λ (g) f ∂ f (v) := f ∂ f (v) n•λ•L Ec (g) d f , ( 4.9)
for all v ∈ V ∂ f and all g ∈ E c . The reconstruction operator has to satisfy some properties stated in Lemma 4.1 below. One possibility is to consider the reconstruction proposed by Codecasa et al. [START_REF] Codecasa | A new set of basis functions for the discrete geometric approach[END_REF], see also [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF][START_REF] Bonelle | Low-order reconstruction operators on polyhedral meshes: Application to Compatible Discrete Operator schemes[END_REF], whereby L Ec (g) is piecewise constant on each diamond d(e) ∩ c, e ∈ E c (see Figure 2, left panel).

Design conditions

More generally, the design conditions on H

∈ E c , c H λ ,c |||g||| 2 2,Ec ≤ g, H (EF)c λ (g) (EF)c ≤ c -1 H λ ,c |||g||| 2 2,Ec . (4.10) (H2) [P 0 -consistency] R Fc (λ•G) = H (EF)c λ (R Ec (G)) for any constant field G in c.
The design conditions on N ∂ f λ are as follows: For all f ∈ F ∂ , with c = c(f ), (N1) [Boundedness] There exists c N , uniform with respect to the mesh and the model parameters, such that, for all g ∈ E c ,

v∈V ∂ f | f ∂ f (v)| -1 N ∂ f λ (g) f ∂ f (v) 2 ≤ c N λ ,c h -1 ,c g, H (EF)c λ (g) (EF)c . (4.11) (N2) [P 0 -consistency] R F ∂ f (n•λ•G) = N ∂ f λ (R Ec (G))
for any constant field G in c. Lemma 4.1 (Design conditions). Let the discrete Hodge and normal flux operators be defined by (4.8) and (4.9), respectively. Assume that the reconstruction operator is such that: Proof. For the proof of (H1)-(H2), see [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF][START_REF] Bonelle | Low-order reconstruction operators on polyhedral meshes: Application to Compatible Discrete Operator schemes[END_REF][START_REF] Codecasa | A new set of basis functions for the discrete geometric approach[END_REF]. To prove (N1), fix f ∈ F ∂ and observe that

v∈V ∂ f | f ∂ f (v)| -1 N ∂ f λ (g) f ∂ f (v) 2 = v∈V ∂ f | f ∂ f (v)| -1 f ∂ f (v) n•λ•L Ec (g) d f 2 ≤ v∈V ∂ f λ ,c ||λ 1/2 L Ec (g)|| 2 L 2 ( f ∂ f (v)) ≤ c tr v∈V ∂ f λ ,c h -1 ,c ||λ 1/2 L Ec (g)|| 2 L 2 (c) = c tr #(V ∂ f )λ ,c h -1 ,c g, H (EF)c λ (g) (EF)c ,
where we have used the Cauchy-Schwarz inequality followed by a discrete trace inequality with c = c(f ) (since λ is constant and L Ec (g) is a piecewise-polynomial) and the definition of the discrete Hodge operator.

This proves (N1) with c N = c tr #(V ∂ f ) (observing that the cardinal number #(V ∂ f ) is uniformly bounded owing to mesh regularity). Finally, letting G be a constant field in c, (N2) follows from

N ∂ f λ (R Ec (G)) f ∂ f (v) = f ∂ f (v) n•λ•L E (R Ec (G)) d f = f ∂ f (v) n•λ•G d f = (R F ∂ f (n•λ•G)) f ∂ f (v) ,
for all v ∈ V ∂ f owing to property (ii) of the reconstruction operator.

Analysis

This section collects the main results concerning the analysis of the CDO scheme with boundary penalty.

To facilitate the reading, the proofs are postponed to Section 7. We define the following norms on E and V ∂ , respectively:

|||g||| 2 λ := g, H E F λ (g) E F , |||q ∂ ||| 2 λ/h = q ∂ , H ∂ λ/h (q ∂ ) (VF) ∂ . (4.12)
Observe that these norms can be localized as

|||g||| 2 λ = c∈C |||g c ||| 2 λ,c with |||g c ||| 2 λ,c = g c , H (EF)c λ (g c ) (EF)c for all g c ∈ E c , and as |||q ∂ ||| 2 λ/h = f ∈F ∂ |||q ∂ f ||| 2 λ/h,f with |||q ∂ f ||| 2 λ/h,f = λ ,c h -1 c v∈V ∂ f | f ∂ f (v)|(q ∂ f,v ) 2 for all q ∂ f ∈ V ∂ f .
The stability of the CDO scheme (4.3) is expressed in the following norm:

|||q||| 2 d,V := |||GRAD(q)||| 2 λ + |||q ∂ ||| 2 λ/h , ∀q ∈ V. (4.13)

Lemma 4.2 (Coercivity and well-posedness). Assume (H1) and (N1). Then, provided η

0 ≥ 1 + 1 2 c 2 N , the following holds: |||q||| 2 d,V ≤ A λ (q, q), ∀q ∈ V, (4.14)
with = 1 2 . Consequently, (4.3) is well-posed.

We now address the consistency of the scheme (4.3). We assume that the exact solution p is in H s (Ω), s > 3 2 (the regularity assumption can be localized to mesh cells), and we consider the (classical) de Rham map R V . We define the following two commutators:

H E F λ , R (∇p) := R F (λ•∇p) -H E F λ (GRAD(R V (p))), (4.15a) N ∂ λ , R (∇p) := R F ∂ (n•λ•∇p) -N ∂ λ (GRAD(R V (p))). (4.15b)
Lemma 4.3 (Consistency). Let p be the unique solution of (4.1) and let p be the unique solution of (4.3).

Assume p ∈ H s (Ω), s > 3 2 . Then, under the assumptions of Lemma 4.2, the following holds:

|||p -R V (p)||| d,V ≤ sup q∈V;|||q||| d,V =1 E λ (p, q), (4.16)
with consistency error

E λ (p, q) := GRAD(q), H E F λ , R (∇p) E F -q ∂ , N ∂ λ , R (∇p) (VF) ∂ . (4.17)
Theorem 4.4 (Convergence rate). Let p be the unique solution of (4.1) and let p be the unique solution of

(4.

3). Assume (H1)-(H2) and (N1)-(N2)

. Assume p ∈ H 2 (Ω). Then, the following holds:

|||p -R V (p)||| d,V   c∈C ρ c λ ,c h 2 c |p| 2 H 2 (c)   1 2 . (4.18)

CDO Scheme for advection-diffusion

Vertex-based CDO schemes for the advection-diffusion problem (1.1) hinge on the discrete bilinear form A β,λ := A β + A λ with A β defined by (3.5) and A λ by (4.4). The discrete problem consists in finding p ∈ V

such that A β,λ (p, q) = q, s VC + q ∂ , φ D (VF) ∂ , ∀q ∈ V, (4.19
)

with s = R C (s) and φ D = R F ∂ ((β•n) -p D ) + η 0 H ∂ λ/h (R V ∂ (p D ))
. The Tonti diagram of the vertex-based CDO scheme (4.19) is presented in Figure 4. Variants, such as using the conservative form of the advective derivative or including a reactive term, can be considered as well; see Remarks 3.2 and 3.3. We define the stability norm on V as |||q||| 2 ad,V := |||q||| 2 a,V + |||q||| 2 d,V with advection-related stability norm defined by (3.14) and diffusion-related stability norm defined by (4.13). Assume p ∈ H s (Ω), s > 3 2 . Then, under the assumptions of Lemma 4.5, the following holds:

p oe V g oe E oe  F s oe  C p ˆoe V  F ˆGRAD Á DIV H E F ⁄ ' H ⁄/h /H (-•n) ≠ I E  C - N ⁄ RHS p oe V g oe E p ˆoe V ˆ GRAD H E F ⁄ ' H ⁄/h /H (-•n) ≠ I V  F - N ⁄

Lemma 4.5 (Coercivity and well-posedness). Assume (β1), (I1)-(I2), and (H ∂ ) for the advection-related terms together with (H1) and (N1) for the diffusion-related terms. Then, provided η

0 ≥ 1 + 1 2 c 2 N , the following holds: |||q||| 2 ad,V ≤ A β,λ (q, q), ∀q ∈ V, ( 4
|||p -R V (p)||| ad,V ≤ sup q∈V;|||q||| ad,V =1 E β,λ (p, q), (4.21
)

with consistency error E β,λ (p, q) = E β (p, q) + E λ (p, q), with E β (p, q) defined by (3.19) with R V in lieu of R V ,
and E λ (p, q) defined by (4.17).

Proof. Combine Lemma 3.3 with Lemma 4.3 (note that H 1 (Ω) ⊂ V β (Ω)).

Example: CDO scheme with Péclet-based upwinding

For all e ∈ E, we define the (algebraic) edge Péclet number as Pe e = λ To write the error estimate, we introduce one last geometric object d(e), for all e ∈ E, which is the so-called diamond around e formed by the two pyramids of apex v ∈ V e and (non-planar) basis f (e), see and let p be the unique solution of (4.19). Assume p ∈ H 2 (Ω). Then, the following holds: which is now applied in the diamonds around edges, is then replaced by

|||p -R V (p)||| ad,V   c∈C ρ c λ ,c h 2 c |p| 2 H 2 (c)   1 2 +   e∈E (τ |∇•β| 2 L ∞ (d(
||q -(R V (q)) v || L 2 (d(e)) + h 1 2 e ||q -(R V (q)) v || L 2 ( f (e)) h e |q| H 1+ (d(e)) ,
for all q ∈ H 2 (Ω), all e ∈ E, and all v ∈ V e . The second adaptation is related to the change in the |||•||| upw,β semi-norm owing to the use of Péclet-based upwinding. We bound again the three terms in the right-hand side of (3.19) for all q ∈ V such that |||q||| ad,V = 1. For the first term, we readily infer that

| q, H VC -∇•β , R (p) VC | (τ -1 2 |||q||| 2,V )   e∈E τ |∇•β| 2 L ∞ (d(e)) h 2 e |p| 2
e∈E >1 (•)   e∈E >1 [[q]] 2 e f (e) |β•n f (e) | d f   1 2   e∈E >1 |β•n| ,e h e |p| 2 H 1+ (d(e))   1 2
.

For all e ∈ E >1 , property (Θ2) implies that Λ e β e ≥ α|β e |. Then, still proceeding as in the proof of Theorem 3.4, we infer that

e∈E >1 (•) |||q||| upw,β + τ -1 2 |||q||| 2,V   e∈E >1 |β•n| ,e h e |p| 2 H 1+ (d(e))   1 2
.

Furthermore, we observe that

e∈E ≤1 (•)   e∈E ≤1 [[q]] 2 e h e λ e   1 2   e∈E ≤1 h -1 e λ -1 e |β•n| ,e |β e |h e |p| 2 H 1+ (d(e))   1 2
.

Owing to mesh regularity, the definition of λ e , and (H1), we infer that the first factor in the right-hand side is bounded by |||GRAD(q)||| λ , while the second factor is bounded by

e∈E |β•n| ,e h e Pe e |p| 2 H 1+ (d(e))
| GRAD(q), I V F β , R (p) E F | |||q||| ad,V   e∈E |β•n| ,e h e min(1, Pe e )|p| 2 H 1+ (d(e))   1 2
.

Finally, the boundary term is bounded as before. Remark 4.2 (Boundary term). It is also possible to modify the discrete boundary Hodge operator so as to enforce the boundary condition using a Péclet-based upwinding; details are omitted for brevity.

Divergence-free advection

In this section, we extend the analysis to the case of a divergence-free velocity field β under assumption (β2);

recall that this assumption provides a real number τ > 0 and a function

ζ ∈ W 1,∞ (Ω) such that ζ ≥ 1 a.e.
in Ω. The advection-related stability norm |||•||| a,V is still defined by (3.14) (where now τ results from (β2)).

The only relevant change in the analysis is that stability (and well-posedness) is now achieved by means of an inf-sup condition instead of a coercivity argument. Since consistency and boundedness hold in the same form as before, inf-sup stability suffices to establish the error upper bounds, so that the convergence rates

Pure advection

Along with (I1)-(I2), we introduce a third property for the discrete contraction operators:

(I3) [Multiplication by ζ]

There are c 1 , c 2 , c 3 , uniform with respect to the mesh-size and the functions β and ζ, such that the following holds for all q ∈ V:

|||ζq||| 2 upw,β + |||ζq||| 2 |β•n| ≤ c 1 ζ 2 (|||q||| 2 upw,β + |||q||| 2 |β•n| ) + c 2 ω β τ -1 |||q||| 2 2,V , (5.1a)
A β (q, ζq) ≥ 1 4 |||q||| 2 a,V -c 3 ω β A β (q, q), (5.1b
)

with ζ := ||ζ|| L ∞ (Ω) and ζq ∈ V such that (ζq) v := ζ(v)q v for all v ∈ V.

Lemma 5.1 (Inf-sup stability). Under hypotheses (β2), (I1)-(I2)-(I3), and (H ∂ ), the following holds:

|||q||| a,V ≤ sup r∈V A β (q, r) |||r||| a,V , ∀q ∈ V, (5.2
)

with = 1 4 max(ζ 2 + c 2 ω β , c 1 ζ 2 ) 1 2 + c 3 ω β -1
.

Proof. We take r = ζq + c 3 ω β q ∈ V. Owing to the triangle inequality, (5.1a), and the obvious bound

|||ζq||| 2,V ≤ ζ |||q||| 2,V , we infer that |||r||| a,V ≤ |||ζq||| a,V + c 3 ω β |||q||| a,V ≤ max(ζ 2 + c 2 ω β , c 1 ζ 2 ) 1 2 + c 3 ω β |||q||| a,V .
Moreover, owing to (5.1b), we infer that

A β (q, r) = A β (q, ζq) + c 3 ω β A β (q, q) ≥ 1 4 |||q||| 2 a,V ,
whence we infer (5.2). We now verify property (I3) in the context of the CDO scheme with upwinding studied in Section 3.2.

Lemma 5.2 ((I3) with upwinding). Assume (β2) and (Λ1). Let the discrete contraction and surface Hodge

operators be given by (3.10)-(3.11)-(3.12). Then, (I3) holds.

Proof. To prove property (5.1a), we observe that, for all q ∈ V, since Λ e β e ≥ 0 for all e ∈ E,

|||ζq||| 2 upw,β = e∈E [[ζq]] 2 e Λ e β e ≤ 2 e∈E { {ζ} } 2 e [[q]] 2 e + { {q} } 2 e [[ζ]] 2 e Λ e β e ,
where , this completes the proof of (5.1a). Proof of (5.1b). The idea of the proof consists of writing A β (q, ζq) in the form A ζβ (q, q) plus a perturbation which can be bounded by the variations of ζ. A straightforward computation proceeding as in the proof of Lemma 3.1 shows that A β (q, ζq) = T 1 + T 2 + T 3 with

{ {ζ} } e = 1 2 v∈Ve ζ(v), [[ζ]] e = v∈Ve ι f (e),c(v) ζ(v), { {q} } e = 1 2 v∈Ve q v ,
T 1 = v∈V e∈Ev 1 2 q 2 v f (e) (-ζβ)•n f (e) d f + e∈E 1 2 [[q]] 2 e ζ e Λ e β e + v∈V ∂ q 2 v f ∂ (v) ζ(β•n) -d f , T 2 = v∈V e∈Ev 1 2 q v g e f (e) (ζ(v) -ζ)β•n f (e) d f + v∈V e∈Ev 1 2 q v g e (ζ e -ζ(v))Λ ve β e , T 3 = v∈V ∂ q 2 v f ∂ (v) (ζ(v) -ζ)(β•n) -d f ,
with g = GRAD(q) and ζ e the mean-value of ζ in e. Since ζ ≥ 1, still proceeding as in the proof of Lemma 3.1 and using now (β2) leads to the bound T 1 ≥ 1 2 |||q||| 2 a,V . Furthermore, using Cauchy-Schwarz inequalities, (Λ2), and mesh regularity, we obtain |T 2 | ω β |||q||| upw,β (τ -1 2 |||q||| 2,V ). Proceeding similarly leads Lemma 5.4 (Multiplication by ζ). Assume that (H1) and (N1) hold. The following holds for all q ∈ V with g = GRAD(q):

|||ζq||| 2 d,V ≤ 2ζ 2 |||q||| 2 d,V + 2c 2 5 ω λ τ -1 |||q||| 2 2,V , (5.5a) A λ (q, ζq) ≥ |||g||| 2 λ + η 0 |||q ∂ ||| 2 λ/h -c N ζ |||g||| λ |||q ∂ ||| λ/h -c 5 ω 1 2 λ |||g||| λ (τ -1 2 |||q||| 2,V ).
(5.5b)

Proof. Proof of (5.5a). The definition of

H ∂ λ/h implies that |||ζq ∂ ||| 2 λ/h = ζq ∂ , H ∂ λ/h (ζq ∂ ) (VF) ∂ ≤ ζ 2 |||q ∂ ||| 2 λ/h
. Furthermore, owing to the cell-wise assembly of the discrete Hodge operator H E F λ and using the triangle inequality, we infer that

|||GRAD(ζq)||| 2 λ = c∈C GRAD c (ζq), H (EF)c λ •GRAD c (ζq) (EF)c = c∈C 2ζ 2 c GRAD c (q), H (EF)c λ •GRAD c (q) (EF)c + c∈C 2δ c (q) 2 ≤ 2ζ 2 |||GRAD(q)||| 2 λ + c∈C 2δ c (q) 2 ,
δ c (q) 2 ≤ c -1 H λ ,c h c e∈Ec   v∈Ve ι v,e (ζ(v) -ζ c )q v   2 ≤ 2c -1 H λ ,c L 2 ζ h 3 c e∈Ec v∈Ve q 2 v ≤ 2c -1 H N V,E λ ,c L 2 ζ |||q||| 2 2,Vc .
Combining the above bounds leads to (5.5a).

Proof of (5.5b). Using (N1), (7.1), and ζ ≥ 1, we infer that

A λ (q, ζq) = GRAD(ζq), H E F λ •GRAD(q) E F -(ζq) ∂ , N ∂ λ •GRAD(q) (VF) ∂ + η 0 (ζq) ∂ , H ∂ λ/h (q ∂ ) (VF) ∂ ≥ GRAD(ζq), H E F λ •GRAD(q) E F -c N |||GRAD(q)||| λ |||(ζq) ∂ ||| λ/h + η 0 |||q ∂ ||| 2 λ/h ≥ GRAD(ζq), H E F λ •GRAD(q) E F -c N ζ |||GRAD(q)||| λ |||q ∂ ||| λ/h + η 0 |||q ∂ ||| 2 λ/h .
Moreover, owing to the cell-wise assembly of the discrete Hodge operator H E F λ and proceeding as above, we infer that

GRAD(ζq), H E F λ •GRAD(q) E F = c∈C GRAD c (ζq), H (EF)c λ •GRAD c (q) (EF)c ≥ |||GRAD(q)||| 2 λ + c∈C GRAD c (w), H (EF)c λ •GRAD c (q) (EF)c
Since H

(EF)c λ is self-adjoint and monotone, we infer that

| GRAD c (w), H (EF)c λ •GRAD c (q) (EF)c | ≤ δ c (q)|||GRAD c (q)||| λ,c , so that GRAD(ζq), H E F λ •GRAD(q) E F ≥ |||GRAD(q)||| 2 λ -   c∈C δ c (q) 2   1 2
|||GRAD(q)||| λ .

Using the above bound on δ c (q) yields (5.5b).

Numerical results

In this section, we investigate numerically CDO advection-diffusion schemes on four families of successivelyrefined, polyhedral meshes of the unit cube Ω = (0, 1) 3 ; see Figure 5 for an example of mesh within each family. These mesh families have been proposed in the FVCA benchmark [START_REF] Eymard | 3D benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF], see also [START_REF] Bonelle | Compatible Discrete Operator schemes on polyhedral meshes for elliptic and Stokes equations[END_REF].

Hexahedral (H)

Linear Pressure and Linear Velocity

This first test case enables us to check the correct implementation of the CDO scheme. The analytic solution of the problem is p

= x + y + z ≠ 3 2 and u = S W W W U x + y + z ≠(x + y + z) 0 T X X X V .
Two sets of BCs are tested: (UnWt) and (UtPr). It is worth mentioning that these BCs are nonhomogeneous.

Expected Results. From the error analysis performed in [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF], we expect that the pressure has a secondorder convergence rate, and that the pressure gradient, the velocity and the vorticity have a first-order convergence rate. Thus, we should capture exactly the pressure field, its gradient and also the vorticity. By exactly, we mean up to machine precision (MP). The velocity should be approximated with a first-order convergence rate.

Numerical Results.

The numerical results are in accordance with the theoretical results derived in [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF]. For the two sets of BCs, the vorticity, the pressure and its gradient are exactly captured by the numerical scheme. The error on the velocity field decreases at least with a first-order convergence rate. Except for the tetrahedral and checkerboard sequences, we observe a superconvergence. In the specific case of (UtPr) BCs and uniform meshes, the velocity is exactly captured. Table 9 summarizes these numerical results. The most accurate results are obtained with the Cartesian sequence and the less accurate with the Kershaw sequence (see Figure 14).

H TU T PrT PrG CB K Er Hfl (u) (UnWt) 1.9 2.6 1.1 1.9 1.6 1.0 1.9 The error with respect to the exact solution p is measured using the following two quantities:

Err 2,V := |||p -R V (p)||| 2,V |||R V (p)||| 2,V , Err ad,V := |||p -R V (p)||| ad,V |||R V (p)||| ad,V .
In our numerical tests, the integrals for the source term and the boundary data are computed using a fourth-order quadrature on elementary sub-simplices of each polyhedral cell.

Anisotropic diffusion and variable advection velocity

We consider the conservative form of the scheme (4. [START_REF] Di Pietro | A discontinuous-skeletal method for advection-diffusionreaction on general meshes[END_REF], where the bilinear form A β is given by (3.8). The exact solution is p(x, y, z) = 1 + sin(πx) sin π y + 1 2 sin π z + 1

3

, and the diffusive tensor λ and the velocity field β are equal to (in the canonical basis of R 3 )

λ =     1 0.5 0 0.5 1 0.5 0 0.5 1     , β =     y -1/2 1/2 -x z     ,
so that the velocity field satisfies hypothesis (β1) for the conservative form (see Remark 3.2). We consider the discrete contraction operator I V F β built using full upwinding as in Section 3.2 and Péclet-based upwinding as in Section 4.3 using the Sharfetter-Gummel map. Figure 6 presents the numerical results, which reflect the theoretical analysis with convergence rates between one and two. The use of Péclet-based upwinding leads to lower errors than full upwinding; the improvement is more pronounced on the SkP mesh sequence than on the other sequences, and is observed on the finer meshes.

ad,V 1
Figure 6: Test case 6.1: Convergence curves for the two error measures on the four mesh families using full upwinding (dashed lines) or Péclet-based upwinding (solid lines); first-and second-order slopes are indicated.

Exponential boundary layer with constant advection velocity

The second test case investigates the behavior of the CDO scheme in the presence of an exponential boundary layer resulting from small diffusion. We consider an isotropic diffusive tensor, i.e. λ = λId, and a constant vector field β with components (2, 3, 0), so that assumption (β2) is satisfied. The exact solution is p(x, y, z) = (xe 2(x-1) λ

)(y 2e 3(y-1) λ

) and exhibits a boundary layer near x = 1 and y = 1 when λ tends to 0.

Figure 7 reports numerical results for λ = 1 (solid lines) and λ = 10 -4 (dashed lines). Note that in this second case, the considered meshes do not resolve the boundary layer. The transition between the two convergence regimes as predicted by Theorem 4.7 is clearly visible. The present test case is also considered by Da Veiga et al. [START_REF] Beirão Da Veiga | A unified approach for handling convection terms in finite volumes and mimetic discretization methods for elliptic problems[END_REF] on the same SkP mesh sequence with a different scheme, where similar convergence rates are reported but with somewhat larger error values.

7 Analysis of CDO schemes for pure diffusion 7.1 Proof of Lemma 4.2

Property (N1) implies that, for all (q ∂ , g) ∈ V ∂ × E, line. As a result, we infer that A λ (q, q) ≥ |||GRAD(q)||| 2 λc N |||GRAD(q)||| λ |||q||| λ/h + η 0 |||q||| 2 λ/h . (7.2)

q ∂ , N ∂ λ (g) (VF) ∂ = f ∈F ∂ q ∂ f , N ∂ f λ (g c ) (VF) ∂ f = f ∈F ∂ v∈V ∂ f q ∂ f,v N ∂ f λ (g c ) f ∂ (v) ≤ c N f ∈F ∂ |||q ∂ f ||| λ/h,
To conclude, we use the quadratic inequality x 2 -2γxy + δy 2 ≥ δ-γ 2 1+δ (x 2 + y 2 ) (valid for any real numbers x, y, γ, δ with δ ≥ 0) with γ = 1 2 c N and δ = η 0 and observe that the choice η 0 ≥ 1 + 1 2 c 2 N implies δ ≥ 1 + 2γ 2 so that δ-γ 2 1+δ ≥ 1 2 . Finally, the well-posedness of (4.3) follows from (4.14).

Proof of Lemma 4.3

Owing to Lemma 4.2, it suffices to show that A λ (p -R V (p), q) = E λ (p, q). To prove this, we observe that q, s VC = -q, R C (∇•(λ∇p)) VC = -q, DIV(R F (λ∇p))

VC -q ∂ , R F ∂ (n•λ•∇p) (VF) ∂ ,
owing to (2.2b), and we use (4.3) and (4.4) to conclude, as well as R V ∂ (p D ) = (R V (p)) ∂ .

Proof of Theorem 4.4

Let T 1 , T 2 be the two terms in the right-hand side of (4.17). Recall that GRAD(R V (p)) = R E (∇p). The term T 1 has already been bounded in [START_REF] Bonelle | Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes[END_REF]; we present here a somewhat simpler proof avoiding the algebraic identity on the inverse of the discrete Hodge operator. Let G c denote the mean-value of ∇p in c. Owing to the local assembly (4.5) and to (H2), we infer that .

T 1 = c∈C GRAD c (q), R Fc (λ•∇p) -H (
Turning to T 2 , we use the local assembly (4.6) and (N2) to infer that, with c = c(f ),

T 2 = f ∈F ∂ q ∂ f , R F ∂ f (n•λ•∇p) -N ∂ f λ (R Ec (∇p)) (VF) ∂ f = f ∈F ∂ q ∂ f , R F ∂ f (n•λ•(∇p -G c )) (VF) ∂ f - f ∈F ∂ q ∂ f , N ∂ f λ (R Ec (∇p -G c )) (VF) ∂ f ,
and we denote by T 2,1 and T 2,2 the two terms in the right-hand side. The Cauchy-Schwarz inequality implies that

|T 2,1 | ≤ |||q ∂ ||| λ/h    f ∈F ∂ v∈V ∂ f λ ,c h -1 c ||∇p -G c || 2 L 1 ( f ∂ f (v))    1 2 ∂ ||| λ/h   f ∈F ∂ λ ,c h 2 c |p| 2 H 2 (c)   1 2
, while using (N1) and proceeding as above, we infer a similar bound on T 2,2 . The proof is complete since ρ c ≥ 1 by definition.

Acknowledgement. The authors thank Jérôme Bonelle (EDF R&D) for fruitful discussions.

1 Figure 1 :

 11 Figure 1: Illustration of primal and dual mesh entities

  is convenient to localize discrete objects to a primal cell or to a boundary face. Let c ∈ C. We define the local subsets V c := {v ∈ V | v ∈ ∂c} (collecting the vertices of the cell c) and E c := {e ∈ E | e ⊂ ∂c} (collecting the edges of the cell c). For all e ∈ E c , we define fc (e) := f (e) ∩ c as the portion of the dual face f (e) inside c (see Figure 2, left panel), and we set F c := { fc (e); e ∈ E c }. The vector space E c is composed of the DoFs of discrete circulation fields g ∈ E attached to E c ; similarly for V c and for F c . The de Rham maps R Ec and R Fc are such that (R Ec (g)) e = e (t e •g) de for all e ∈ E c , and (R Fc (φ)

  as the portion of the dual face f ∂ (v) inside f (see Figure2, right panel), and we set F ∂

Figure 2 :

 2 Figure 2: Illustration of local mesh entities for a cell (left) and a boundary face (right)

. 5 )

 5 and where we have set s := R C (s) and φ D := R F ∂ (φ D ). A synthetic presentation of the scheme (3.4) is the so-called Tonti diagram shown in Figure 3.

Figure 3 :

 3 Figure 3: Tonti diagram of the vertex-based CDO scheme for pure advection

Remark 3 . 1 (

 31 Discrete contraction operator). Recalling from Remark 2.1 that the discrete dual divergence operator DIV does not involve faces on the boundary ∂Ω, property (I1) is the discrete counterpart of the Leibniz formula -Ω p(∇•β)q dΩ = Ω p(β•∇q) dΩ -Ω p∇•(βq) dΩ, where the two rightmost terms in (3.6) form together the discrete counterpart of Ω p∇•(βq) dΩ. Furthermore, property (I2) is the discrete

1 2

 1 and V e := {v ∈ V | v ∈ e} for all e ∈ E, and we use the notation n f (e),c(v) = ι f (e),c(v) n f (e) for the unit normal to f (e) pointing outward c(v). For all e ∈ E and all v ∈ V e , we fix a real number Λ ve ∈ [-1, 1] (the algebraic upwinding parameter) such that the following holds: For all e ∈ E, (Λ1) v∈Ve Λ ve = 0, and setting Λ e := v∈Ve ι f (e),c(v) Λ ve , β e Λ e ≥ 0 holds. (Λ2) There exists c Λ > 0, uniform with respect to the mesh and the model parameters, such that β e Λ e ≥ c Λ |β e |.

Lemma 3 . 1 (

 31 Stability, (I1)-(I2)). Let the discrete contraction and surface Hodge operators be given by (3.10)-(3.11)-(3.

Lemma 3 . 2 (

 32 Coercivity and well-posedness). Under hypotheses (β1), (I1)-(I2), and (H ∂ ), the following holds:

  and the discrete adjunction property (2.1).

Theorem 3 . 4 (

 34 and the discrete adjunction property (2.1). The conclusion is straightforward. Convergence rate). Assume (β1). Let L β be the Lipschitz constant of β and assume that L β τ -1 . Let the discrete contraction and surface Hodge operators be given by (3.10)-(3.11)-(3.12). Assume (Λ1)-(Λ2).

2 .

 2 Let c ∈ C e ; the local dual face fc (e) consists of two triangles, say ff,c (e), each touching one of the two facesf of c sharing e. Set δ f,c (e) := ff,c (e) |β•n ff,c (e) | d fff,c (e) β•n ff,c (e) d f , so that 0 ≤ δ f,c (e) = 2 ff,c (e) (β•n ff,c (e) ) -d f . If (β•n ff,c(e) ) -takes positive values on ff,c (e), then 0 ≤ δ f,c (e) ≤ 2|β e |; otherwise, β•n f f,c (e) vanishes at some point in ff,c (e). Then, using the fact that β•n ff,c (e) is Lipschitz in ff,c (e) together with mesh regularity leads to 0 ≤ δ f,c (e) L β h e | ff,c (e)|. Since f (e) (•) ≥ f (e) (•), summing these bounds over c ∈ C e and the faces f leads to 0 ≤ f (e) |β•n f (e) | d ff (e) β•n f (e) d f ≤ f,c δ f,c (e) |β e | + L β h e | f (e)|.

  ) p = p D a.e. on ∂Ω. (4.1b)

  which is the discrete counterpart of the map g → λ•g, and the discrete boundary operators N ∂ λ : E → F ∂ (normal flux) and H ∂ λ/h : V ∂ → F ∂ (boundary penalty), which weakly enforce boundary conditions à la Nitsche and which are the discrete counterparts of the maps g → n•λ•g and p → (λ/h)p at the boundary, respectively. The discrete problem consists in finding p ∈ V such that

(

  EF)c λ are as follows: For all c ∈ C, (H1) [Stability] H (EF)c λ is self-adjoint and monotone, and there exists c H > 0, uniform with respect to the mesh and the model parameters, such that, for all g

L |||g||| 2 2 ,

 2 (i) [Stability] L Ec (g) is a piecewise-polynomial function in c, and there exists c L > 0, uniform with respect to the mesh-size, such that c L |||g||| 2 2,Ec ≤ ||L Ec (g)|| 2 L 2 (c) ≤ c -1 Ec for all g ∈ E c . (ii) [Partition of unity] L Ec (R Ec (G)) = G for any constant field G in c. (iii) [Dual consistency] c L Ec (g) dc = e∈Ec ( fc(e) n fc(e) d f )g e for all g ∈ E c . Then, (H1)-(H2) and (N1)-(N2) hold.

Figure 4 :

 4 Figure 4: Tonti diagram of the vertex-based CDO scheme for advection-diffusion weakly enforced boundary conditions

- 1 e

 1 |f (e)| -1 β e h e with λ e = max c∈Ce λ ,c , C e = {c ∈ C | e ⊂ ∂c}, and β e defined in Section 3.2. We then use (3.10)-(3.11) to define the discrete contraction operators I E C β and I V F β with Péclet-dependent upwinding parameter Λ ve = Θ(ι f (e),c(v) Pe e ), where the function Θ : R → R is such that (Θ1) Θ(x) + Θ(-x) = 0 and Θ(x) ≥ 0 for all x ∈ R ≥0 . (Θ2) There exists α > 0 such that Θ(x) ≥ α for all x ≥ 1 (the lower bound on x is arbitrary; changing its value only changes the constants in the error bounds). Note that (Θ1) implies (Λ1) since β e Λ e = 1 2 λe| f (e)| he v∈Ve ι f (e),c(v) Pe e Θ(ι f (e),c(v) Pe e ) ≥ 0. Since (Λ1) holds, Lemma 3.1 implies that (I1)-(I2) hold; hence, stability and well-posedness hold owing to Lemma 4.5. An example for the function Θ is the Sharfetter-Gummel map Θ(x) = coth x 2 -2 x , see Roos et al. [35] for further insight and examples. The function Θ is related to the function |A| introduced in [19] in the context of high-order face-based discretizations by the relation |A|(x) = xΘ(x).

Figure 2 (

 2 Figure 2 (left panel). Note that ∪ e∈E d(e) = Ω.

Theorem 4 . 7 (

 47 Convergence rate). Assume (β1). Let L β be the Lipschitz constant of β and assume that L β τ -1 . Let the discrete contraction and surface Hodge operators be given by (3.10)-(3.11)-(3.12). For the diffusion-related operators, assume (H1)-(H2) and (N1)-(N2). Let p be the unique solution of (1.1)

2 .

 2 Consider now the second term. Let E >1 := {e ∈ E | |Pe e | > 1} and E ≤1 := {e ∈ E | |Pe e | ≤ 1}. The summation in the right-hand side of (3.22) is split as e∈E >1 (•) + e∈E ≤1 (•). Proceeding as in the proof of Theorem 3.4, we infer that

  e | f (e)| -1 |β e |h e = |Pe e | ≤ 1 and | f (e)| h 2 e . Collecting the bounds on e∈E >1 and e∈E ≤1 leads to

Remark 4 . 1 (

 41 Limit regimes). In the advection-dominant regime with |Pe e | ≥ 1 for all e ∈ E, the error bound (4.22) behaves as h 1/2 (see Theorem 3.4), while, in the diffusion-dominant regime with |Pe e | ≤ h e for all e ∈ E, it behaves as h (see Theorem 4.4). The case where h e ≤ Pe e ≤ 1 corresponds to transition regimes and intermediate orders of convergence.

derived in Theorem 3 . 4

 34 for pure advection and in Theorem 4.7 for advection-diffusion still hold. In what follows, we consider the non-dimensional numbers ω β = L 2 ζ |β| hτ and ω λ = L 2 ζ λ τ , with L ζ the Lipschitz constant of ζ, |β| := ||β|| L ∞ (Ω) , and λ := max c∈C λ ,c .

Remark 5 . 1 (

 51 Factor ω β ). An upper bound on ω β yields a lower bound on . A simple upper bound is to replace h by a global length scale associated with Ω (i.e., h can be replaced by a global length scale in (5.1a) and (5.1b)). A sharper bound is ω β ≤ L ζ |β| τ under the mild assumption L ζ h ≤ 1 (meaning that h resolves the scale of spatial variations of ζ).

  where ζ c is the value of ζ at the barycenter of c and δ c (q) := |||GRAD c (w)||| λ,c and w := (ζζ c )q. The upper bound in (H1), the definition of GRAD c and that of the |||•||| 2,Ec -norm yield

Figure 5 :

 5 Figure 5: Polyhedral meshes

2 Figure 7 :

 27 Figure7: Test case 6.2: Convergence curves for the two error on the four mesh families using Péclet-based upwinding for λ = 1 (solid lines) and λ = 10 -4 (dashed lines); half-, first-, and second-order slopes are indicated.

2 L 1 (ρ c λ ,c h 2 c |p| 2 H 2 2 . 2 L 1 h 2 c |p| 2 H 2

 212222122 EF)c λ (R Ec (∇p) (EF)c = c∈C GRAD c (q), R Fc (λ•(∇p -G c )) (EF)c -c∈C GRAD c (q), H (EF)c λ (R Ec (∇p -G c )) (EF)c ,and we denote by T 1,1 , T 1,2 the two terms in the right-hand side. The Cauchy-Schwarz inequality, mesh regularity, and the lower bound in (H1) imply that|T 1,1 | |||GRAD(q)||| λ   c∈C e∈Ec ρ c λ ,c h -1 c ||∇p -G c ||Similarly, the Cauchy-Schwarz inequality for H(EF)c λ (i.e., g 1 , H (EF)c λ (g 2 ) (EF)c ≤ |||g 1 ||| λ,c |||g 2 ||| λ,c for all g 1 , g 2 ∈ E c) and the upper bound in (H1) imply that|T 1,2 | |||GRAD(q)||| λ   c∈C e∈Ec λ ,c h c ||∇p -G c ||

  .2b) Remark 2.1 (Boundary term). We have chosen the above definition of DIV since it is naturally associated with the dual mesh where dual cells attached to boundary vertices are, by definition, not closed by dual faces. An alternative choice is to modify the definition of DIV to include dual boundary faces attached to primal boundary vertices; by doing so, the boundary term appears in (2.1) and no longer in (2.2b).

10 Graphe bilan

	Advection Scalaire	Instationnaire	S1
	Elliptique	Eq. Di usion-Advection scalaire Inst.
				S2

  and it depends linearly and monotonically on α (i.e., α ≥ α a.e. in ∂Ω implies that H ∂ α ≥ H ∂ α in the sense of quadratic forms), so that whenever

  Let p be the unique solution of (1.1) and let p be the unique solution of (4.19).

	.20)
	with = 1 2 . Consequently, (4.19) is well-posed.
	Proof. Combine Lemma 3.2 with Lemma 4.2.
	Lemma 4.6 (Error bound).

  e)) h e + |β•n| ,e )h e min(1, Pe e )|p| 2 with |β•n| ,e := ||β•n f (e) || L ∞ ( f (e)) and |p| H 1+ (d(e)) = |p| H 1 (d(e)) + h e |p| H 2 (d(e)) .Proof. The bound on the diffusion-related terms derived in Theorem 4.4 still holds. For the advectionrelated terms, there are two adaptations from the proof of Theorem 3.4. The first one is that we consider R V (p) in lieu of R V (p) since we are now bounding the error (p -R V (p)). The approximation property (3.21),

			1 2
	H 1+ (d(e))		,	(4.22)

  and [[q]] e is defined in Lemma 3.1. Since 2{ {ζ} } 2 e ≤ 2ζ 2 and [[ζ]] 2 e ≤ (L ζ h e ) 2 , we infer that

	|||ζq||| 2 upw,β ≤ 2ζ 2 |||q||| 2 upw,β + 2	e∈E	L 2 ζ h 2 e { {q} } 2
	to (3.12), |||ζq||| 2 |β•n| ≤ ζ 2 |||q||| 2 |β•n|	e ≤ 1 2	v∈Ve q 2 v , and mesh regularity. Since, owing

e Λ e β e , and we conclude using 0 ≤ Λ e β e ≤ |β| | f (e)|, { {q} } 2

Table 6 :

 6 Features of prism meshes with polygonal basis

	Version 1.0

Table 7 :

 7 Features of checkerboard meshes

	(a) CB4 Mesh
	(b) CB8 Mesh

Table 8 :

 8 Features of Kershaw meshes

Conservative form with anisotropic di usion and variable vector field

	0.1	10 ≠4 10 ≠3 10 ≠2 0.1 10 2 Err 2,V and Err ad,V	H sequence 10 3 10 4 ord er 1 or de r 2 #V	Err 2,V and Err ad,V	10 ≠4 10 ≠3	SkP sequence 10 4 ord er 1 or de r 2 #V	10 5
			10 ≠1					10 ≠1
		Err 2,V and Err ad,V	10 ≠3 10 ≠2	or de r 1 o rd e r 2			Err 2,V and Err ad,V	10 ≠2	or de r 1 o rd er 2
			10 2	10 3	10 4 #V	10 5		10 ≠3	10 3	10 4 #V	10 5
						1

Conservative form with anisotropic di usion and variable vector field

	0.1
	Err 2,V and Err ad,V	10 ≠3 10 ≠2	ord er 1 or de r 2			Err 2,V and Err ad,V	10 ≠3	ord er 1 or de r 2
		10 ≠4					10 ≠4
			10 2	10 3	#V	10 4			10 4	#V	10 5
			10 ≠1	CB sequence		10 ≠1	K sequence
		Err 2,V and Err ad,V	10 ≠3 10 ≠2	or de r 1 o rd e r 2				Err 2,V and Err ad,V	10 ≠2	or de r 1 o rd er 2
			10 2	10 3	10 4 #V	10 5		10 ≠3	10 3	10 4 #V	10 5
							1

Conservative form with anisotropic di usion and variable vector field

	0.1
	Err 2,V and Err ad,V	10 ≠3 10 ≠2	ord er 1 or de r 2			Err 2,V and Err ad,V	10 ≠3	ord er 1 or de r 2
		10 ≠4					10 ≠4
			10 2	10 3	#V	10 4			10 4	#V	10 5
			10 ≠1						10 ≠1
		Err 2,V and Err ad,V	10 ≠3 10 ≠2	or de r 1 o rd e r 2				Err 2,V and Err ad,V	10 ≠2	or de r 1 o rd er 2
			10 2	10 3	10 4 #V	10 5		10 ≠3	10 3	10 4 #V	10 5
							1

Conservative form with anisotropic di usion and variable vector field

	Err 2,V and Err ad,V	10 ≠3 10 ≠2	ord er 1 or de r 2			Err 2,V and Err ad,V	10 ≠3	ord er 1 or de r 2
		10 ≠4					10 ≠4
			10 2	10 3	#V	10 4			10 4	#V	10 5
			10 ≠1						10 ≠1
		Err 2,V and Err ad,V	10 ≠3 10 ≠2	or de r 1 o rd e r 2				Err 2,V and Err	10 ≠2	or de r 1 o rd er 2
			10 2	10 3	10 4 #V	10 5		10 ≠3	10 3	10 4 #V	10 5

  Schwarz inequalities for the summations v∈V ∂

	H sequence 0.2 Exponential boundary layer 0.2 Exponential boundary layer			SkP sequence
	0.2 Exponential boundary layer 0.2 Exponential boundary layer		10 ≠1	10 ≠1
	10 2 10 ≠1 10 ≠3 10 ≠2 10 ≠1 10 2 10 ≠1 Err 2,V and Err ad,V 10 ≠3 10 ≠2 10 ≠1 10 2 10 ≠3 10 ≠2 10 ≠1 Err Err 2,V and Err ad,V 2,V and Err ad,V 10 2 10 ≠3	10 3 ord er 1 or d er 2 #V CB sequence 10 4 #V 10 3 10 4 ord er 1 or d er 2 10 3 10 4 ord er 1 or d er 2 #V 10 3 10 4	10 ≠3 10 ≠2 10 ≠1 Err 2,V and Err ad,V 10 ≠3 10 ≠2 10 ≠2 ad,V Err 2,V and Err ad,V 2,V and Err Err 10 ≠3	10 4 order .5 ord er 1 or de r 2 #V order .5 10 4 ord er 1 or de r 2 #V K sequence 10 4 order .5 ord er 1 or de r 2 #V	10 5 10 5	10 5
	10 ≠2 10 3 10 3 10 ≠1 Err 2,V and Err ad,V 10 ≠2 10 ≠2 Err Err 2,V and Err ad,V 2,V and Err ad,V	10 3	10 4 orde r .5 or de r 1 o rd er 2 #V #V 10 4 10 4 #V orde r .5 or de r 1 o rd er 2 orde r .5 or de r 1 o rd er 2	10 5 10 5	10 5	10 3 10 3 10 ≠2 10 ≠2 10 ≠1 Err 2,V and Err ad,V 10 ≠2 10 ≠1 10 ≠1 Err Err ad,V 2,V and Err 2,V and Err ad,V	10 3	10 4 #V #V 10 4 10 4 #V ord er .5 or de r 1 o rd e r 2 ord er .5 or de r 1 o rd e r 2 ord er .5 or de r 1 o rd e r 2	10 5 10 5	10 5
						2	2
						2

f |||g c ||| λ,c ≤ c N |||q ∂ ||| λ/h |||g||| λ , (7.1)

where we have used the local assembly of N ∂ λ on the first line (with c = c(f )), the discrete Cauchy-

f and f ∈F ∂ , and the fact that f ∈F ∂ |||g c ||| 2 λ,c ≤ |||g||| 2 λ on the second

). Since A β (q, q) ≥ 1 2 (|||q||| 2 upw,β + |||q||| 2 |β•n| ) owing to (3.16), we infer that

and the conclusion follows using Young's inequality.

Advection-diffusion

As in Section 4.3, we consider the Péclet-dependent upwinding parameters Λ ve = Θ(ι f (e),c(v) Pe e ) under assumptions (Θ1)-(Θ2). Recall that (Θ1) implies (Λ1). Assume (H1) and (N1) for the diffusion-related terms. Then, provided η 0 ≥

being the maximum number of edges touching a mesh vertex, the following holds:

)

Proof. Set r := ζq + c 4 q. Since (Λ1) holds, we infer from Lemma 5.2 that (I3) holds. Moreover, since c 4 ≥ c 3 ω β and A β (q, q) ≥ 0, Lemma 5.1 implies that

(5.4)

Moreover, owing to (7.2) and to Lemma 5.4 below, we infer that

where we have set g = GRAD(q). Using Young's inequality for the third term on the right-hand side yields

Using the same quadratic identity as in Lemma 4.5, this time with γ = 1 2 c N (ζ + c 4 ) and δ = η 0 (1 + c 4 ), and observing that the choice η 0 ≥

Combining this bound with (5.4) yields A β,λ (q, r) ≥ 1 8 |||q||| 2 ad,V . We conclude using |||r||| ad,V ≤ |||ζq||| ad,V + c 4 |||q||| ad,V and |||ζq||| ad,V ≤ max(

Remark 5.3 (η 0 ). The lower bound for η 0 obtained in Lemma 5.3 slightly differs, up to a numerical factor, from that obtained in Lemma 4.2 for zero advection; the reason is that both proofs have not been optimized regarding the lower bound in the quadratic identity.