Manuel Bernard

Ghislain Lartigue
email: ghislain.lartigue@coria.fr

Guillaume Balarac

Vincent Moureau

Guillaume Puigt

A framework to

Keywords: Deconvolution, Node-Centered Finite Volume, Simplical mesh, High-order approximation, Polynomial reconstruction, Compact stencil

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Motivation and objectives

The Finite Volume Method (FVM) has become a major tool in CFD in the last decades, mostly because of its unique feature of naturally enforcing exact conservation at the discrete level. In this framework, the quantities of interest are the average of the transported variables over Control Volumes (CV). As a consequence, pointwise quantities are not directly available. Traditionally, the direct use of the cell-averaged quantities in the discretization leads to a first-order accurate scheme. But such a scheme is not accurate enough for most industrial and research applications. Following the procedure proposed by Van Leer [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method[END_REF], a second-order accuracy can be recovered if the flux integral over the control-volume surface is not directly expressed from the cell-centered data: a piece-wise linear extrapolation, eventually coupled with a slope limiter, is mandatory. Obtaining any other scheme accuracy is possible if the error between the exact solution and the numerical approximation at any point (or integrated over a CV) is of sufficient order. This numerical approximation can be estimated by combining and integrating several Taylor expansions. This error estimation is illustrated hereafter by considering the transport equation of a scalar φ at a constant velocity u in conservative form:

∂φ ∂t + ∇ • (φu) = 0. (1)
The finite-volume formulation of the problem is obtained by integrating this partial derivative equation (PDE) over a generic control-volume Ω I and by commuting the volume averaging and time differentiation operators (only time-independent meshes are considered here):

dφ

Ω I dt + 1 V I Ω I ∇ • (φu) dV = 0 (2)
The measure of the control volume is defined as:

V I = Ω I dV, (3)
which represents either a volume in 3D, an area in 2D or a length in 1D. The control-volume averaging operator is then:

φ Ω I = 1 V I Ω I φ dV (4)
The divergence of the convective flux in Eq. 2 may be rewritten as an integral over the control-volume boundary ∂Ω I thanks to the Gauss theorem: dφ

Ω I dt + 1 V I ∂Ω I φu • dS = 0 (5)
The exact value of the flux on the boundary ∂Ω I is then estimated numerically and the quality of this approximation depends on the choice of the numerical scheme. The important fact here is that Eq. 5 is an evolution equation for the control-volume average φ Ω I rather than for the nodal value φ I or for any other pointwise value.

However, higher-order schemes require higher-order surface integrals than need either higher-order quadrature rules or moments. In any case, pointwise values of φ or moments based on Taylor series expansions must be defined.

For example, the second-order integral approximation only needs the evaluation of φ at the face center. Even this very simple formulation highlights the necessity to recover pointwise values of φ from the integrated quantities φ Ω I .

This process is known as deconvolution.

Deconvolution is an important step in building of high-order Finite Volume approximations. With the standard k-exact high-order reconstruction, the computation of the polynomial coefficients for approximating the unknowns is possible through the introduction of an extended stencil: the stencil must contain enough control volumes to recover an invertible system where polynomial coefficients are the unknowns. This process is not new: introduced by Barth & Frederickson [START_REF] Barth | Higher order solution of the euler equations on unstructured grids using quadratic reconstruction[END_REF], it was for instance considered by Abgrall [START_REF] Abgrall | On essentially non-oscillatory schemes on unstructured meshes: analysis and implementation[END_REF], Ollivier-Gooch et al. [START_REF] Ollivier-Gooch | Quasi-eno schemes for unstructured meshes based on unlimited data-dependent least-squares reconstruction[END_REF][START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Ollivier-Gooch | Obtaining and verifying high-order unstructured finite volume solutions to the Euler equations[END_REF][START_REF] Sharbatdar | Smoothed truncation error in functional error estimation and correction using adjoint methods in an unstructured finite volume solver[END_REF], Mandal et al. [START_REF] Mandal | High resolution finite volume computations on unstructured grids using solution dependent weighted least squares gradients[END_REF], ... Generally, the stencil used to compute the scheme contains more control volumes than the required number of polynomial coefficients and the computation of the polynomial coefficients is performed by means of a least-squares method over a larger stencil. Recently, Sejekan et al. [START_REF] Sejekan | Improving finite-volume diffusive fluxes through better reconstruction[END_REF] proposed an improvement of the least-squares method based on the H 1 projection instead of the L 2 projection in order to derive an accurate gradient usable for the viscous flux computation. The Monotone Optimal Order-Detection (MOOD) technique [START_REF] Clain | A high-order finite volume method for systems of conservation laws -Multidimensional Optimal Order Detection (MOOD)[END_REF][START_REF] Diot | Improved detection criteria for the Multi-dimensional Optimal Order Detection (MOOD) on unstructured meshes with very high-order polynomials[END_REF][START_REF] Diot | The Multidimensional Optimal Order Detection method in the threedimensional case: very high-order finite volume method for hyperbolic systems[END_REF] follows the same principle for defining local polynomial approximation and mostly differs with k-exact approximation by the polynomial basis.

The definition of a large stencil in a massively parallel framework is the main bottleneck to maintain high performance and an alternative approach is the Successive-Correction k-exact approximation. Introduced by Brenner, Haider, Courbet, Croisille and co-authors for cell-centered finite-volume discretization in several papers [START_REF] Haider | Discrétisation en maillage non structuré général et applications LES[END_REF][START_REF] Haider | Efficient implementation of high order reconstruction in finite volume methods[END_REF][START_REF] Haider | Parallel implementation of k-exact finite volume reconstruction on unstructured grids[END_REF][START_REF] Pont | Self adaptive turbulence models for unsteady compressible flows Modèles de turbulence auto-adaptatifs pour la simulation des écoulements compressibles instationnaires[END_REF][START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible RANS-LES simulations on fully unstructured grids[END_REF][START_REF] Haider | A high-order interpolation for the finite volume method: The coupled least squares reconstruction[END_REF][START_REF] Menasria | Improving the treatment of near-wall regions for multiple-correction k-exact schemes[END_REF], the procedure only needs the smallest possible stencil (local cell and direct neighbors by the faces). Starting from the definition of a discrete gradient operator, the procedure consists in analysing the truncation error of the gradient in order to correct it and to increase its accuracy. Then, the discrete gradient operator is applied to the gradient in order to obtain second-order derivatives. Here again, an analysis of the truncation error of the second-order derivatives enables to define a more accurate polynomial representation of the unknowns. In practice, the stencil is reduced to the minimum but the number of fields to exchange in a parallel environment increases. To the author's knowledge, such a correction procedure was never applied to node-centered formulations.

The Compact Least-Squares Finite-Volume (CLSFV) method [START_REF] Wang | Compact high order finite volume method on unstructured grids I: basic formulations and one-dimensional schemes[END_REF][START_REF] Wang | Compact high order finite volume method on unstructured grids II: extension to two-dimensional Euler equations[END_REF][START_REF] Wang | Compact high order finite volume method on unstructured grids III: Variational reconstruction[END_REF][START_REF] Zhang | Compact high order finite volume method on unstructured grids IV: Explicit multi-step reconstruction schemes on compact stencil[END_REF] gives a framework to define a high-order polynomial of the unknown but it differs with k-exact and MOOD methods by the procedure to build the approximation. Here, starting from the smallest stencil (current cell and its direct neighbors by faces), it is required that the cell averages of the reconstruction polynomial and its various orders of spatial derivatives are conserved on the stencil. This leads to inverting a linear system of equations.

In the context of cell-centered formulations, Betchen and Straatman [START_REF] Betchen | An accurate gradient and hessian reconstruction method for cell-centered finite volume discretizations on general unstructured grids[END_REF] define a new reconstruction method to deal with accurate gradient and Hessian. Their idea is simple: a third-order accurate representation of unknowns needs a second-order accurate gradient and a first-order accurate Hessian. To obtain these operators, Betchen and Straatman write a system of equations valid over the whole computational mesh and invert it using an underrelaxed block-Jacobi method. Such a method does not cope with our requirements for HPC since inverting the system cannot be performed locally, by a cell-wise formulation.

The aim of this work is to present a framework that can provide a third-order polynomial reconstruction of a field from CV-averaged quantities on any simplicial mesh. From this polynomial reconstruction, it is then straightforward to find accurate approximations of the field, but also of its gradients and its Hessian. It will be shown that it is possible to achieve a second-order accuracy on the field by only considering the direct neighbors of the mesh vertices. Moreover, with this same connectivity, it is possible to obtain a third-order accuracy in most situations, with the noticeable exception of boundaries. Obtaining third-order accuracy with such a compact stencil is clearly an advantage in the context of HPC applications, where parallel performance is an issue. Moreover, it is possible to extend the present method to recover third-order accuracy near the interfaces. Finally, the extension of the present method to higher-order reconstructions is straightforward from a theoretical point of view. However, the computational cost and the mathematical complexity associated to higher-order methods can become difficult to handle and for this reason, only second-and third-order methods will be presented here.

Another important issue that will be addressed in this work is the quality of the approximations of the classical FV discrete operators on a generic (i.e. not regular) mesh. We will show that in node-centered codes, most discrete gradient operators are not consistent or first-order only in the best case. Similarly, most discrete Hessian operators will usually diverge when the mesh is refined. The present framework proposes to overcome these limitations by exhibiting discrete operators that provide a second-order accuracy on gradient and a first-order accuracy for Hessian on any mesh.

The applications of this deconvolution technique for Node-Centered Finite-Volumes codes are numerous and only a few of them will be presented here to illustrate the ability of this framework to achieve high-order accuracy in any situation:

• High-precision interpolation

• Curvature computation of a levelset function for two phase flows simulations • High-order numerical schemes for convection.

2 A framework for deconvolution of node-centered FVM

Taylor series expansion and grid moments

Assuming a sufficiently smooth scalar field, any value φ(x) can be approximated anywhere in the domain by performing a Taylor series expansion around any node I:

φ(x) = φ I + (∇φ I ) • ∆ I (x) + 1 2 (∇∇φ I ) : (∆ I (x) ⊗ ∆ I (x)) + O(∆ 3 ) (6)
where ∆ I (x) = xx I . The ⊗ and ∇ symbols have their usual meaning: i.e. the tensor product and the gradient operator, respectively. Thus ∇φ I is the gradient of φ estimated at x I and ∇∇φ I is the Hessian of φ, also estimated at x I .

In the context of the FVM, the averaging operator defined in Eq. 4 is then applied to this expression to provide a first relation between the nodal and the averaged value of φ:

φ

Ω I = φ I + (∇φ I ) • ∆ I Ω I + 1 2 (∇∇φ I ) : ∆ I ⊗ ∆ I Ω I + O(∆ 3 ) (7)
where

∆ I Ω I = 1 V I Ω I (x -x I ) dV = O(∆) (8)
∆ I ⊗ ∆ I Ω I = 1 V I Ω I (x -x I ) ⊗ (x -x I ) dV = O(∆ 2 ) (9)
are the first-and second-order moments of the control volume, respectively. These quantities are of key importance in the present method and thus deserve a few comments. First, it can be noticed that ∆ I Ω I is simply the vector between the vertex I and the mass center of the control-volume. It is thus zero when I is the center of the control volume. Similarly, the ∆ I ⊗ ∆ I Ω I term is a symmetric rank 2 tensor which is never zero: its eigenvalues are the moments of inertia of the control volume with respect to I. Then, it should be emphasized that in Eq. 7, both the Taylor expansion and the averaging operation are performed around the same position I. This is indeed a very special case and by no means a requirement: for example, one can also define ∆ I Ω J by integrating the function

∆ I (x) over a different control volume Ω J .
As already mentioned, it is straightforward to extend Eq. 7 to higher order by defining tensors

∆ ⊗3 I Ω J , ∆ ⊗4 I Ω J ,
and so on. Of course, the amount of data associated to these tensors increases substantially for higher orders but it can also be noted that these tensors are geometric quantities which can be pre-computed and stored once and for all: the CPU cost associated to their evaluation is thus not an issue, provided that they can be stored in memory.

Moreover, these high-rank tensors have the maximal symmetry property: the amount of data which must be stored to describe them remains tractable for practical applications.

Some technical details presenting the computation of these moments are gathered in Appendix A.

Control-Volumes, mesh quality and discrete differential operators

A simple inspection of Eq. 7 shows that the evaluation of the discrete differential operators at node I and the choice of the control-volume around this node will have a major impact on the deconvolution procedure. The aim of this section is to clarify these two points.

In the context of the Node-Centered approach, the CV are the cells of the dual mesh. Several methods have been described extensively in the literature to build a dual mesh from the primal mesh (barycentric, Voronoï, . . . ) [START_REF] Mullen | Hot: Hodge-optimized triangulations[END_REF] so this choice must be clearly stated. Here, the barycentric method is used, but the results can be reformulated for other dual meshes without any major difficulty. The case of non-simplicial meshes will not be discussed here as it is always possible to decompose any element into a finite number of simplices in a pre-processing step: the proposed framework can thus be applied to any grid, provided that such a decomposition has been applied beforehand.

Once the primal and dual meshes have been given, it is then possible to introduce the concept of mesh quality:

for a given number of elements, a regular mesh will usually provide a better accuracy than a grid containing highly distorted elements. This concept is detailed in Appendix B where some notations and some objective criteria to quantify the lack of regularity of a given grid are introduced.

Many discretizations of the usual differential operators have been proposed in the context of FVM [START_REF] Barth | A 3-d upwind euler solver for unstructured meshes[END_REF][START_REF] Turner | A comparison of gradient approximations for use in finite-volume computational models for two-dimensional diffusion equations[END_REF][START_REF] Hyman | High order finite volume approximations of differential operators on nonuniform grids[END_REF].

Among them, the following pair-based formulations for the nodal gradient and Hessian are easily derived by applying the Gauss Theorem to the dual mesh and are often used in CFD codes, either in the simple form presented hereafter or in more complex forms to define the interface quantities [START_REF] Puigt | Discretisation of diffusive fluxes on hybrid grids[END_REF]. Indeed, the Gauss theorem needs the definition of an interface quantity that is assumed here to be an average of left and right contributions:

G I (φ) = 1 V I J∈N I φ I + φ J 2 × S IJ (10)
or an interface gradient that is approximated by a standard finite difference:

H I (φ) = 1 V I J∈N I φ J -φ I IJ × 1 2 (n IJ ⊗ S IJ + S IJ ⊗ n IJ ) (11)
where the notations refer to Fig. 1 and where N I is the set of vertices connected to I by an edge of the mesh, i.e.

the direct neighbors of I. Both formulations can be found in [START_REF] Barth | Aspects of unstructured grids and finite-volume solvers for the euler and navier-stokes equations[END_REF], except that we have slightly modified Eq. 11 to ensure the mixed partial derivatives theorem. By noting V the set of all vertices of the mesh, φ = (φ I ) I∈V is a vector containing the values evaluated at each node of the mesh. This quantity can either be pointwise (e.g. φ = (φ I ) I∈V or ∆ J = (∆ J (x I )) I∈V ) or an average value (e.g. φ Ω = φ

Ω I I∈V or ∆ Ω J = ∆ Ω I J I∈V
). A special attention must be paid when dealing with moments such as ∆ J and ∆ Ω J : the index J refers to the center of the function and must be kept constant in the gradient/Hessian operators. For the sake of simplicity, everything is presented here for a scalar φ but the method can be extended trivially to any tensor. IJ and s 2 IJ . In 3D, the situation is similar but the number of facelets based on an edge is not fixed (generally 5 to 8 facelets by edge). The vector n IJ represents the unit vector direction on edge IJ.

I J K L E C 1 C 2 F 2 F 1 n IJ s 2 IJ s 1 IJ S IJ
The discrete nodal gradient of Eq. 10 can achieve second-order accuracy on regular meshes, but is only first-order on generic meshes as will be demonstrated below. It can also be highlighted that the adaptation of Eq. 10 to the cell-centered Finite-Volume approach is not straightforward: the mesh interface being generally not located exactly between cell centers, the interface quantity cannot be defined by a simple averaging procedure [START_REF] Sozer | Gradient calculation methods on arbitrary polyhedral unstructured meshes for cell-centered CFD solvers[END_REF]. Moreover, the approximation of the interface gradient cannot be considered in standard cell-centered Finite-Volume approach since the line passing between cell centers does not generally cross the interface in its center. Similarly, it will also be shown that the approximation of the Hessian given in Eq. 11 usually does not converge on distorted meshes. As a consequence, these operators can not be used in Eq. 7 to achieve a third-order reconstruction of the solution on distorted meshes. This clearly highlights the necessity to improve the calculation of these two differential operators to perform the high-order deconvolution of the average solution.

Third-order deconvolution procedure

In this section, we detail the methodology that enables to express high-order approximations of φ I , ∇φ I and ∇∇φ I from φ Ω . To achieve this goal, some discrete operators for the nodal gradient G I and the nodal Hessian H I must already be available, even if they are not high-order. Their exact expression is actually not very important and we will thus consider generic linear operators written as:

G I φ Ω = J∈N I φ Ω J P 1 IJ (12)
H I φ Ω = J∈N I φ Ω J P 2 IJ ( 13
)
where N I = N I {I} is the set of direct neighbors of I augmented with I itself and P 1 IJ and P

2

IJ are generic rank 1 and rank 2 tensors. One can notice that both operators are compact, in the sense that they only imply the direct neighbors of a given node. In the present paper, we used the approximations presented in Eq. 10 and Eq. 11 but other choices would have provided similar results.

To be consistent, it is important that the discrete Hessian satisfies the Schwartz theorem on mixed partial derivatives which simply requires that P 2 IJ should be symmetrical. Moreover, in what follows, we will assume that:

P n = O ∆ -n for n ∈ {1, 2} (14)
which is rather natural to ensure the operators consistency. Additional consistency conditions will be discussed thoroughly hereafter.

Applying the volume averaging, discrete gradient and discrete Hessian operators to the Taylor expansion of φ (Eq. 6), and using the linearity of these operators, leads to:

φ Ω I = φ I 1 Ω I + (∇φ I ) • ∆ I Ω I + 1 2 (∇∇φ I ) : ∆ I ⊗ ∆ I Ω I + O ∆ 3 , (15)
G I φ Ω = φ I G I 1 Ω + (∇φ I ) • G I ∆ I Ω + 1 2 (∇∇φ I ) : G I ∆ I ⊗ ∆ I Ω + O ∆ 2 , (16)
H I φ Ω = φ I H I 1 Ω + (∇φ I ) • H I ∆ I Ω + 1 2 (∇∇φ I ) : H I ∆ I ⊗ ∆ I Ω + O (∆) , (17)
It can be seen that some new high-rank tensors are involved in Eq. 16 and Eq. 17. These tensors are built by applying the discrete gradient and Hessian operators to the vectors ∆ I Ω and ∆ I ⊗ ∆ I Ω . From a practical point of view, and thanks to Eq. 12 and Eq. 13, it is sufficient to evaluate ∆ I Ω J and ∆ I ⊗ ∆ I Ω J for J ∈ N I , and then to compute some linear combination of these quantities at each node of the mesh.

Eq. 15, Eq. 16 and Eq. 17 can be rewritten all together as:

L I = C O3 I • H I + E (18)
where:

L I =     φ Ω I G I φ Ω H I φ Ω     , H I =   φ I ∇φ I ∇∇φ I   and E =   O ∆ 3 O ∆ 2 O (∆)   (19)
are the vectors of low-order approximations, exact nodal values and errors respectively. The 3 rd -order convolution matrix C O3 I is then:

C O3 I =     1 Ω I ∆ I Ω I 1 2 ∆ I ⊗ ∆ I Ω I G I 1 Ω G I ∆ I Ω 1 2 G I ∆ I ⊗ ∆ I Ω H I 1 Ω H I ∆ I Ω 1 2 H I ∆ I ⊗ ∆ I Ω     (20)
The previous notations rely on a consistent unfolding of the various high-rank tensors that appear in Eq. 19 and

Eq. 20 to form a classical linear system. This unfolding is actually trivial and will not be commented any further here, except for the fact that it must take into account the various symmetries of the tensors to ensure that the matrix C O3 I has the maximal-rank property.

This matrix is of key importance in the present work and deserves some additional comments. First it must be noted that the ideal situation is the case where C O3 I = I: in this case the discrete and continuous variables are identical up to the theoretical order indicated by the vector E. In this situation, no deconvolution is needed as the discrete operators can be used directly to perform a polynomial representation of the solution in a given control volume. The main point of the present article is to realize that such an ideal situation is never achieved, even on regular meshes. Moreover, on distorted meshes, this matrix can depart greatly from identity and the discrepancy between L and H can become very important. In any case, one has to solve Eq. 18 to obtain high-order approximations of H from the available data, i.e. L.

Considering a single line of the matrix, it can be noticed that the consistency between the discrete and the continuous part is achieved when the elements located at the left of the diagonal of C O3 I are zero and when the diagonal term is one. Assuming that these consistency constraints are fulfilled, the order of the approximation is then given by the position of the first non-zero term on the right of the diagonal.

Several comments can be formulated at this point. First, from the definition of the averaging operator, 1

Ω I = 1,
and thus φ Ω I is always a consistent approximation of φ I . Another desirable (but not mandatory) property is that the discrete gradient and Hessian satisfy the following constraints:

G I 1 Ω = 0 and H I 1 Ω = 0 (21)
This is for example the case for operators given in Eq. 10 and Eq. 11. These conditions state that the discrete gradient and Hessian operators applied to a constant field must be zero, which are rather natural requirements.

Moreover these two requirements imply a very nice property: the reconstructed gradient and hessian are then decoupled from the reconstructed nodal value. As a consequence, it is sufficient to consider a reduced linear system, i.e. without the first line and column of the convolution matrix, to perform the numerical reconstruction procedure.

Finally, the nodal value is obtained in a post-processing step to ensure that the reconstruction conserves exactly the average value. In what follows, the full system is considered for the sake of clarity but numerically, only the reduced linear system is inverted, following the procedure that is described below.

The important special case of regular meshes is now considered. In the present paper, a mesh is said to be regular around vertex I if the condition ∆ J Ω J = 0 holds for all J ∈ N I . In this situation it can be shown that for the discrete operators considered in Eq. 10 and Eq. 11, many terms cancel in the deconvolution matrix and then:

H I φ Ω = ∇∇φ I + O ∆ h (22)
G I φ Ω = ∇φ I + O ∆ 2 (23)
φ

Ω I = φ I + 1 2 H I φ Ω : ∆ I ⊗ ∆ I Ω I ( 24
)
where h = 2 in 1D and h = 0 otherwise. Equation 24shows that φ Ω I is a second-order approximation of φ I in the case of regular meshes.

On the other hand, it can also be shown that if the mesh is not regular those discrete operators are usually not consistent. More specifically, in this situation, G I ∆ I Ω = I 2 and H I ∆ I Ω = 0 3 and finally φ Ω I is only a first-order approximation of φ I .

It is now possible to define the following new operators, simply by inverting the system defined in Eq. 18,

provided that the matrix C O3 I has the full rank property:

  φ O3 I G O2 I H O1 I   = C O3 I -1 •     φ Ω I G I φ Ω H I φ Ω     (25)
where φ O3 Of course it is possible to downgrade the quality of this deconvolution to obtain a second-order approximation of φ and a first-order approximation ∇φ:

φ O2 I G O1 I = C O2 I -1 • φ Ω I G I φ Ω ( 26
)
where the second-order convolution matrix C O2 I is simply:

C O2 I = 1 Ω I ∆ I Ω I G I 1 Ω G I ∆ I Ω (27)
As already shown above, the first-order convolution matrix C O1 I simply reduces to the scalar 1 and thus:

φ O1 I = φ Ω I (28)
is a first-order approximation on a generic mesh.

To conclude, if the convolution matrix C O3 I is invertible, the method described above provides a third-order polynomial expansion of the solution around every CV of the computational domain:

φ O3 I (x) = φ O3 I + G O2 I • ∆ I (x) + 1 2 H O1 I : (∆ I (x) ⊗ ∆ I (x)) = φ I + O ∆ 3 . ( 29
)
It must be noted that this procedure can be carried out analytically in some simple situations: the corresponding results are given in Appendix C.

Invertibility of the convolution matrix

As already mentioned, the convolution matrix C I and the deconvolution matrix C -1 I only depend on geometrical quantities and can thus be computed and stored once for all in a pre-processing step. Their computation is thus transparent from a performance point of view. The only additional cost associated to the deconvolution is thus a matrix-vector product of very limited size, which can be done efficiently on all modern processors.

The main advantage of the present method should now be clear: building a third-order approximation of the solution even on very distorted meshes only requires to compute two compact operators and to invert a small matrix on each node of the mesh. There is no need to build long-distance stencils, which can be both difficult to implement and numerically inefficient for massively parallel solvers. On the other hand, each node must have enough neighbors so that the convolution matrix is invertible. Indeed, the connectivity of a given node needs to be at least equal to the dimension of the convolution matrix (see Tab. 1 that summarizes the size of the various convolution matrices).

For example, the third-order deconvolution matrix in 3D has 10 independent lines and columns: 1 for the field, 3 for the gradient and 6 for the Hessian (symmetric 3 × 3 matrix). It means that the set N I must contain 10 elements, i.e. each node must be connected to 9 neighbors.

If this condition is not respected (for example, the nodes lying on the boundaries usually do not have that many neighbors), several solutions can be proposed:

• it is always possible to perform a 2 nd -order deconvolution at the concerned node: on a simplicial mesh of dimension d, the cells have d + 1 vertices and each vertex is connected to the other d vertices which are thus direct neighbors. This implies a local loss of accuracy, but this is usually not a problem since it concerns a limited number of nodes. Moreover, most numerical methods have an order decrease near the boundaries due to the necessary upwinding.

• another possibility is to use the polynomial reconstruction performed on some neighboring nodes, if they have themselves a sufficient connectivity.

• a last technique could consist in blending low-order local reconstructions into a higher-order one, as performed in the P n -P m method (see [START_REF] Dumbser | A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes[END_REF] as the first reference on this technique).

The reason why the present framework can only be applied to node-centered Finite Volumes should now be clear: in the cell-centered method, the number of direct neighbors is simply the number of faces of the simplex which is always equal to d + 1. As indicated by Tab. 1, this connectivity only allows to achieve a second-order dimension (d) 1D 2D 3D

C O1 I 1 1 1 C O2 I 2 3 4
C O3 I 3 6 10
C On I n n(n+1) 2 n(n+1)(n+2) 6
Table 1: Size of the convolution matrices in 1D, 2D and 3D (taking the various symmetries into account).

accuracy on the scalar field.

Link with the Least-Squares approach

It must be pointed out that there are many similarities between our approach and the Weighted Least-Squares method proposed by Ollivier-Gooch [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF]. Indeed, let us introduce the matrix of the moments around node I :

R I =       1 Ω I ∆ I Ω I 1 2 ∆ I ⊗ ∆ I Ω I 1 Ω J 1 ∆ I Ω J 1 1 2 ∆ I ⊗ ∆ I Ω J 1 . . . . . . . . . 1 Ω J N ∆ I Ω J N 1 2 ∆ I ⊗ ∆ I Ω J N       (30)
It is a matrix of dimension (N + 1) × M where N stands for the number of neighbours of node I and M is the size of the convolution matrix as detailed in Tab. 1. This matrix can be used to write the node-averaged Taylor expansion around the node I:

      φ Ω I φ Ω J 1 . . . φ Ω J N       = R I   φ I ∇φ I ∇∇φ I   + O(∆ 3 ) (31)
We also introduce the low-order operator matrix (notations are taken from Eq. 12 and Eq. 13:

P I =    1 0 . . . 0 P 1 II P 1 IJ1 . . . P 1 IJ N P 2 II P 2 IJ1 . . . P 2 IJ N    (32)
which is used to build the low order operators from neighbouring nodes as:

    φ Ω I G I φ Ω H I φ Ω     = P I       φ Ω I φ Ω J 1 . . . φ Ω J N       (33)
With these notations, our convolution matrix is simply:

C I = P I R I ( 34
)
while the matrix that arises in the Weighted Least-Squares formalism is:

S I = R T I Λ I R I ( 35
)
where Λ I is a diagonal matrix introduced by Barth [5] to lessen the weight of farthest nodes in the reconstruction.

There is thus a clear link between the two approaches: the initial statement (Eq. 31) is identical but the way to solve this over-determined linear system is different. In our framework, we state that:

  φ I ∇φ I ∇∇φ I   = (P I R) -1 P I       φ Ω I φ Ω J 1 . . . φ Ω J N       (36)
while the LS formalism solves for:

  φ I ∇φ I ∇∇φ I   = (R T I Λ I R I ) -1 R T I       φ Ω I φ Ω J 1 . . . φ Ω J N       (37)
A major difference between the two frameworks is the condition number of the matrices C I and S I that need to be inverted to build the reconstruction. By a simple dimensional analysis, it can be seen that the condition number of our convolution matrix C I is dimensionless: it only depends on the mesh quality but not on the mesh size. This is not the case of the S I matrix: its eigenvalues have physical dimensions 1, ∆ 2 , ∆ 4 , . . . This point is of paramount importance from a numerical point of view and, to the authors' understanding, it is a clear advantage of the present method over the classical Least-Square approach.

As an illustration, the Fig. 2 presents the condition number of the convolution matrix for various 1D, 2D and 3D meshes (with very different physical sizes and with large amount of perturbation). It clearly appears that this condition number is better on regular meshes than on distorted meshes but is always of the order of unity even on quite irregular grids.

Both regular and distorted meshes are used to evaluate the performance of the method: those meshes are identified thanks to the notation M d b , where d is the dimension of the problem and b is a parameter that indicates the perturbation of the grid, as explained in Appendix B (b = 0 indicates a regular mesh and the more b increases, the more the mesh is deformed).

In both 2D and 3D, some meshes have also been generated by using dedicated commercial software. Those meshes are identified thanks to the notation M d sof t .

Extension to higher orders

First, the extension of Eq. 6 to the order k is straightforward. Then, some discrete differentials operators also up to order k must be available. Those operators can be obtained by composing the low-order first and second derivatives as many times as needed, for example:

D (odd) = G • H • . . . • H (38)
D (even) = H • H • . . . • H (39)
Of course, this implies that larger and larger stencils will be involved but this is anyway needed to ensure the invertibility of the convolution matrix. Then, one has to build and to invert the convolution matrix C Ok I . For this, the grid moments up to ∆ ⊗k-1 I Ω J must also be pre-computed and the previously introduced differential operators must be applied to them. When k is greater than 3, even if still feasible, this can become challenging from a computing point of view, as it implies to deal with large amount of data for each vertex of the mesh. Moreover, achieving third-order accuracy is usually sufficient for many practical applications.

Numerical verification of the accuracy of the deconvolution

The present section describes the tests that are used to assess the accuracy of the deconvolution procedure described above. In these tests, a 1D bump function is used, i.e. a compactly supported C ∞ function, which depends only on x:

φ(x) = exp 1 R 2 0 exp 1 x 2 -R 2 0 if |x| < R 0 0 else (40)
In all presented results, the domain length is L x = 1 and R 0 = 0.4. The convergence order is studied by varying the number of cells N x in the principal direction. In 2D and in 3D, the domain width/height is always set to 30 cells to keep a reasonable CPU cost on the finest meshes: this is not an issue since the solution mostly remains 1D.

Moreover, additional tests with full 2D and 3D solutions have not shown any noticeable difference.

To apply the deconvolution procedure, the average value of the scalar φ Ω I over each CV of the mesh must be available at the initial time:

• A first possibility to obtain this average field is to compute it analytically by integrating the bump function on each CV. However, depending on the signal expression, this can be cumbersome, especially in 2D and 3D.

• Another option is to compute this field from Eq. 15: this requires to know the analytical expression of φ, ∇φ and ∇∇φ and the first and second order moments. This method will provide a third-order accurate approximation of the analytical average value and is used in the present situation.

• In the general case, when the nodal data is known but not its successive derivatives, it is not a trivial task to compute the nodal gradient up to the second order and the nodal Hessian up to the first order, to ensure a third order reconstruction of the average data.

In this last situation, it is still possible to apply the discrete gradient and Hessian operators to Eq. 6 at point I:

G I (φ) = ∇φ I • G I (∆ I ) + ∇∇φ I : G I (∆ I ⊗ ∆ I ) + O(∆ 2 ) ( 41
)
H I (φ) = ∇φ I • H I (∆ I ) + ∇∇φ I : H I (∆ I ⊗ ∆ I ) + O(∆) (42)
For the previous equations to hold, the assumptions G I (1) = 0 and H I (1) = 0 must be verified (Eq. 21) but this is not an issue as already explained. Following the same path as for the averaged quantities (see Eq. 25), this can be re-written as a linear system:

G O2 I H O1 I = C O3 I -1 • G I (φ) H I (φ) (43)
where the pointwise convolution matrix

C O3 I C O3 I = G I (∆ I ) 1 2 G I (∆ I ⊗ ∆ I ) H I (∆ I ) 1 2 H I (∆ I ⊗ ∆ I ) (44)
does not require the knowledge of the CV moments but only evaluation of the function ∆ I (x) = xx I for all J ∈ N I . The evaluation of this pointwise convolution matrix C O3 I is actually very similar to the evaluation of C O3 I and will not be commented any further.

Once φ Ω I is available (by any of the aforementioned method), the discrete differential operators Eq. 10 and

Eq. 11 are applied to this averaged field and then the high-order approximations φ O3 I , G O2 I and H O1 I are obtained by solving linear system of Eq. 25.

The accuracy of the resulting quantities is then evaluated by computing L 1 -, L 2 -and L ∞ -norms of errors between the theoretical values and the approximations. As expected, the convergence order is independent of the specific norm chosen to measure the error. Thus, for the sake of conciseness, only results corresponding to L 2 -norm are represented hereafter. For example, the error on the deconvoluted scalar field reads:

E(φ) 2 = 1 V Ω NΩ I=1 φ ex I -φ O3 I 2 V I ( 45
)
where N Ω is the number of nodes in the domain Ω. For the differential operators, the L 2 -norm of error is evaluated for each component and then the Euclidean norm and the Frobenius norm are computed for gradient and Hessian, respectively.

Fig. 3a shows that the classical Hessian operator converges in 1D whereas it is either O (1) on regular meshes or even O ∆ -1 on distorted 2D and 3D meshes. More precisely:

• the error on the diagonal terms is O(∆ 2 ) on regular meshes and O(∆) in the general case.

• the error on the off-diagonal terms is O(1) on regular meshes and O(∆ -1 ) in the general case.

That explains the specific behavior observed on the 1D results, as there is no off-diagonal terms in this situation.

Regarding the deconvoluted Hessian, Fig. 3b shows that it is O(∆ 2 ) on regular meshes and O(∆) on distorted meshes, whatever the problem dimension.

Concerning the approximation of the gradient, Fig. 4 shows that, even if the classical operator is O(∆ 2 ) on regular meshes, its convergence decreases to O(∆) in the general case. Here again, the deconvoluted gradient remains O(∆ 2 ), no matter the grid quality.

Finally, as illustrated in Fig. 5, the averaged scalar φ Ω I approximates the pointwise value up to O(∆ 2 ) on regular meshes and O(∆ 1 ) on distorted meshes. By using the new framework, the pointwise reconstruction of the scalar is O(∆ 4 ) and O(∆ 3 ) on regular and distorted meshes respectively.

Application to high-order interpolation

The availability of a high-order polynomial reconstruction (Eq. 29) offers the possibility to perform an accurate interpolation of a field at any location in the domain. It permits for example to interpolate precisely any fluid properties from an Eulerian phase towards Lagrangian particles. Another obvious application of this high-order deconvolution framework is to perform accurate mesh-to-mesh interpolation in the context of dynamic mesh adaptation [START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF] for example. This last point is of particular importance: dynamic mesh adaptation implies a very frequent remeshing and any spurious numerical noise can accumulate over time and produce unexpected discrepancies if low-accuracy interpolation is used. In the context of dynamic mesh adaptation, a conservative and monotonicity preserving interpolation was shown to overcome standard interpolation technique [START_REF] Alauzet | p 1 -conservative solution interpolation on unstructured triangular meshes[END_REF] but is still limited to second-order of accuracy and work is still needed to consider higher accuracy. For completeness, it must also be mentioned that an alternative to interpolation with polynomials is interpolation using the Radial Basis Functions. The literature is large but recent applications concern the treatment of non conforming grid interface for high-order finite-volume schemes [START_REF] Bras | A technique of flux reconstruction at the interfaces of non-conforming grids for aeroacoustic simulations[END_REF] or the coupling computations using non conforming grids and potentially different discretization techniques [START_REF] Gervasio | The INTERNODES method for non-conforming discretizations of PDEs[END_REF].

M 1 0.1 M 1 0.2 M 1 0.4 8 4 0 -4 -8 -12 log ( E ) 2D M 2,F K 0 M 2,F K 0.01 M 2,F K 0.05 M 2,F K 0.2 M 2,iso 0 M 2,iso 0.35 M 2 sof t
M 1 0.1 M 1 0.2 M 1 0.4 0 -4 -8 -12 -16 log ( E ) 2D M 2,F K 0 M 2,F K 0.01 M 2,F K 0.05 M 2,F K 0.2 M 2,iso 0 M 2,iso 0.35
M 1 0.1 M 1 0.2 M 1 0.4 -8 -12 -16 -20 -24 -28 -32 -36 log
( E ) 2D M 2,F K 0 M 2,F K 0.01 M 2,F K 0.05 M 2,F K 0.2 M 2,iso 0 M 2,
The interest of the proposed framework is illustrated by applying the high-order interpolation method proposed by Délèze et al. [START_REF] Délèze | Finite elements of c1-class on a tetrahedron[END_REF][START_REF] Délèze | Eléments finis tétraédriques de classe C1 et de degré deux[END_REF]. This method is an extension of the Hermite interpolation method to 2D and 3D simplices which relies on the construction of a local finite-element basis in each cell of the mesh. This leads to interpolation functions which are C 1 -differentiable everywhere, and C 2 -differentiable almost everywhere and that can achieve third-order convergence order. It is particularly interesting as it ensures the continuity of the fluid properties (and of their first derivative) when Lagrangian particles cross an element boundary. To achieve these desirable properties, the Délèze method requires to have the value of the field and its gradient at each vertex of the mesh with respectively 3 rd and 2 nd order accuracy.

To illustrate this application, the bump function presented in section 3.1 is interpolated from a source mesh consisting of a manufactured distorted grid using perturbation parameter b = 0.2 (see Appendix B) towards a regular destination grid, i.e. b = 0. Both grids contain the same number of nodes (which do not coincide) and the convergence order is studied by refining the source and the destination grid simultaneously.

Fig. 6 shows the L 2 -norm of the error on the solution obtained on the destination grid by the interpolation method proposed by Délèze. This figure shows three different datasets:

• The curve labeled analytical is obtained by using the analytical values for the data and its gradient in the formula of Délèze. As expected, a 3 rd order accuracy is obtained.

• The low-order curve corresponds to the results obtained when replacing the analytical values of φ I and ∇φ I by the averaged data φ Ω I and its low-order gradient G I φ Ω . It can be seen than in this situation, the method degenerates to 1 st -order.

• Finally, the high-order set is obtained by using the deconvoluted values φ O3 I and G O2 I : in this situation, the 3 rd -order theoretical convergence is recovered, with however a lower precision than the first method which uses the analytical values.

Application to levelset curvature computation

In many multiphase applications, the interface is treated by the mean of the levelset method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamiltonjacobi formulations[END_REF]: in this case, an accurate evaluation of geometrical quantities is requested to ensure the consistency of the simulation. For example, computing the pressure drop between the two phases with Young-Laplace equation requires to evaluate the curvature of the interface. However, computing a consistent curvature on a distorted mesh is a common issue in multiphase flows [START_REF] Desjardins | An accurate conservative level set/ghost fluid method for simulating turbulent atomization[END_REF][START_REF] Chiodi | A reformulation of the conservative level set reinitialization equation for accurate and robust simulation of complex multiphase flows[END_REF][START_REF] Marchandise | A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics[END_REF].

In what follows, the levelset function φ can represent either the hyperbolic tangent or the signed distance function depending on the considered levelset method. In both cases, the normal to the interface is defined as: and the the total curvature of the interface, i.e. the sum of the main curvatures, reads:

n = ∇φ ∇φ . ( 46
)
K = -∇ • n. (47)
As summarized by Goldman [START_REF] Goldman | Curvature formulas for implicit curves and surfaces[END_REF], various formulations are available to evaluate numerically this curvature. All these formulations are equivalent when considering continuous differential operators but can show significant differences when evaluated numerically. In this paper we focus on two formulations.

The first one is a two-step method, based on the preliminary computation of the normal following Eq. 46. Then, its divergence is computed as the trace of the gradient of the normal:

K nc = -tr(∇n) = -tr ∇ ∇φ ∇φ (48)
This method requires two applications of the gradient operator to the levelset function and is thus non-compact (as indicated by the subscript nc). The second formulation is expressed as a function of the gradient and the Hessian of the levelset function.

K c = ∇∇φ ∇φ : ∇φ ⊗ ∇φ ∇φ 2 -I 2 (49)
which requires to evaluate the discrete gradient and Hessian to compute the curvature but only implies the direct neighbors: this method is thus named compact, as indicated by the subscript c.

Let us point out that with both formulations, due to the discrete operators definition, the curvature is known at the mesh vertices. On the other hand, most numerical methods that solve the Laplace equation on the pressure need to evaluate the pressure drop at interface position, which generally does not coincide to a node position.

In the case of classical FVM framework, regardless the formulation employed (Eq. 48 or Eq. 49), the curvature is linearly interpolated at interface location Γ with the usual formula:

K Γ = θ IJ K I + (1 -θ IJ )K J , (50)
where:

θ IJ = ∆ I (x Γ ) ∆ I (x J ) , (51)
In this last equation, ∆ I (x Γ ) is the vector between node I and interface Γ and ∆ I (x J ) is the vector between node I and node J.

When using the new high-order FVM framework, the deconvoluted gradient and Hessian are first evaluated on both sides of each edge [IJ] that crosses the interface, i.e. at node I and J, with Eq. 25.

Then their values at the interface position Γ are computed by performing a Taylor series expansions with respect to I and J:

   G O2 Γ,I = G O2 I + ∆ I (x Γ ) • H O1 I + O(∆ 2 ) (52a) G O2 Γ,J = G O2 J + ∆ J (x Γ ) • H O1 J + O(∆ 2 ) (52b)    H O1 Γ,I = H O1 I + O(∆) (53a) H O1 Γ,J = H O1 J + O(∆) (53b)
Finally, the values G O2 Γ and H O1 Γ that are used in Eq. 48 or in Eq. 49 to evaluate the curvature at the interface are computed with the following linear interpolation:

   G O2 Γ = θ IJ G O2 Γ,I + (1 -θ IJ )G O2 Γ,J (54a)
H O1 Γ = θ IJ H O1 Γ,I + (1 -θ IJ )H O1 Γ,J (54b)
To illustrate the accuracy of both formulations, the curvature of a spherical bubble of radius R is computed for various mesh resolutions both in 2D and in 3D, for regular (b = 0) and distorted (b = 0.2) grids. Note that the coarser mesh N x = 8 corresponds to R 2.66 ∆ x . In each situation, the error is computed with the classical and the deconvoluted discrete operators presented in section 2 and the results are gathered in Fig. 7a and Fig. 7b. More precisely, these figures present the L 2 -norm of the relative error on curvature which is evaluated by comparing the discrete value along each edge to the analytical result K 0 = (d -1)/R:

E(K) 2 = 1 N e Ne e=1 K i -K 0 K 0 2 (55)
where N e is the number of edges of the mesh that intersect the interface.

As illustrated in Fig. 7a, the non-deconvoluted differential operators can provide a O(∆ 2 ) convergence order on regular meshes when combined with the divergence formula (Eq. 47). But when the same methodology is applied to a distorted grid, the numerical evaluation of the curvature fails to converge and is O(1). The compact formulation is even worse: it only provides an O(1) error on regular grids and a O(∆ -1 ) error on distorted meshes (the error increases as the mesh is refined).

Alternatively, the Fig. 7b shows that, when using the new deconvoluted operators, both formulations lead to an O(∆ 2 ) convergence on regular meshes and an O(∆) convergence on distorted grids.

-16

log 2 ( E 2 ) 2D K nc M 2,F K 0 M 2,F K 0.2 K c M 2,F K 0 M 2,F K 0.2 4 6 8 10 4 0 -4 -8 -12
-16

log 2 (N x ) log 2 ( E 2 ) 4 6 8 10 log 2 (N x ) 3D K nc M 3 0 M 3 0.2 K c M 3 0 M 3 0.2 (

Application to flux computation for convection problems

Thanks to the deconvolution procedure exposed in section 2.3, it is possible to compute a high-order approximation of the solution everywhere in a given CV and more particularly at its boundaries. The objective of this section is to show that it is possible to use this polynomial expansion to compute the convection fluxes on the CV boundaries. This is simply achieved by performing the exact integration of this polynomial over the boundary faces. This can be done either by a Gauss quadrature or by using the surface moments on boundary faces. Both methods give the same result: the later method is less common but more straightforward here as it involves moments, just like the deconvolution technique.

Starting from Eq. 5, we are looking for an accurate evalutation of the flux ∂Ω I φu • dS. This integral is performed over the full CV boundary ∂Ω I which can always be split into a set of smaller simplicial pieces, called facelets hereafter:

∂Ω I = f ∂ f (56)
Thanks to the linearity of the integration, it is thus sufficient to be able to compute this integral on a generic simplicial facelet and then to sum up all the contributions to have to the full flux across the CV.

We will restrict the discussion to the evaluation of the convective flux with a constant velocity to avoid any problems arising from the discretization of the velocity field. The case of variable velocity can be addressed by two different level of approximations:

• the easiest and most straightforward solution is to use the value of the velocity at the center of each facelet.

If a linear interpolation is used, this method will provide a second-order approximation of the velocity.

• another option is to build a polynomial expansion of the velocity, just as for the scalar, and to compute the flux from the product of these two polynomials. Only the terms up to order 3 are needed to provide a third-order approximation of the flux.

First, we express that the velocity is constant over the facelet f and that this facelet has a constant normal vector n f everywhere (because it is a simplex) so that the contribution of the facelet f to the global residual at node I writes:

F I,f = ∂ f φu • dS = u • n f ∂ f φ dS (57)
Then we can replace φ by its third-order approximation (Eq. 29):

F I,f = u • n f ∂ f φ O3 I + G O2 I • ∆ I (x) + 1 2 H O1 I : (∆ I (x) ⊗ ∆ I (x)) dS + O ∆ 2+d (58)
We can now introduce the facelet surface (both scalar and vector):

S f = ∂ f dS and S f = S f n f (59)
and the various facelet moments with respect to I:

1 f = 1 S f ∂ f 1 dS = 1 (60) ∆ I f = 1 S f ∂ f ∆ I (x) dS (61)
∆ I ⊗ ∆ I f = 1 S f ∂ f ∆ I (x) ⊗ ∆ I (x) dS (62)
Similarly to the control-volume moments, all these quantities can be computed exactly as a function of the coordinates of the simplex vertices. Now Eq. 58 becomes:

F I,f = u • S f φ O3 I + G O2 I • ∆ I f + 1 2 H O1 I : ∆ I ⊗ ∆ I f + O ∆ 2+d (63)
It can be seen that the integrated flux is of order 2 + d, as the integration over the CV boundary introduces an additional d -1 order to the flux evaluation order. The residual will then be computed by dividing this integrated flux by the measure V I of the CV and the residual will thus be of second order (for a third-order deconvolution).

However, in the context of unsteady simulations, this residual is multiplied by the timestep to obtain the solution at the next discrete time. The error on the solution at the next iteration is thus O ∆ 2 × ∆t . If this timestep is subject to a stability constraint such as the CFL condition, it is proportional to the grid size and the error on the solution is finally O ∆ 3 .

Three different flux evaluations will be examined to assess the performance of the proposed method.

The first one, named S (for Simple), approximates the convective flux on facelet f from node I with the formula:

F S I,f = u • S f φ Ω I (64)
where the usual non-deconvoluted value of the transported variable is used in the integration.

The second one uses the deconvolution limited to order 2 and is named D 2 :

F D2 I,f = u • S f φ O2 I + G O1 I • ∆ I f (65)
Finally, the method D 3 uses the third-order deconvolution presented above:

F D3 I,f = u • S f φ O3 I + G O2 I • ∆ I f + 1 2 H O1 I : ∆ I ⊗ ∆ I f ( 66
)
The purpose of the present work is not to design a specific numerical scheme but to show that it is possible to build a high-order approximation of the flux at the control volume interface. The convective residuals associated to the aforementioned fluxes will be built in very simple ways, by considering either fully upwind (Eq. 67) or fully centered (Eq. 68) fluxes at the facelets:

F u f =    F I,f if u • IJ > 0 0 if u • IJ = 0 F J,f if u • IJ < 0 (67) F c f = 1 2 F I,f + F J,f ) (68)
where I and J are the two vertices that share the facelet f .

The 6 methods introduced above (S c , S u , D c 2 , D u 2 , D c 3 and D u 3 ) have been applied to the convection of the bump function described in section 3.1. The results presented below show the discrepancies between the numerical solution and the analytical solution after the signal has been convected for one period. Figures 8, 9 and 10 show the L 2 -norm of the errors for the same 1D, 2D and 3D domains than those described in section 3.1.

All three figures share the same trend except for one specific feature that will be discussed in detail below. It can be first noticed that the S schemes (first lines of the figures) exhibit their usual features: when used in conjunction with a centered formulation (S c ), it is second-order on regular meshes with a loss of accuracy on distorted meshes.

On the other hand, the upwind formulation (S u ) always produces a first-order convergence, whatever the mesh quality. Moreover, as it will be illustrated below, it introduces a large amount of numerical dissipation.

The D 2 scheme (second lines of the figures) shows a second-order convergence in all situations, whatever the mesh quality or the amount of upwinding. This is the expected behaviour: the method was demonstrated to be second-order in all situations. Moreover, it appears that just like in the S schemes, the D 2 scheme is stable even in the centered case and does not require any flux limiting technique.

Finally, the D 3 scheme (third lines of the figures) shows a specific behavior that was not observed with the S and D 2 schemes. When used with a centered approach (D c

3 ), it is unstable in 2D and in 3D (on distorted grids only): this is illustrated by the error which increases when the mesh is refined, on the lower left graph of Fig. 10. However, it can be observed that these instabilities totally disappear when an upwind formulation (D u

3 ) is used. This last choice still provides a third-order convergence on all meshes for 1D, 2D and 3D domains. It can also be observed that before becoming unstable, the D c 3 approach offers a fourth-order convergence. As shown recently by [START_REF] Guermond | Invariant domain preserving discretization-independent schemes and convex limiting for hyperbolic systems[END_REF], once a high-order evaluation of the flux is available, it is always possible to apply a convex-limiting technique to mix low-order and high-order approximations to guarantee both the precision and stability of the numerical scheme: such a technique could be used to mix F D c 3 and F D u 3 to have fourth-order convergence where the mesh is regular and third-order on distorted meshes.

To get a better idea of the improvement provided by the new formulation, the solution after one turn-over time are represented in Fig 11 (for the 1D case only). The left side of the figure is obtained with the centered formulation while the right side is obtained with the upwind scheme. The upper part corresponds to a regular mesh and the lower part represents the solution on a distorted mesh (this can be seen directly on the position of the nodes). Both meshes are composed of 32 cells.

The first observation is that the S schemes behaves poorly in all situations except for the S c formulation on regular grids. On a distorted mesh, large spurious oscillations are observed near the gradients. On the other hand, when using the S u formulation, the signal is strongly dissipated: the peak value of the bump has been divided by a factor of 2.

The second major observation is that both D 2 and D 3 schemes are very insensitive to either the mesh quality or the amount of upwinding. In all cases, the transported signal is mush closer to the exact solution, more particularly for the D 3 scheme. All these observations are confirmed by the values of the L 2 -norm of the errors associated to each signal, which are reported in the legend box of each figure.

As already mentioned, the deconvolution and flux reconstruction procedure presented here has a computational overhead compared to simpler methods such as the S scheme. However, this additional CPU remains low, and it is

M 1 0.2 M 1 0.3 M 1 0.4 -4 -8 -12 -16 -20 -24
-28 log ( E )

D 2 scheme M 1 0 M 1 0.1 M 1 0.2 M 1 0.3 M 1 0.4
log 2 ( E 2 ) S scheme M 2,F K 0 M 2,F K 0.1 M 2,F K 0.2 M 2,iso 0 M 2,iso 0.35 M 2 sof t -4 -8 -12 -16 -20 -24
-28

log 2 ( E 2 ) D 2 scheme M 2,F K 0 M 2,F K 0.1 M 2,F K 0.2 M 2,iso 0 M 2,iso 0.35 M 2 sof t 4 6 8 10 -4 -8 -12 -16 -20
-24 rapidly counterbalanced by the better accuracy of the method. This is illustrated by Fig. 12 which shows the error evolution with respect to the CPU time needed to obtain the solution for the 3D convection problem. Each point of the figure has been obtained with an M 3 0.2 grid but with different resolutions: the leftmost points corresponds to coarse grids while the rightmost points have been obtained on fine meshes. It can be observed that the D u 3 scheme can achieve a given accuracy much faster than the low order S c scheme.

-28 log 2 (N x ) log 2 ( E 2 ) 4 6 8 10 log 2 (N x ) D 3 scheme M 2,F K 0 M 2,F K 0.1 M 2,F K 0.2 M 2,iso 0 M 2,iso 0.35 M 2 sof t (

Conclusion

Increasing the order of convergence of Finite Volume schemes while keeping an efficient and compact formulation is highly challenging. The development of such schemes requires a rigorous framework which enables to control all the discretization errors. One of these errors comes from the fact that FV methods solve for the averaged quantities, which need to be deconvoluted to build local higher-order polynomial reconstructions of these quantities.

Several methods were proposed in the past to deal with this issue, such as the k-exact formulation (standard and successive correction formulations), the MOOD procedure or the Compact Least-Squares Finite-Volume method. These procedures were mainly studied in the context of the cell-centered FV formulation and to the authors' knowledge, less attention was paid to the node-centered FV approach.

In this paper, a deconvolution method is presented for the node-centered FV approach. It was shown that the node-centered formulation gives a general framework to attain a third-order reconstruction, using only the direct neighbors of a control volume in an unstructured grid composed of simplices. The reconstruction needs to solve a linear system Ax = b of size 10x10 (in 3D) or 6x6 (in 2D) or 3x3 (in 1D) on each node of the mesh, where A depends on local metrics only and thus can be computed and inverted in a pre-processing step. As a consequence, the CPU overhead associated to the proposed deconvolution is very limited. Near a boundary, the number of mesh vertices linked to a boundary node is not large enough to attain a third-order of accuracy but a second-order reconstruction can still be obtained. If needed, third-order can also be achieved at boundaries by using the polynomial expansion from neighbouring nodes.

Several applications from different physical domains were proposed and analyzed, namely high-order interpolation, levelset curvature computation and flux computation for hyperbolic transport equations. It is shown that the proposed technique is efficient and gives results of better quality than the standard scheme. Focusing on the Navier-Stokes equations, another analysis should be performed in the near future to treat the diffusion scheme as well as the case of variable transport coefficients.

A Grid moments computation

The technical details to compute the moments of a generic simplex are briefly explained hereafter. Let us consider a generic d-dimensional simplex, whose vertices coordinates are noted x 1 , . . . , x d+1 . We can first compute its volume with the classical formula:

V = Ω dV = 1 d! |X| ( 69
)
where X is a d × d matrix formed by the coordinates of d edges of the simplex that all share a common vertex, for example:

X =   x 2 -x 1 x 3 -x 1 x 4 -x 1 y 2 -y 1 y 3 -y 1 y 4 -y 1 z 2 -z 1 z 3 -z 1 z 4 -z 1   for d = 3 (70)
By adapting the formulas given by [START_REF] Rathod | Integration of polynomials over an arbitrary tetrahedron in euclidean three-dimensional space[END_REF] for a tetrahedron, it is then possible to derive the first two moments of any simplex, with respect to the origin.

M 1 O = 1 V Ω x dV (71) = 1 d + 1 d+1 i=1 x i (72) M 2 O = 1 V Ω x ⊗ x dV (73) = 1 (d + 1)(d + 2)     d+1 i=1 d+1 j=1 x i ⊗ x j   + d+1 i=1 x i ⊗ x i   (74)
Deriving similar expressions for higher order moments, although tedious, is not a problem and can be conducted for any simplex.

The next ingredient that is needed to compute the moments of a simplex with respect to any point I is the translation formulas, which are given here only for orders 1 and 2 but which can be derived to any order without any difficulty with the binomial formula in a non-commutative ring:

M 1 I = M 1 O -x I (75)
M 2 I = M 2 O -x I ⊗ M 1 O -M 1 O ⊗ x I + x I ⊗ x I (76)
Finally, thanks to the linearity of the integration operator, the full control-volume moments can be computed by first performing a simplicial decomposition of the CV and then by summing up the moments of the elementary simplices which compose the control-volume (as long as they are computed with respect to the same position, I in the present case).

B Mesh Quality

In the present work, numerical experiments are performed to assess the quality of various operators, i.e. their order of convergence, both on regular and distorted simplicial meshes. The specific procedure that was used to build these meshes is detailed below.

In 1D, the procedure is straightforward: the 1D domain of size L is divided into N segments of identical size to produce a regular mesh whose cells have an equal size of ∆ = L N . The 1D distorted meshes are then constructed by applying a perturbation to each node of this regular mesh. For the node I, the amplitude of this perturbation is equal to r I ∆ where r I is a random variable with a uniform distribution on [-b/2; +b/2], where b is a perturbation parameter which is fixed for a given mesh. Those meshes will be noted M 1 b . In 2D, a regular triangular mesh is built by performing a Friedrichs and Keller triangulation of a Cartesian mesh (see Fig. 13-a). Following the same procedure as in 1D, this grid is then perturbed by randomly shifting its vertices inside a square box of size b∆ surrounding each node, as depicted in Fig. 13-b. Those meshes will be noted M 2,F K b .

Just by stretching and shifting the M 2,F K 0 mesh described above (but without modifying the connectivity), it is also possible to build a mesh only composed of equilateral triangles (Fig. 13-c). In this situation, the control volume is a regular hexagon: this is the most isotropic CV that can be obtained with this connectivity. This mesh can also be perturbed by the technique described above and those meshes will be noted M 2,iso b . An example of M 2,iso 0.35 mesh is shown in Fig. 14 to illustrate the type of mesh which is concretely used.

Based on the same idea, the 3D tetrahedral meshes that are used in the present paper are built in two steps.

First, a Cartesian hexahedral grid is built and its vertices are shifted randomly into a cube of size b∆. Then, each hexahedron is divided into six tetrahedra following the method proposed by Dompierre et al. [START_REF] Dompierre | How to subdivide pyramids, prisms and hexahedra into tetrahedra[END_REF] (see Fig. 13-d).

Those meshes will be noted M 3 b . It must be noticed that the perturbation parameter b can not exceed a certain value (which depends on the initial mesh) to ensure that all the cells have a positive volume.

C Analytical evaluation of the convolution matrix

In some simple cases, the convolution / deconvolution procedure can be performed analytically. This is particularly interesting for a distorted 1D mesh and for the mesh referred to as M 2,iso 0 . The aim of this appendix is to present briefly these two results.

With respect to the notations introduced in Fig. 15, the convolution matrix on a distorted 1D mesh reads:

C O3 I =      1 h + -h - 4 1 24 (h + ) 2 -h + h -+ (h -) 2 0 G I ∆ I Ω 1 2 G I ∆ I ⊗ ∆ I Ω 0 H I ∆ I Ω 1 2 H I ∆ I ⊗ ∆ I Ω      , (77)
with

G I ∆ I Ω = 1 4(h + + h -) h ++ + 3h + + 3h -+ h --, (78)
G I ∆ I ⊗ ∆ I Ω = 1 12(h + + h -) (h ++ ) 2 + 5h + h ++ + 7(h + ) 2 -(h --) 2 -5h --h --7(h -) 2 , ( 79
)
H I ∆ I Ω = 1 2h -h + (h + + h -) h -h ++ + h --h + h --+ h + , (80)
H I ∆ I ⊗ ∆ I Ω = 1 6h -h + (h + + h -) h -(h ++ ) 2 + 5h -h + (h ++ + h --) + 7h -h + (h + + h -) -(h -) 3 + h + (h --) 2 -(h + ) 3 . (81)
Its determinant,

det C O3 I = (h -+ 2h + + h ++ ) (h --+ 2h -+ h + ) (h --+ 2h -+ 2h + + h ++ ) 48h -h + (h -+ h + ) , (82)
x i-2

x i-1 x i x i+1 x i+2 x i-1 2 x i+ 1 2 x i-1 x i x i+1 h -- h - h + h ++ ∆ i-1
Ωi-1 On 1D uniform mesh of size h, the convolution matrix Eq. 77 reduces to the well-known link between the nodal value and the cell averaged value:

∆ i Ωi ∆ i+1 Ωi+1
C O3 I =   1 0 h 2 24 0 1 0 0 0 1   . ( 83
)
For a 2D regular grid composed of equilateral triangles (Fig. 16, the same case that in [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF]), the convolution matrix that arises in our framework is: and the resolution of this system leads to:

C O3 I =              1 0 0 5h 2
             φ I ∂ x φ I ∂ y φ I 1 2 ∂ xx φ I 1 2 ∂ yy φ I ∂ xy φ I              =                 
246φ I -5(φ J 1 +φ J 2 +φ J 3 +φ J 4 +φ J 5 +φ J 6 ) 216 2(φ J 1 -φ J 4 )+φ J 2 +φ J 6 -φ J 3 -φ J 5 6h √ 3(φ J 2 +φ J 3 -φ J 5 -φ J 6 )

6h φ J 1 -2φ I +φ J 4 2h 2
2(φ J 2 +φ J 3 +φ J 5 +φ J 6 )-6φ I -φ J 1 -φ J 4 ) 6h 2 √ 3(φ J 2 +φ J 5 -φ J 3 -φ J 6 )

3h 2                  . ( 85
)
The components of the gradient are identical to the equation 13 of [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF] which is obtained with the Least-Square method in the linear reconstruction case.

Figure 1 :

 1 Figure 1: Notations for a 2D control-volume based on the barycentric dual mesh. Here, N I = {J, K, L}. Considering the boundary face ∂Ω IJ between nodes I and J (red line), the full face vector S IJ is the sum of the two facelets surface vectors s 1IJ and s 2 IJ . In 3D, the situation is similar but the number of facelets based on an edge is not fixed (generally 5 to 8 facelets by edge). The vector n IJ represents the unit vector direction on edge IJ.

I

  is a third-order approximation of φ I , G O2 I a second-order approximation of ∇φ I and H O1 I a first-order approximation of ∇∇φ I , respectively.

Figure 2 :

 2 Figure 2: Probability Density Function of the condition number of the convolution matrix for various meshes

IFigure 3 :

 3 Figure 3: Errors on Hessian for different mesh quality, without deconvolution (a) and with deconvolution (b). Meshes quality notations are detailed in Appendix B. Black dotted lines correspond to convergence order 1, 2, 3 and 4 respectively.

IFigure 4 :

 4 Figure 4: Errors on gradient for different mesh quality, without deconvolution (a) and with deconvolution (b). Meshes quality notations are detailed in Appendix B. Black dotted lines correspond to convergence order 1, 2, 3 and 4 respectively.

IFigure 5 :

 5 Figure 5: Errors on scalar for different mesh quality, without deconvolution (a) and with deconvolution (b). Meshes quality notations are detailed in Appendix B. Black dotted lines correspond to convergence order 1, 2, 3 and 4 respectively.

Figure 6 :

 6 Figure 6: Relative L 2 -norm of error on interpolation by using Délèze methodology presented in section 3.2. Black dotted lines correspond to convergence order 1, 2 and 3 respectively.

  a) Low Order (b) High Order

Figure 7 :

 7 Figure 7: Errors on curvature computation using both formulations presented in section 3.3. Black dotted lines correspond to convergence order 1, 2.

1 Figure 8 :

 18 Figure 8: Errors after one turn-over in 1D using the S, D 2 and D 3 schemes. Dotted lines correspond to convergence order 1, 2, 3 and 4 respectively.

2 Figure 9 :

 29 Figure 9: Errors after one turn-over in 2D using the S, D 2 and D 3 schemes. Dotted lines correspond to convergence order 1, 2, 3 and 4 respectively.

3 Figure 10 :

 310 Figure 10: Errors after one turn-over in 3D using the S, D 2 and D 3 schemes. Dotted lines correspond to convergence order 1, 2, 3 and 4 respectively.

Figure 12 :

 12 Figure 12: Error after one turnover plotted as a function of the CPU time needed to obtain the solution. Each point corresponds to a different mesh resolution. The data were collected on a 3D domain with distorted grids (M 3 0.2 ).

∆

Figure 13 :

 13 Figure 13: Examples of manufactured meshes used for the convergence order tests. (a), (b) and (c) correspond to 2D situations with the same coloring conventions as above while (d) represents the procedure used to partition an hexahedron into six tetrahedrons.

Figure 14 :

 14 Figure 14: Example of M 2,iso 0.35 manufactured mesh with N x = 16.

Figure 15 :

 15 Figure 15: Notations for a distorted 1D barycentric dual mesh.

Acknowledgement

Dr. Guillaume Puigt is partially supported by LMA2S (Laboratoire de Mathématiques Appliquées à l'Aéronautique et au Spatial), the Applied Mathematics Lab of ONERA. We thank PSPC Innov'Hydro for funding M. Bernard, GENCI (Grant 2A00611) and CRIANN for providing the computing ressources for the various test cases.


Bau Plinio
email: plinio.bau@irt-saintexupery.com

Marc Cousineau

Bernardo Cogo

Frédéric Richardeau

Vinnac Sébastien

Didier Flumian

Nicolas C Rouger
email: nicolas.rouger@laplace.univ-tlse.fr

Plinio Bau

Bernardo Cougo

Frederic Richardeau

Sebastien Vinnac

Sub-nanosecond delay CMOS Active Gate Driver for Closed-Loop dv/dt Control of GaN Transistors

Keywords: Active gate driver, GaN, switching analysis, dv/dt, EMI, power electronics, ASIC for power ic

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Different solutions have been previously demonstrated to control separately the dv/dt and di/dt sequences with silicon, SiC and GaN FETs in the view of improving the loss versus EMI tradeoff. Different strategies to improve this tradeoff are: a variable impedance output stage, an open-loop control with previously programmed, an adjusted impedance sequence [START_REF] Takamiya | Power electronics 2.0: IoT-connected and Al-controlled power electronics operating optimally for each user[END_REF][START_REF] Rose | A GaN HEMT driver IC with programmable slew rate and monolithic negative gate-drive supply and digital current-mode control[END_REF][START_REF] Yu | A smart gate driver IC for GaN power transistors[END_REF][START_REF] Dymond | A 6.7-GHz Active Gate Driver for GaN FETs to Combat Overshoot, Ringing, and EMI[END_REF][START_REF] Park | Flexible dv/dt and di/dt control method for insulated gate power switches[END_REF] or a closedloop control to reduce the gate current during the dv/dt sequence [START_REF] Sun | Active dv/dt control of 600V GaN transistors[END_REF][START_REF] Grezaud | A Gate Driver With Integrated Deadtime Controller[END_REF]. The open-loop solutions rely on pre-optimization and fine tuning of the variable impedance switching sequence, which can be sensitive to parameter dispersion [START_REF] Takamiya | Power electronics 2.0: IoT-connected and Al-controlled power electronics operating optimally for each user[END_REF][START_REF] Rose | A GaN HEMT driver IC with programmable slew rate and monolithic negative gate-drive supply and digital current-mode control[END_REF][START_REF] Yu | A smart gate driver IC for GaN power transistors[END_REF][START_REF] Dymond | A 6.7-GHz Active Gate Driver for GaN FETs to Combat Overshoot, Ringing, and EMI[END_REF][START_REF] Park | Flexible dv/dt and di/dt control method for insulated gate power switches[END_REF]. Closed-loop controls are either based on discrete components [START_REF] Sun | Active dv/dt control of 600V GaN transistors[END_REF][START_REF] Grezaud | A Gate Driver With Integrated Deadtime Controller[END_REF] or integrated solutions [START_REF] Yu | A smart gate driver IC for GaN power transistors[END_REF]. The discrete solutions typically have a large response time, above few nanoseconds or tens of nanoseconds. Our approach consists to get a full CMOS integration for all the required functions. The size of the high voltage capacitor required to sense the dvDS/dt has been reduced to a few pico-Farads, which can be integrated on-chip in our work. This method can be used both for turn-on and turn-off transients. However, the first developments concern only the active control of the turn-on transient. The turn-off sequence is typically controlled by the output load current in the case of high switching speed (low switching losses), and the Zero Voltage Switching (ZVS) condition is more critical than the turn-off switching speed, in the case of synchronous buck converter [START_REF] Bau | A CMOS gate driver with ultra-fast dV/dt embedded control dedicated to optimum EMI and turnon losses management for GaN power transistors[END_REF]. However, it should be mentioned that the proposed active gate driver can be easily modified to actively control the turn-off of a fast switching GaN transistors using the exact same principle as the one detailed hereafter.

II. PRINCIPLE OF THE ACTIVE GATE DRIVER

Fig. 1 shows the transistor level schematics of the proposed CMOS AGD. Additionally, to the required output buffer (M5-M6), a high voltage sense capacitor CS and a two stage current mirror M1-M2 and M3-M4 are integrated. During the dv/dt sequence, the dv/dt sense current iCs is amplified and participates to reduce the gate current ig during the Miller plateau (Fig. 2).

Consequently, the dv/dt is reduced, while keeping a fast di/dt sequence unchanged. The trade-off between dv/dt and the turn-on energy loss EON is then optimized, comparatively to a change of the gate resistor RG, as already demonstrated [START_REF] Sun | Active dv/dt control of 600V GaN transistors[END_REF][START_REF] Grezaud | A Gate Driver With Integrated Deadtime Controller[END_REF]. https://doi.org/10.1109/ISPSD.2019.8757693 As presented earlier, the transistors M1-M6 are 5V transistors which are particularly suited to drive efficiently EPC GaN transistors, and can also drive GaN systems transistors, albeit with non-optimal driving voltages. The transistors M5 and M6 are designed with a +/-3A source/sink current rating, with integrated pre-amplification and short-circuit protection. DP is a protection diode to prevent the current mirror input voltage vinM from achieving values outside the safe operating region of the transistor M1. ibias is a current source improving the response time of the feedback loop, which will be further discussed in section IV. Fig. 2 shows the qualitative current and voltage waveforms during a turn-on of the GaN power device. This time interval can be divided into the di/dt (between t1 and t2) and the dv/dt (between t2 and t3) sequences. The time duration tD-t2 is the feedback loop delay, further analyzed in section IV. The main idea of this method is to emulate a virtual capacitor added between the gate and drain nodes of the power device only during the dv/dt phase. The value of such an equivalent capacitor can be controlled by the gain G of the feedback loop (gain of the cascaded current mirrors M1-M2 and M3-M4). One should note that this gain G results from the first stage PMOS transistors and the second stage NMOS transistors G = GP × GN. The emulated capacitance effect occurs when a feedback current iFB is sinking from the gate of the power device during the dv/dt phase. This feedback current is generated by the current mirror gain G and the CS sense capacitor value. This added emulated capacitance equal to GCS can be easily computed using a simplified model for the power transistor where the dv/dt of the power device is expressed by (1):

d𝑣 DS d𝑡 = - 𝑣 DRV -𝑣 M 𝑅 G 𝐶 GD + 𝑖 FB 𝐶 GD = - 𝑣 DRV -𝑣 M 𝑅 G (𝐶 GD +𝐺𝐶 S ) (1)
where vDRV is the power supply voltage of the driver and vM is the Miller plateau voltage, shown in Fig. 2.

This principle has been previously and successfully applied to the control of Silicon IGBT [START_REF] Sun | Active dv/dt control of 600V GaN transistors[END_REF] and high voltage GaN FET [START_REF] Grezaud | A Gate Driver With Integrated Deadtime Controller[END_REF] for which the turn-on time is significantly high. However, this closed-loop approach is particularly difficult to achieve when GaN

B. Implementation in 48V and 400V DC bus voltages

Two different boards have been developed to demonstrate the active control of GaN transistors in 48V and 400V applications.

In both experiments, only the low side (LS) GaN transistor is driven by our CMOS AGD, while the high side (HS) GaN transistor has its gate and source shorted. The 4-layer PCB (48V) and 6-layer PCB (400V) have been designed to minimize the gate and power loops, hence only voltage probes are used. In this first implementation, an external CS ceramic capacitor is used, 2pF/250V for 48V (J0603D2R0BXPAJ) and 1pF/500V for 400V (MC0805N1R0C501CT).

C. Experimental Results

In the case of 48V bus voltage, the dv/dt is attenuated from 15V/ns to 6.3V/ns with the CMOS AGD (Fig. 5). Simulation results shows the switching losses are increased from 1.93µJ to 2.34µJ during turn-on. Cadence ™ simulation results fitting experimental data shows the switching losses are reduced by 15.7% comparing to the increase of the gate resistance RG, while keeping the same attenuation in dv/dt. For the 400V application, the dv/dt is attenuated from 45V/ns to 6.6V/ns, (Fig. 6). Both Author version of the accepted article published in the IEEE proceedings ISPSD2019https://doi.org/10.1109/ISPSD.2019.8757693

IV. RESPONSE TIME OF THE FEEDBACK LOOP

A. Reducing the delay time by pre-biasing the current mirrors

One way to reduce the delay of the feedback loop (tD-t2) is to pre-bias the current mirrors M1-M2 and M3-M4 by precharging their capacitors with a low biasing current ibias. Doing so, the pre-biased transistors, biased in their linear region, are able to provide a very fast response for any external transient. Table I shows the results obtained with Cadence ™ -Spectre ™ transient simulations to determine the value of the delay defined in Fig. 2 for different biasing currents ibias. Different total gain values (5, 10, 20 and 50) for the current mirrors are simulated and implemented within the prototype. One has to note that all transistors M1-M4 are designed differently for each case of total gain G value. If different gain modifies the dv/dt in closedloop and consequently the input current iCs, it also implies different parasitic capacitance values which impact the overall response time. In these simulations, for comparison purpose, a constant dv/dt equal to 15V/ns with a sense capacitor of 2pF is considered. The design must be done carefully, splitting the total gain G into two different gain GP and GN to minimizing the parasitic capacitors of the transistors, leading to an optimized size reducing the several time constants of the feedback loop. * shown in Fig. 7 The time intervals tD-t2 and tE-tD are determined on Fig. 7 showing a simulation of the response time with the current mirror G = 20 pre-biased with a 2mA ibias current. The input current goes from 0 to 26mA that correspond to a dv/dt equal to either 26V/ns with CS = 1pF or 13V/ns with CS = 2pF.

After tens of picoseconds the feedback loop provides a current at the output (iFB) and therefore is already acting and affecting the dv/dt.

The delay decreases as the bias current ibias increases. This improvement is at the expense of an additional power consumption. Indeed, using a pre-bias current, extra losses in transistors M1 to M4 are generated, which can reach high levels due to the high value of the gain GP. It should be noted that, this biasing current can be provided only during the switching.

Then, a trade-off between extra losses and bandwidth improvement has to be considered.

Author version of the accepted article published in the IEEE proceedings ISPSD2019https://doi.org/10.1109/ISPSD.2019.8757693

B. Bandwidth analysis

If a non-null current ibias is used to pre-bias the two current mirrors, a small signal analysis can be performed to determine the bandwidth of the feedback loop. Fig. 8 shows the equivalent circuit required for the study in the frequency domain. The input signal is the current coming from the capacitor CS during a turn-on event assuming that a constant dv/dt occurs during the transient. The output of the second stage is connected to the resistance RG and the gate of the power device behaving like a voltage-dependent capacitor CISS during the Miller's plateau. Small signal models of transistors M1 to M4 are used in Fig. 8:

where CDP is the capacitive contribution of the reverse biased protection diode DP.

A Cadence™ AC simulation (Fig. 9) demonstrates a bandwidth close to 550MHz (-3dB gain) for the two-step current mirror with a current gain G = 20 and considering an input current made of a biasing current ibias = 2mA increased by a 12mA current generated by a constant dv/dt equal to 6v/ns in a sensing capacitor CS = 2pF.

Fig. 1 .

 1 Fig. 1. CMOS AGD topology in a transistor level schematic.

Fig. 2 .

 2 Fig. 2. Qualitative waveforms for gate voltage and currents with the AGD during turn-on (closed-loop). (a) Current and voltage at the power gate and (b) currents being supplied and sinked simultaneously by the AGD.

Fig. 3

 3 Fig.3shows a microscope photograph of the fabricated AGD in AMS 180nm CMOS technology. This first prototype has been designed as a test chip with several current mirror designs and embedded high voltage capacitors. Even though quite simple, the design of the current mirrors is subject to important compromises such as the feedback loop gain splitting between the CS value, the first and second current mirror stage gain values, the large maximum output current in M4, the dynamic bias current ibias and the M2-M3 biasing during the active driving. Key constraints are current consumptions, total silicon area (embedded capacitor + current mirrors) and bandwidth. These design considerations will be illustrated in section IV. After a quick characterization in probe station [9], the AGDs are packaged in QFN24 6mm×6mm. The packaged AGDs are then implemented to drive commercially available low voltage and high voltage GaN FETs.

Fig. 3 .

 3 Fig. 3. Optical microscope photograph of the CMOS AGD test chip built in 0.18µm technology.

Fig. 4 (Fig. 4 .

 44 Fig.4(a) shows the system implementation in a 48V-36A power commutation cell with EPC2001 eGaN ™ FET. Fig.4 (b)shows a similar implementation in a 400V-30A application, with GaN systems e-mode GS66508T transistors. Additional components are limited to decoupling capacitors, protection and configuration buffers and gate resistors.

  cases use RG = 4.4Ω. This choice of RG has been made for the EPC2001C because characterization tests in open-loop with lower values showed no more significant increase in peak dv/dt with lower RG values. For 400V tests, the same value is used for comparative purposes.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Experimental results with the CMOS AGD in the 48V DC bus voltage. With the closed-loop, the dv/dt is attenuated from 15V/ns to 6.3V/ns.

Fig. 7 .

 7 Fig. 7. Simulation of the transient response time with G = 20 (vDRV = 4V, ibias = 2mA, iCs = 30mA).

Fig. 8 . 1 =

 81 Fig. 8. Small signal analysis of the AGD current mirrors during dv/dt phase. The current gain transfer function is expressed by (2).

  Author version of the accepted article published in the IEEE proceedings ISPSD2019https://doi.org/10.1109/ISPSD.2019.8757693FET are switching within a few nanoseconds. This work aims to demonstrate experimentally the feasibility and the limitations of such an AGD with both low voltage and high voltage GaN FETS, thanks to a full CMOS integration.

TABLE I .

 I COMPARISON OF PERFORMANCE FOR THE SAME METHOD OF DV/DT CLOSED-LOOP CONTROL

		Published by [6]	Published by [7]	This work
	Reduction in | dV/dt |	1.8 to 0.5V/ns	27 to 8V/ns	45 to 6.6V/ns
	vHVdc	600V	300V	400V and 48V
	vDS fall
	time during	333ns	11ns	20ns and 4ns
	turn-on
	Embedded in an ASIC	No	No	Yes
		III. IMPLEMENTATION
	A. CMOS AGD

Author version of the accepted article published in the IEEE proceedings ISPSD2019https://doi.org/10.1109/ISPSD.2019.8757693 With this large bandwidth, the step response of the cascaded mirrors provides a settle time equal to (tE-tD)50% = 220ps that is close to the result provide on table II.

V. EMBEDDED CAPACITOR FOR