A Bazile
email: alban.bazile@mines-paristech.fr

E Hachem

J C Larroya-Huguet

Y Mesri

Variational Multiscale error estimator for anisotropic adaptive fluid mechanic simulations: application to convection-diffusion problems

Keywords: CFD, VMS, Error Estimator, Mesh Adaptation, Convection-Diffusion

In this work, we present a new a posteriori error estimator based on the Variational Multiscale method for anisotropic adaptive fluid mechanics problems. The general idea is to combine the large scale error based on the solved part of the solution with the sub-mesh scale error based on the unresolved part of the solution. We compute the latter with two different methods: one using the stabilizing parameters and the other using bubble functions. We propose two different metric tensors H iso and H new aniso . They are both defined by the recovered Hessian matrix of the solution and rely on the sub-grid scale error estimator. Thus, we develop a new anisotropic local error indicator and we test it for mesh adaptation on convection-dominated benchmarks in 2D and 3D. The results show that the proposed error indicator leads to enhanced and accurate solutions while using a drastically reduced number of elements.

Introduction

The use of Computational Fluid Dynamics (CFD) for industrial applications has been in constant increase for the last decades. Researchers are continuously developing new techniques to reach higher level of precision. Nevertheless, to comply with industrial expectations, a trade-off has to be found between high precision results and high computational costs [START_REF] Slotnick | CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences[END_REF]. Different strategies can be found in the literature. Most of them are related to high order elements (see [START_REF] Brogniez | A high-order discontinuous Galerkin method with Lagrange multipliers for advection-diffusion problems[END_REF][START_REF] Wang | High-order CFD methods: current status and perspective[END_REF][START_REF] Sevilla | A high-order stabilised ALE finite element formulation for the Euler equations on deformable domains[END_REF]), parallel computing (see [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF][START_REF] Meng | Scalable large-scale fluid-structure interaction solvers in the Uintah framework via hybrid task-based parallelism algorithms: Uintah hybrid task-based parallelism algorithms[END_REF][START_REF] Ibanez | Hybrid MPI-thread parallelization of adaptive mesh operations[END_REF][START_REF] Panitanarak | A parallel log barrier-based mesh warping algorithm for distributed memory machines[END_REF]) or, in particular, adaptive methods (see [START_REF] Mesri | Continuous mesh adaptation models for CFD[END_REF][START_REF] Jannoun | Fully space-time metric based anisotropic mesh adaptation for unsteady problems[END_REF][START_REF] Hassan | An anisotropic adaptive finite element algorithm for transonic viscous flows around a wing[END_REF][START_REF] Hachem | Unified adaptive Variational MultiScale method for two phase compressible-incompressible flows[END_REF]). Indeed, adaptive methods make it possible to improve the accuracy and the efficiency of numerical methods. In particular, anisotropic mesh adaptation has proved to be powerful in capturing dynamically the heterogeneities that can appear in numerous physical applications including those having boundary or inner layers [START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Aubry | Boundary layer mesh generation on arbitrary geometries[END_REF]. In these cases, gradients of the solution are highly directional and can be captured with a good accuracy using fewer additional elements. These mesh adaptation techniques are based on local modifications of an existing mesh. Usually, it consists in a local stretching of the elements which is defined by a metric field. This metric field is built from an error analysis on the mesh. In [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF][START_REF] Hachem | On the stabilized finite element method for steady convection-dominated problems with anisotropic mesh adaptation[END_REF][START_REF] Jannoun | Anisotropic meshing with time-stepping control for unsteady convection-dominated problems[END_REF] for example, the error analysis is done on the edges of the elements.

In fact, theories of anisotropic error estimation have been well developed, leading to some standardization of the adaptation process. Error estimation of the discretization error and in particular, of the interpolation error have been performed in a number of works such as in [START_REF] Almeida | Adaptive finite element computational fluid dynamics using an anisotropic error estimator[END_REF][START_REF] Mesri | On optimal simplicial 3d meshes for minimizing the Hessian-based errors[END_REF][START_REF] Loseille | Continuous Mesh Framework Part I: Well-Posed Continuous Interpolation Error[END_REF]. From the interpolation error analysis, several recent results [START_REF] Agouzal | Minimization of gradient errors of piecewise linear interpolation on simplicial meshes[END_REF][START_REF] Labbé | Verification of threedimensional anisotropic adaptive processes[END_REF][START_REF] Loseille | Continuous Mesh Framework Part I: Well-Posed Continuous Interpolation Error[END_REF][START_REF] Mesri | Automatic coarsening of three dimensional anisotropic unstructured meshes for multigrid applications[END_REF][START_REF] Kamenski | How a Nonconvergent Recovered Hessian Works in Mesh Adaptation[END_REF] have brought renewed focus on metric-based mesh adaptation where the underling metric is derived from a recovered Hessian. Indeed, the Hessian based metric mesh adaptation has several advantages from which we note: (i) the general computation framework, (ii) the relatively easy way of implementation and above all, (iii) the robustness.

Despite the practical construction of this kind of anisotropic error estimators, the information derived from them is only an indicator for the mesh adaptation mechanics. Indeed, the fine scale features related to some dynamic flow solutions are still difficult to capture without considering the PDE-dependent approximation error. The estimation of the latter involves the development of multiscale a posteriori error estimators that can be computationally prohibitive.

However, the Variational Multiscale method (VMS) that we use to stabilize our continuous finite element scheme provides, by construction, a cost free PDE-dependent a posteriori error estimator. The VMS approach introduced in [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF][START_REF] Hughes | A new finite element formulation for computational fluid dynamics: I. Symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics[END_REF][START_REF] Hughes | A new finite element formulation for computational fluid dynamics: VIII. The galerkin/leastsquares method for advective-diffusive equations[END_REF][START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF] consists in the splitting of the solution into a resolved part (i.e. coarse scales) and an unresolved part (i.e. subscales). The resolution of unresolved part gives a direct access to the sub-mesh scale information of the solution and allows us to compute an approximation error estimator without solving any additional equation. Recently, several works related to variational multiscale (or subscales) error estimators in the VMS framework have been published. Granzow et al. in [START_REF] Granzow | Output-based error estimation and mesh adaptation for variational multiscale methods[END_REF] developed an error representation for output quantities based on a dual enrichment technique. Also, in [START_REF] Baiges | Variational Multiscale error estimators for solid mechanics adaptive simulations: An Orthogonal Subgrid Scale approach[END_REF], Baiges et al. proposed a general error estimator for the finite element solution of solid mechanics problems. John et al. proposed in [START_REF] John | A robust SUPG norm a posteriori error estimator for stationary convection-diffusion equations[END_REF] a robust residual-based a posteriori estimator for the SUPG finite element method applied, in particular, to stationary convection-diffusion-reaction equations.

In this work, we compute the subscales error estimator using two different methods. The first one uses stabilizing parameters derived from the VMS analysis and local norms defined on the elements. It has been developed by Hauke et al. in [32,[START_REF] Hauke | Variational multiscale a-posteriori error estimation for multi-dimensional transport problems[END_REF] and provides an element-wise computation of the error . Referring to Irisarri et al., the second method uses a linear combination of bubble functions to establish a pointwise computation of the error. It has been developed for one-dimensional transport equation in [START_REF] Irisarri | Pointwise Error Estimation for the One-Dimensional Transport Equation Based on the Variational Multiscale Method[END_REF] and for 2D transport equation in [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF].

The first contribution of this paper is then to propose a new isotropic mesh adaptation technique based on the previous subscales error estimator. To do so, we define a new isotropic metric tensor H iso and we solve an optimization problem under the constraint of a fixed number of elements. Until now, the variational multiscale error estimators were only derived for isotropic mesh adaptation. In our work, we also propose to use the sub-mesh scale information for anisotropic mesh adaptation.

The second contribution is therefore to combine both the coarse scales interpolation error indicator and the subscales error estimator for anisotropic mesh adaptation. To the best of our knowledge, such strategy has never been tempted. To do so, we derive a new anisotropic metric tensor H new aniso that allows to take into account the anisotropic variations of the solution on the mesh and also relies upon the sub-grid information of the solution. This new anisotropic mesh adaptation technique takes into account the finite element VMS framework used for the numerical resolution of the convection-diffusion equation. In fact, the VMS approach allows to quantify the a posteriori subscales errors that can be used as weight to enrich the anisotropic error indicator of the conventional Hessian based approach. Unlike existing anisotropic a posteriori error estimators like [START_REF] Papastavrou | A posteriori error estimators for stationary convection-diffusion problems: a computational comparison[END_REF][START_REF] Formaggia | An anisotropic a-posteriori error estimate for a convection-diffusion problem[END_REF][START_REF] Picasso | Adaptive finite elements with large aspect ratio based on an anisotropic error estimator involving first order derivatives[END_REF] to mention a few, this hybrid approach takes advantages from both a priori and a posteriori methods. Indeed, we build an isotropic a posteriori error estimator targeting the subscales and we combine it with an anisotropic a priori interpolation based error indicator targeting the coarse scales.

The paper is structured as follows. In Section 2, we introduce the Variational Multi-Scale framework. From this mathematical background, we present, in Section 3, the two different computation methods of the subscales error estimator. Then, mesh adaptation is presented in Section 4. As a reminder, we introduce the Hessian based anisotropic mesh adaptation in Section 4.1. Then, we define the isotropic metric tensor H iso built from the subscales error estimator in Section 4.2. Finally, we present the new combination of the interpolation error indicator and the subscales error estimator in Section 4.3. This combination is tested on 2D and 3D convection-diffusion benchmarks in Section 5. In this last section, an error analysis of the exact error is made to highlight the two contributions of the paper.

Convection-diffusion equation in the Variational Multi-Scale framework

To begin this paper, let's introduce some notations. Referring to Irissarri et al. in [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF], we express the partial differential equation as:

   Lu = f in Ω u = g on Γ g Bu = h on Γ h (1)
where B is an operator acting on the boundary that defines the natural boundary condition; and L is a generic differential operator. g is the value of the Dirichlet boundary condition and h determines the Neumann boundary condition. In this work, we consider only Dirichlet boundary conditions. For the case of a convection-diffusion equation, we have:

Lu = -a∆u + v • ∇u ( 2
)
where a is the diffusion coefficient and v the velocity field.

We define the solution and test function spaces as standard Sobolev spaces:

S = {u ∈ H 1 (Ω) | u = g on Γ g } V = {w ∈ H 1 (Ω) | w = 0 on Γ g } (3)
The variational formulation of Eq. 1 is:

Find u ∈ S such that:

a(w, u) = (w, f ), ∀w ∈ V (4)
where a(., .) is a bilinear form, (., .) the L 2 (Ω) inner product.

The Galerkin finite element formulation

Applying the Finite Element Method (FEM), we mesh the domain into n el non-overlapping elements Ω e . We write Ω and Γ as:

Ω = n el e=1 Ω e Γ = n el e=1 Γ e \ Γ (5)
We define the Finite Element solution and weighting function spaces S h ⊂ S and V h ⊂ V such that:

S h = {u h ∈ H 1 (Ω) | u h|Ωe ∈ P k , u h|Γg = g ∀Ω e ∈ Ω} V h = {w h ∈ H 1 (Ω) | w h|Ωe ∈ P k , w h|Γg = 0 ∀Ω e ∈ Ω} (6)
with P k the space of polynomial of degree k. In this paper, we will only use linear elements, hence k = 1.

With these definitions, we apply the standard Galerkin method:

Find u h ∈ S h such that, a(w h , u h ) = (w h , f ), ∀w h ∈ V h (7)
It is well known that this formulation is unstable and leads to spurious oscillations when the convective term of the equation is dominant. For this reason, we stabilize the formulation using the Streamline Upwind Petrov-Galerkin (SUPG) numerical scheme.

Streamline Upwind Petrov-Galerkin (SUPG) method

The Streamline Upwind Petrov-Galerkin (SUPG) method consists in introducing stabilizing terms in the standard Galerkin formulation. To do so, for all terms in Eq. ( 7), we replace the weighting function w h by a new weighting function w h + τ e v • ∇w h . This modification of the formulation is usually interpreted as adding more weight to the node upstream, reducing the weight on the node downstream. It adds an artificial weighted diffusion along the streamline direction. This stabilization is done locally and on each element we have:

a(w h + τ e v • ∇w h , u h ) = (w h + τ e v • ∇w h , f ) ∀w h ∈ V h (8)
Thanks to the bi-linearity of a(., .) and of the L 2 (Ω) inner product, the formulation can be written as:

a(w h , u h ) + a(τ e v • ∇w h , u h ) -(τ e v • ∇w h , f ) = (w h , f ) ∀w h ∈ V h (9)
And so,

a(w h , u h ) + (τ e v • ∇w h , Lu h -f ) = (w h , f ) ∀w h ∈ V h (10)
Finally, summing on each element, we have a new bi-linear form a τ (., .) and we can write the following formulation:

a τ (w h , u h ) = a(w h , u h ) + Ωe∈ Ω(v • ∇w h , τ e (Lu h -f )) Ωe = (w h , f ) ∀w h ∈ V h (11)
Concerning the choice of the stabilizing parameter τ e , we refer to Hughes et al. in [START_REF] Hughes | Advances in Stabilized Methods in Computational MechanicsThe variational multiscale method-a paradigm for computational mechanics[END_REF] and we choose:

τ e = h e 2||v|| 2
coth(P e Ωe ) -1 P e Ωe [START_REF] Hachem | Unified adaptive Variational MultiScale method for two phase compressible-incompressible flows[END_REF] where h e is the local mesh size, computed as described in [START_REF] Hachem | On the stabilized finite element method for steady convection-dominated problems with anisotropic mesh adaptation[END_REF]. P e Ωe is the local Peclet number defined as:

P e Ωe = ||v|| 2 h e 2a (13)
Remark 1. In order to study the error of this numerical scheme, we use the Variational Multiscale analysis. In fact, it has to be noted that the SUPG scheme can be considered as a particular form of the generalized VMS formulation. Indeed, the stabilizing term of Eq. ( 11) can also be seen as the effect of the subscales on the coarse scales. Hughes gives more details about this concordance in [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF].

The Variational Multi-Scale formulation

The VMS formulation consists in decomposing the solution and test functions spaces into two sub-spaces: a mesh scale subspace (or coarse scales) (S h , V h ) and an under-mesh scale subspace (or subscales) (S , V ) such that S = S h ⊕ S and V = V h ⊕ V . Therefore, we can decompose the solution and test functions as follow:

u = u h + u , u h ∈ S h , u ∈ S w = w h + w , w h ∈ W h , w ∈ S (14)
Thanks to the orthogonality between the coarse scales subspace and the subscales subspace, the variational form can be split into a coarse scales sub-problem and a subscales sub-problem [START_REF] Hughes | Advances in Stabilized Methods in Computational MechanicsThe variational multiscale method-a paradigm for computational mechanics[END_REF]:

a(w h , u h ) + a(w h , u ) = (w h , f ) ∀w h ∈ V h a(w , u h ) + a(w , u ) = (w , f ) ∀w ∈ V (15)
We start by solving the subscales sub-problem (second equation). For smooth functions on the element interior but rough across the inter-element boundaries, the integration by parts leads to the following equation:

a(w , u ) = -a(w , u h ) + (w , f ) a(w , u ) = -(w , Lu h -f ) -(w , [Bu h ]) Γ -(w , Bu h ) Γ h ∀w ∈ V (16)
where the jump term [.] represents the difference of the fluxes on both sides of the element boundaries (see [START_REF] Hauke | The multiscale approach to error estimation and adaptivity[END_REF] for more details). An analytic solution of problem ( 16) can be found. This solution will be developed in the next subsection.

A posteriori error estimation on solution's subscales

In fact, the VMS analysis gives access to a solution of the subscales u . Considering the relation u = u h + u , this part of the solution can also be considered as the subscales error. However, u is not explicitly computed during the simulation. To obtain it, we use an a posteriori computation.

In this paper, we use two different methods for the a posteriori computation of the subscales error estimator. Both methods are residual based and rely on convection-dominated regime assumptions:

• The first method consists in computing the error estimator as the multiplication of the stabilizing parameter's norm and the residual's norm.

The computation is done element-wise because of the local definition of the stabilizing parameter and the local error norms used.

• The second method does not use the stabilizing parameter. The error estimator is computed thanks to the pointwise error estimation of Irisarri et al. in [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF]. It uses a set of bubble functions as a substitution of the subscales Green's functions.

In this section, we give the key features of these two computation methods.

Computation of the error estimator with stabilizing parameters

Going back to Section 2.3, we had the following sub-problem for the subscales:

a(w , u ) = -(w , Lu h -f ) -(w , [Bu h ]) Γ -(w , Bu h ) Γ h ∀w ∈ V
According to [START_REF] Hughes | Advances in Stabilized Methods in Computational MechanicsThe variational multiscale method-a paradigm for computational mechanics[END_REF], this problem can be solved analytically using Green's function g (x, y). We have the following paradigm:

u (x) = - Ωy g (x, y)(Lu h -f )(y)dΩ y - Γy g (x, y)([Bu h ])(y)dΓ y - Γ hy g (x, y)(Bu h )(y)dΓ y (17)
This error estimator includes:

• Interior residuals, Lu h -f in Ω • Inter-elements residuals, [Bu h ] on Γ • Natural boundary condition residual Bu h on Γ h
From here, multiple cases can be considered. In this work, we use linear shape functions with a piecewise constant source term f . Therefore, we can make the assumption that the residual Lu h -f is P0. Furthermore, we only have Dirichlet boundary conditions. Thus, we neglect the effects of the boundary condition residual. Finally, as for all convection-dominated regimes, we will first consider the smooth case, i.e. that the discontinuity of the subscales is not preponderant and that the inter-elements residuals can be neglected. With these assumptions, we write:

u (x) ≈ - Ωy g (x, y)(Lu h -f )(y)dΩ y (18)
We will also make the assumption that the error has a local influence. Thus, we can replace the subscales Green's function by a local Green's function that vanishes at the element boundary, i.e., g e (x, y) = 0 on Γ e . We can write:

g (x, y) ≈ g e (x, y) ≈ τ e δ(y -x) (19)
with,

τ e = 1 meas(Ω e ) Ωe,x Ωe,y g e (x, y)dΩ e,x dΩ e,y (20)
In [START_REF] Hauke | The multiscale approach to error estimation and adaptivity[END_REF], the author shows that the above assumptions work well for convectiondominated regimes. In fact, numerical experiments show that for high Peclet number flows and stabilized solutions, the contribution of jump terms is negligible.

From here, we can therefore write a local error estimator like the following:

u (x)| Ωe ≈ - Ωe g e (x, y)(Lu h -f )(y)dΩ y = - Ωe τ e δ(y-x)(Lu h -f )(y)dΩ y = -τ e (Lu h -f )(x) (21)
Taking the L 2 norm, we can finally write the following local error estimator for the solution subscales:

||u (x)|| L 2 ,Ωe = |τ e | ||Lu h -f || L 2 ,Ωe (22)
We know the expression of τ e from the application of the SUPG method and we can a posteriori compute the P0 residual in our finite element computation domain. As said in the introduction, the sub-scale information computed here is element-wise and can be used as such in mesh adaptation.

Computation of the error estimator with high order bubbles functions

In this section, we use the work of Irisarri et al. in [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF]. This time, we apply a pointwise computation of the error estimator. This computation method consists in decomposing the error into two components according to the nature of the residuals:

u (x) = u bub (x) + u poll (x). ( 23
)
The first term u bub is the internal residual error and it is related to the local internal residual, f -Lu h , inside the elements. As we will see later, this part of the error is modeled locally thanks to a set of bubble functions. The second term u poll is the inter-element error. It represents the pollution error due to sources of errors outside the element. As said in the previous section, it is negligible when considering convection-dominated regime. Consequently, in this paper we will consider only the internal residual error and suppose that:

u (x) ≈ u bub (x). (24)

Practical aspects of computation

As proposed by Irisarri et al. in [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF], the way of obtaining a numerical solution of u bub is by solving the following discrete problem:

Find u bub ∈ S h bub such that a(w bub , u bub ) = (w bub , f -Lu h ), ∀w bub ∈ S h bub ( 25
)
This time, the error component is expressed with a combination of bubble functions:

u bub (x) = n bub i=1 c b i b i (x) (26)
Considering bubbles functions of order 3, we have:

u bub (x) = c b 1 b 1 (x) + c b 2 b 2 (x) + c b 3 b 3 (x) (27)
with c b i unknown constant to be determined.

Referring to [START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF], the definition of the first bubble function b 1 (x) is the following:

b 1 (x) = (d + 1) d+1 d+1 i=1 λi ( 28
)
where d is the dimension of the problem and λi are the barycentric coordinates in the reference element.

The next bubble functions b 2 (x) and b 3 (x) are built by adding the monomials of the Pascal triangle with center in the barycenter c e = (ξ e , η e ) of the element. For example, in 2D, in the reference element:

Ω ref = {(ξ, η) : 0 ≤ ξ ≤ 1; 0 ≤ η ≤ 1 -ξ}, we choose the following bubble functions: b 1 (ξ, η) = 27 × ξη(1 -ξ -η) b 2 (ξ, η) = 27 × ξη(1 -ξ -η)(ξ -ξ b ) b 3 (ξ, η) = 27 × ξη(1 -ξ -η)(η -η b ) (29)
with c e = (ξ b , η b ), ξ b = 1/3 and η b = 1/3. Approximating u bub (x) by Taylor series and neglecting the second order terms, we have an expression of u bub (x) close to the centroid c i of the element [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF]: As said before, we make the assumption that the residual f -Lu h is P0. Therefore, we have:

u bub (x) ≈ b 1 (x)(f -Lu h )(c i )+b e y 1 (x) d(f -Lu h ) dy 1 | y=c i +b e y 2 (x) d(f -Lu h ) dy 2 | y=c i +HO (
d(f -Lu h ) dy 1 | y 1 =c i = d(f -Lu h ) dy 2 | y 2 =c i = 0 (32)
Consequently, the internal residual can be simply expressed as follows:

u bub (x) = b 1 (x)(f -Lu h )(c i ) (33)
Developing the residual with the convection-diffusion generic operator L, we have:

u bub (x) = b 1 (x)(f (c i ) + a∆u h (c i ) -v • ∇u h (c i )) (34)
Then, u h is P1, therefore, ∇u h (c i ) is a constant inside the element and ∆u h (c i ) = 0. For this reason, we finally get:

u bub (x) = b 1 (x)(f (c i ) -v • ∇u h (c i )) (35)
The above computation of the error estimator is pointwise. In fact, the error estimator is given at each point x of the domain. However, to include the error information in the mesh adaptation, we need an information of the error inside each element. To get this information, we can compute 3 different types of norm: L 1 , L 2 and L ∞ . We recall here the definition of these 3 norms:

||u bub || L 1 (Ω e ) = |Ω e | × 1 i N interp |u bub (x i )| ( 36
)
||u bub || L 2 (Ω e ) = |Ω e | × 1 i N interp |u bub (x i )| 2 (37)
||u bub || L ∞ (Ω e ) = max |u bub (x)| | x = (x i ) 1 i N interp (38)
where (x i ) 1 i N interp are the interpolation points defined by the integration points of the bubble function b 1 (x) and N interp is the number of interpolation points.

Validation of the error estimator computation with bubbles functions

Now that we have developed practical computation aspects, we validate the computation on a benchmark of a convection-diffusion equation. To do so, we take the example given in [START_REF] Irisarri | A posteriori pointwise error computation for 2-D transport equations based on the variational multiscale method[END_REF] of an error estimation where the FEM is stabilized by the SUPG scheme. We recall that the differential operator is given by: Lu = -a∆u + v • ∇u with a = 0.03, f = 1, v = (1, 1) and homogeneous Dirichlet boundary conditions.

We compute the SUPG solution on P1 triangular elements instead of the Q2 bilinear quadrilateral elements used in the referred paper. Therefore, we choose a number of elements in order to get approximately the same SUPG solution in comparison with Irisarri et al.. We choose to divide each square cells in 8 triangles and we obtain a mesh of 512 triangular elements. We compare our SUPG solution with the reference on Figure 1. Then, we compute the error estimate ||u bub || L ∞ (Ω e ) as explain in the previous subsection. We compare our results with the referred paper on Figure 2. Despite the difference in graphic representations, we observe a good concordance between the computation of Irisarri et al. and our computation. In fact, we reach approximately the same maximal values inside the elements and the error is located in the same area for both computations. Therefore, we validate the computation of the error estimator with bubble functions. We will see in the next section how to use this information in mesh adaptation.

Mesh adaptation with the subscales error estimator

In this section, we propose 3 different types of mesh adaptation:

1. An isotropic mesh adaptation: This technique uses the previously described subscales error estimator. From the latter, we build an isotropic metric tensor H iso to adapt the mesh. 2. An anisotropic mesh adaptation (described in [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF]): This technique uses the anisotropic local interpolation error indicator. From the latter, we derived an anisotropic metric tensor H aniso to adapt the mesh. 3. A new anisotropic mesh adaptation: This technique uses a new anisotropic local error indicator that takes into account (i) the interpolation error indicator and (ii) the subscales error estimator. From the latter, we derive a new metric tensor H new aniso . As a reminder, we start by introducing the principles of anisotropic mesh adaptation and the construction of the anisotropic metric tensor H aniso . Then, the isotropic tensor H iso based on the subscales error estimator is defined in Section 4.2. Finally, a combination between the anisotropic local error indicator and the subscales error estimator is proposed and the new metric tensor H new aniso is defined in Section 4.3.

Principles of anisotropic mesh adaptation

To discretize our computational domain, we use anisotropic and unstructured meshes. It implies that we can stretch the elements in certain directions according to the solution features. To do so, we start by performing an error analysis on the mesh. Then, to correlate the error with the geometry, a metric field is defined. From this metric field, an anisotropic error indicator is defined and used as a functional for a re-meshing optimization problem.

Let us consider a certain triangulation Ω h . We can derive an upper bound of the approximation error using an interpolation error analysis in the L p norm. Referring to Almeida et al. in [START_REF] Almeida | Adaptive finite element computational fluid dynamics using an anisotropic error estimator[END_REF], this upper bound is expressed thanks to the recovered Hessian of the approximated solution u h . In fact, using P1 linear elements, we usually cannot compute directly the Hessian of the solution. Instead, we compute an approximation called the recovered Hessian matrix:

||u -u h || L p (Ω) ≤ C N -α Ω h ||H R (u h )(x)|| p L p (Ωe) (39)
where α ≥ 0, N Ω h the number of elements of the mesh, H R (u h )(x) the recovered Hessian matrix and C is independent of the element size.

To apply the re-meshing strategy, we build an equilateral tetrahedron in the metric space. As in [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF], it is defined at an arbitrary point P by the local metric field M:

M(P ) = 1 h 1 (P ) e 1 ⊗ e 1 + ... + 1 h d (P ) e d ⊗ e d (40)
with (e i ) i=1,d the eigenvectors of H R (u h (x)) and h i (P ) the mesh sizes in the e i directions.

However, the recovered Hessian matrix is not a metric because it is not positive definite. Therefore, we define the following metric tensor:

H aniso = RΛR T ( 41
)
where R is the orthogonal matrix built with the eigenvectors (e i ) i=1,d of H R (u h (x)) and Λ = diag(|λ 1 |, ..., |λ d |) is the diagonal matrix of absolute value of the eigenvalues of H R (u h (x)). This metric tensor can also be written as follow:

H aniso = RΛR T = |λ 1 |e 1 ⊗ e 1 + ... + |λ d |e d ⊗ e d (42)
Here, we want to align the mesh with the solution field u. It means that we want the error to be equi-distributed on each direction of the domain. To do so, the shape of each element have to be such that the local error is equal in any direction. It is equivalent to say that the local error is constant per element in the principal directions of curvature. Therefore, we have: In [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF], the authors define a minimization problem where the functional is the error indicator of Eq. ( 44). This optimization problem is expressed as follow:

|λ 1 |h 2 1 = ... = |λ d |h 2 d = cte (43)
                          
Find h Ωe = {h 1,Ωe , ..., h d,Ωe }, Ω e ∈ Ω h that minimizes the cost function:

F (h Ωe ) = Ωe∈Ω h (η Ωe ) p
under the constraint:

N Ω h = C -1 0 Ωe∈Ω h Ωe d i=1 1 h i,Ωe dΩ e
where C 0 is the volume of a regular tetrahedron and Ω h is the new triangulation.

(45) Then, this multidimensional optimization problem is replaced by a onedimensional optimization problem. The unknown is no more h Ωe but the mesh size corresponding to the max of the (h i,Ωe ) 1 i d . The way of dealing with anisotropy is via stretching factors (s i,Ωe ) 1 i d-1 defined between the mesh sizes (h i,Ωe ) 1 i d . The solution is given by Mesri et al. in [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF]. In 3D and L p norm, it is expressed thanks to the following theorem: Theorem 1. For d = 3, the optimization problem (45) has a unique solution and is given by:

               h 3,Ωe = β (2p+3) 3 C 1,Ωe Ωe C 2,Ωe dΩ e 1 2(p+3) h 2,Ωe = s 2,Ωe h 3,Ωe h 1,Ωe = s 1,Ωe s 2,Ωe h 3,Ωe with          s 1,Ωe = h 1 h 2 = |λ 2 | |λ 1 | 1 2 s 2,Ωe = h 2 h 3 = |λ 3 | |λ 2 | 1 2
where

C 1,Ωe = 3 p C 0 s 1,Ωe s 2 2,Ωe × ||u || L ∞ (Ωe) u T OL p × |λ 3 | p , C 2,Ωe = C -1 0 1 s 1,Ωe s 2 2,Ωe and β 1 2 3 (p+3) = N -1 Ω h Ωe∈Ω h    1 Ωe C 2,Ωe dΩ e 1 2 3 (p+3) Ωe C 2,Ωe 2p + 3 3 C 1,Ωe 1 2 3 (p+3) dΩ e   
Finally, the above solution gives the mesh sizes in the 3 directions that define the metric field on each element. These mesh sizes are computed with respect to a fixed number of elements N Ω h .

Isotropic mesh adaptation with the subscales error estimator

The first use of the subscales error estimator computed in Section 3 is for an isotropic mesh adaptation. Referring to [START_REF] Babuška | The Finite Element Method and Its Reliability[END_REF] and from a theoretical convergence point of view, we know that the error is linked to the local mesh size h according to the following relation:

||u || L ∞ (Ω) ≈ C.h k+1 ( 46
)
where C is independent of the element size. In our case, using linear elements, we have:

||u || L ∞ (Ω) ≈ C.h 2 (47)
To control the mesh size of our adaptive process, the user has to give a tolerance value ||u tol || L ∞ so that the local mesh size stays above an acceptable bound. In fact, this tolerance corresponds to the desired error that the user wants to obtain on the mesh. For example, in a pointwise error estimation, ||u tol || L ∞ is the tolerated error at each control point inside the element. Referring to [START_REF] Hauke | Mesh adaptivity for the transport equation led by variational multiscale error estimators[END_REF], we can write the tolerated pointwise error as a scalar positive value u T OL :

||u tol || L ∞ = u T OL (48)
Starting from here, we can consider three strategies to adapt the mesh. Indeed, we can decide to either:

1. uniformly distribute the desired local error norm over the old mesh or, 2. uniformly distribute the desired local error norm over the new mesh or, 3. uniformly distribute the desired pointwise error.

In this paper, we will use the 2 nd mesh adaptation strategy for the element-wise computation of the error estimator (i.e. the one using stabilizing parameters); and the 3 rd mesh adaptation strategy for the pointwise computation of the error estimator (i.e. the one using bubble functions).

For the 3 rd strategy, still referring to [START_REF] Hauke | Mesh adaptivity for the transport equation led by variational multiscale error estimators[END_REF], we write the relationship between the error and the tolerance in a uniform pointwise error distribution:

h new h 2 = u T OL u ( 49
)
where h new is the size of the new mesh (after mesh adaptation), u is the estimated error.

Then, we build the isotropic metric tensor:

H iso = RΛR T = |λ|e 1 ⊗ e 1 + ... + |λ|e d ⊗ e d ( 50
) with |λ| = 1 h 2 new = ||u || L ∞ (Ω e ) u T OL × 1 h 2 (51)
where R is the orthogonal matrix built with the eigenvectors (e i ) i=1,d of H R (u h (x)).

Here, the eigenvalues λ of the metric are equal in all directions. Therefore, the mesh will be adapted isotropically. However, we keep the eigenvectors of H R (u h (x)) in the definition of H iso . Thus, we keep the element's orientation prescribed by the recovered Hessian matrix. Finally, the re-meshing strategy is the same that the one described in Section 4.1. In particular, as we apply an isotropic re-meshing, the stretching factors (s i,Ωe ) 1 i d-1 between the mesh sizes (h i,Ωe ) 1 i d are equal to 1.

Combination of subscales error estimator with anisotropic mesh adaptation

The second use of the subscales error estimator computed in Section 3 is for an anisotropic mesh adaptation. In this section, we propose to combine both the coarse scales error indicator and the subscales error estimator. To do so, we derive a new metric that allows to take into account the anisotropic variations of the solution on the mesh but also relies on the subscales error estimator previously computed.

We first consider the previous anisotropic local error indicator define in Eq. ( 44) and recalled here:

η Ωe = d|Ω e | 1 p |λ d (x 0 )|h 2 d
Then, going back to the previous section, we write:

h 2 d,new h 2 d = u T OL ||u || L ∞ (Ω e )
Now, the unknown of the re-meshing problem is h d,new . In fact, we want the new mesh size to take into account the error estimator of the subscales. Thus, we propose a new anisotropic local error indicator:

η Ωe,new = d|Ω e | 1 p × |λ d (x 0 )| × ||u || L ∞ (Ω e ) u T OL × h 2 d,new (52)
From here, we can define the new anisotropic metric tensor as:

H new aniso = RΛR T = ||u || L ∞ (Ω e ) u T OL |λ 1 |e 1 ⊗ e 1 + ... + ||u || L ∞ (Ω e ) u T OL |λ d |e d ⊗ e d (53)
Doing so, we keep the anisotropic effects from the solution variations but we isotropically scale this effect by the subscales error estimator.

With this new error indicator, we solve the same optimization problem of Eq. 45 with the unknown h d,new . The optimal mesh is obtained exactly in the same way but using η Ωe,new as functional. Consequently, the optimization problem becomes:

                          
Find h Ωe,new = {h 1,Ωe,new , ..., h d,Ωe,new }, Ω e ∈ Ω h that minimizes the cost function:

F (h Ωe,new ) = Ωe∈Ω h (η Ωe,new ) p
under the constraint:

N Ω h = C -1 0 Ωe∈Ω h Ωe d i=1 1 h i,Ωe,new dΩ e
where C 0 is the volume of a regular tetrahedron and Ω h is the new triangulation.

(54) Following the same proof that Mesri et al. in [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF], the solution of this new optimization problem becomes: Theorem 2. For d = 3, the optimization problem (54) has a unique solution and it is given by:

               h 3,Ωe,new = β (2p+3) 3 C 1,Ωe,new Ωe C 2,Ωe dΩ e 1 2(p+3) h 2,Ωe,new = s 2,Ωe h 3,Ωe,new h 1,Ωe,new = s 1,Ωe s 2,Ωe h 3,Ωe,new with          s 1,Ωe = h 1 h 2 = |λ 2 | |λ 1 | 1 2 s 2,Ωe = h 2 h 3 = |λ 3 | |λ 2 | 1 2
where

C 1,Ωe,new = 3 p C 0 s 1,Ωe s 2 2,Ωe × ||u || L ∞ (Ωe) u T OL p × |λ 3 | p , C 2,Ωe = C -1 0 1 s 1,Ωe s 2 2,Ωe and β 1 2 3 (p+3) = N -1 Ω h Ωe∈Ω h    1 Ωe C 2,Ωe dΩ e 1 2 3 (p+3) Ωe C 2,Ωe 2p + 3 3 C 1,Ωe,new 1 2 3 (p+3) dΩ e   
By re-defining the optimization problem, we integrate the effect of the subscales error estimator with respect to the constraint on the fixed number of 20 elements. This new re-meshing strategy is tested on 2D and 3D benchmarks in the next section.

Numerical examples

In this section, we test the effectivity and the accuracy of our new local error indicator on 2D and 3D convection-diffusion benchmark cases in convection-dominated regime. To do so, we start by computing the Finite Element solution u h thanks to the SUPG scheme described in this paper. The solution is computed on 6 different types of mesh:

• Fixed mesh: Using no mesh adaptation,

• Isotropic mesh 1: Using the isotropic metric tensor H iso computed with the subscales error estimator using the stabilizing parameter (see Section 3.1),

• Isotropic mesh 2: Using the isotropic metric tensor H iso computed with the subscales error estimator using the bubble functions (see Section 3.2),

• Anisotropic mesh 1: Using the anisotropic metric tensor H aniso computed without the subscales error estimator,

• Anisotropic mesh 2: Using the new anisotropic metric tensor H new aniso computed with the subscales error estimator using the stabilizing parameter,

• Anisotropic mesh 3: Using the new anisotropic metric tensor H new aniso computed with the subscales error estimator using the bubble functions.

For the benchmark cases 5.1 and 5.3, we have a known analytic solution u. Thus, we can compute the exact discrete error e h = u-u h on each element of the mesh. With this local discrete error, we compute the global L 2 norm of the error such that:

||e h || L 2 = Ω e 2 h dΩ 1/2 (55)
For the case 5.1, we define an efficiency index I ef f for our subscales error estimator as follow:

I ef f = ||u || L 2 ||e h || L 2 (56)

Case 1: regular boundary layers in 2D

This first case has a continuous solution and regular boundary layers. It has been studied by several authors like Zhang et al. in [START_REF] Zhang | Finite element superconvergence on Shishkin mesh for 2-D convection-diffusion problems[END_REF] or Hachem et al. in [START_REF] Hachem | On the stabilized finite element method for steady convection-dominated problems with anisotropic mesh adaptation[END_REF]. We consider the domain Ω = (0, 1) 2 and the velocity field v(x, y) = (1, 1) T . Instead of a source term f = 1 like in Section 3.2.2, we choose the source term corresponding to the following analytic solution:

u(x, y) = xy 1 -e -1-x a 1 -e -1-y a (57)
Thus, we have the following source term:

f (x, y) = (x + y) 1 -e -1-x a e -1-y a + (x -y) e -1-y a -e -1-x a ( 58
)
We compute the Finite Element SUPG solution u h for two different diffusion coefficients: a = {10 -3 ; 10 -4 }. In addition, for practical purpose, we define the following stretching factor:

S Ωe = h 2 new h 2 = u T OL ||u || L ∞ (Ω e ) (59)
To keep an acceptable mesh even where ||u || L ∞ (Ω e ) is close to zero, we have to impose a lower bound for the stretching factor S Ωe . We decide that, from an iteration to the other, the size of the element should not increase of more than 10 times the previous one. Therefore, the condition: S Ωe 0.01 is applied.

The effectivity index for the two ways of computation of the error estimator is given on Figure 3. The first observation that can be made concerning the efficiency index I ef f is that for the two different Peclet numbers and for the different mesh sizes considered, our efficiency index I ef f stays close to 1. Therefore, we can say that our error estimator is effective for all configurations. Secondly, we observe that the error estimation computed with bubble functions is more effective than the one computed with stabilizing parameters which overestimates the error. Finally, we observe that the effectivity index I ef f gets closer to 1 when the diffusion coefficient is reduced (i.e. when the element Peclet number increases). This finding is in accordance with the convectiondominated regime assumption. Now that we have confirmed the effectivity of our error estimation, we can use this estimation in the mesh adaptation. To analyze the mesh convergence of u h , we compute u h for 6 different mesh sizes: 2 000, 20 000, 40 000, 60 000, 80 000 and 100 000 elements.

The distributions of the discrete error e h are presented on Figure 4. The corresponding meshes for each error distribution are given below. On this figure, the case corresponds to a diffusion coefficient of a = 10 -3 . Mesh To begin the analysis of the results on the error distributions, we compare the Figure 4.(a), (b) and (c). We notice that the isotropic mesh adaptation according to our subscales error estimator reduces the error in the boundary layer. In fact, we notice from Figure 4.(b) and (c) that the mesh is isotropically adapted in the boundary layer. However, we observe that the error stays above the prescribed tolerance u T OL = 0.01 when we get close to the wall. To counter this effect, anisotropic mesh adaptation is a really efficient technique. On Figure 4.(d), we observe that the anisotropic mesh adaptation makes it possible for the error to go under the tolerance. However, we still keep an error approximately equals to the tolerance on the right top angle.

To eliminate this last part of the error, the new local error estimator proposed in Section 4.3 is well adapted as one can see on Figure 4.(f). In fact, we notice from the comparison of these 6 figures that our new local anisotropic error indicator η Ωe,new is the best driver of mesh adaptation. We observe on Figure 4.(e) and 4.(f) a drastic reduction of the error in the boundary layers compared to the other mesh adaptation techniques. Therefore, the use of the sub-scale information into the new anisotropic error indicator allows to improve the locality of the mesh nodes with respect to the equi-distribution of different error scales. Furthermore, we see that the equi-distribution is optimal when the error estimator is computed with the bubble functions (see Figure 4.(f)). To finish, we notice that the constraint of a fixed number of elements is well respected with a deviation of less than 19% with the targeted number of elements N .

The reduction of the error can also be observed on the results for the global L 2 norm of the error ||e h || L 2 . This is true for the case a = 10 -3 on Figure 5 and the case a = 10 -4 on Figure 6. As expected, the isotropic mesh adaptation reduces the error for both cases of a. The error estimator computed with bubble functions shows to be the best driver of isotropic mesh adaptation; in particular when the number of element increases.

To increase the slope of convergence, we use anisotropic mesh adaptation. In fact, it has been shown in a number of works (for example in [START_REF] Mesri | Dynamic parallel adaption for three dimensional unstructured meshes: Application to interface tracking[END_REF]) that the Hessian based anisotropic mesh adaptation leads to a convergence of order 2. This fact is highlighted in our results as we observe the increase of the slope going from isotropic to anisotropic mesh adaptation. On Figure 5 and 6, we observe that our new anisotropic local error indicator η Ωe,new behaves as expected. This is the first important result of the present work. In fact, we observe that for coarse meshes (i.e. N = 2 000 elements), the effect of the subscales error estimator is important. This reduction is even slightly better for the case with a = 10 -4 . On the latter figure, we observe that, with only 20 000 elements, we are below the error of an interpolation based anisotropic mesh adaptation of 100 000 elements.

As the number of elements increases, we notice that the effect of the subscales error estimator becomes less important. Indeed, the slopes of the Anisotropic mesh 2 and Anisotropic mesh 3 cases decreases. This behavior is expected when we acknowledge that the subscales error estimator can also be considered as an estimation of the subscales error. Thus, when the number of element increases, the local mesh size decreases and, for this reason, the modeled part of the solution becomes smaller.

The second important result of the present work concerns the last (Anisotropic mesh 3 ) case corresponding to anisotropic mesh adaptation taking into account the error estimator computed with bubble functions. Again, we can state that this error estimator is the best choice to drive anisotropic mesh adaptation in comparison with the one computed with stabilizing parameters. In addition, we notice that the subscales error estimator computed with stabilizing terms has no effect when we reach an important number of elements (for example: N = 100 000).

Case 2: parabolic boundary layers with recirculating convective field

For this benchmark, we consider a domain Ω = (-1, 1) 2 with a nonconstant flow. The diffusion coefficient is a = 0.001. We apply a zero source term inside the domain. The right side wall has a Dirichlet boundary condition equals to 1. The other walls have a homogeneous Dirichlet boundary conditions. We apply a recirculating convective field defined by v = (2y(1 -x 2 ), -2x(1 -y 2 )) that will form boundary layers in the domain. The SUPG solution is given on Figure 7. This example represents the propagation of the right wall heat inside the domain through the convective field. For this complex problem, it is difficult to capture the structure of the boundary layer. Our goal here is to show the efficiency of our new anisotropic error indicator in the capture of the parabolic layer. To illustrate the latter, we plot the two different error estimators (one computed with stabilizing parameter and the other with bubble functions) on Figure 8. The aim of this figure is to compare the two isotropic components of the two different error estimators. Therefore, it is sufficient to plot only the subscales part of the error estimation. As a result, we observe that the localization of the error is improved by the bubble functions method.

Case 3: regular boundary layers in 3D

This 3D benchmark case is inspired by the one in Section 5.1. It takes the same configuration concerning velocity field, source term, boundary conditions and analytic solution with an extrusion in the z direction. However, the convection-diffusion equation is now solved on a cubic domain Ω = (0, 1) × (0, 1)×(0, 1). The analytic solution is given on Figure 10. We compute the Finite Element SUPG solution u h with the 6 mesh adaptation techniques previously describe for 9 mesh sizes: N = 40 000, 80 000, 120 000, 160 000, 200 000, 300 000, 400 000, 50 elements.

As for the 2D benchmark, the error distribution corresponding of each mesh is given on Figure 11. b) and (c), we observe the instabilities due to the highly convective regime. Having isotropic elements, the local mesh size in the direction of the velocity field is too important. As said in the introduction of this paper, we observe that anisotropic mesh adaptation is a really efficient technique to eliminate the spurious oscillations in the boundary layer. This is highlighted on Figure 11.(d). However, there is still a part of the error that is not reduced on this figure. To counter this effect, we see that our new anisotropic local error indicator η Ωe,new is, again, efficient. Also, we observe that the constraint of a fixed number of elements is again well respected with a deviation of less than 6% with the targeted N . Finally, the analysis on error distribution in 2D still holds in 3D and moreover, it seems to be enforce by adding a new dimension to the problem. The above statements on error distribution are also to be found in the global L 2 norm of the error ||e h || L 2 on Figure 12. Moreover, the effects of the subscales error estimator on anisotropic mesh adaptation is enforced in 3D. The important result here is that, with only 40 000 elements in 3D, the error is below an interpolation based anisotropic mesh adaptation on 600 000 elements with our new anisotropic local error estimator η Ωe,new . Again, taking into account the sub-scale information in mesh adaptation makes it possible to reduce drastically the error on the 3D coarse meshes. With about 20 times less elements, we obtain the same global error ||e h || L 2 .

Conclusions

We proposed a new a posteriori error estimator based on the Variational Multiscale method for anisotropic adaptive fluid mechanics problems. This new error estimator is defined with the combination of both (i) an interpolation based anisotropic error indicator and (ii) a subscales error estimator. The results show that this combination allows to obtain highly precise solutions with much less elements in comparison with other mesh adaptation techniques. Furthermore, the definition of a new re-meshing optimization problem allows us to include the sub-grid information in mesh adaptation with respect to the constraint of a fixed number of elements.

The subscales error estimator is computed with two different methods and the results show that one using bubble functions is better suited to drive both isotropic and anisotropic mesh adaptation. In fact, the localization of the residual subscales error is better established with the latter.

This work shows that this new anisotropic mesh adaptation strategy is capable of dealing with boundary layers of convection-diffusion problems associated with high Peclet number flows in 2D and 3D. It is therefore a promising idea to simulate complex CFD aerothermal problems.

  Ωe g e (x, y)(y 1 -c i,1 )dΩ y and b e y 2 = Ωe g e (x, y)(y 2 -c i,2 )dΩ y[START_REF] John | A robust SUPG norm a posteriori error estimator for stationary convection-diffusion equations[END_REF]

Figure 1 :

 1 Figure 1: Comparison of the SUPG solution between Irisarri et al. [35] (left) and our solution (right)

Figure 2 :

 2 Figure 2: Comparison of error estimate between Irisarri et al. [35] (left) and our error estimation (right)

Following 1 p

 1 the work of Mesri et al., we introduce the following local error indicator in the L p norm: η Ωe = d|Ω e | |λ d (x 0 )|h 2 d (44) where |λ d (x 0 )| is the maximum eigenvalues of H R (u h (x)) corresponding to direction d, |Ω e | is the volume of the element and h d is the length of the element in direction d.

Figure 3 :

 3 Figure 3: Effectivity indexes for a = 10 -3 and a = 10 -4

Figure 4 :

 4 Figure 4: Error distributions and meshes of 2D regular boundary layers for different mesh adaptation techniques

Figure 6 :

 6 Figure 6: Error in L 2 norm with a = 10 -4 for regular boundary layers in 2D

Figure 7 :

 7 Figure 7: Numerical SUPG solution for parabolic layers with a = 10 -3

Figure 9 :

 9 Figure 9: Resulting meshes for different types of mesh adaptation for the parabolic boundary layers.

Figure 10 :

 10 Figure 10: Analytic solution for 3D convection-diffusion equation with a = 10 -3

Figure 12 :

 12 Figure 12: Error in L 2 norm with a = 10 -3 for regular boundary layers in 3D

Acknowledgments

This work was supported by Safran Aircraft Engines [CIFRE 2015/1348].

adaptation is under a constraint of 20 000 elements and the figures are respectively: The resulting adapted meshes are presented on Figure 9. The cases (a), (b), (c), (d), (e) and (f) correspond to the same mesh adaptation techniques than in the previous example (Case 1). On Figure 9 (b) and (e) we observe that the isotropic and anisotropic mesh adaptation techniques follow the localization prescribed by the error estimator based on the stabilizing parameter (see Figure 8). Consequently, the mesh is adapted only on a restraint region of the parabolic boundary layer. On Figure 9(d), we observe that the conventional Hessian based anisotropic mesh adaptation already capture the 4 boundary layers. However, it does not take into account the sub-grid information that can improve the anisotropic mesh adaptation. In fact, the mesh adaptation technique using the new anisotropic local error indicator taking into account the subscales error indicator based on bubble functions (see Figure 9(f)) shows to be the best choice to accurately capture the 4 parabolic boundary layers.


Maxim Kolomeets
email: kolomeec@comsec.spb.ru

Andrey Chechulin
email: chechulin@comsec.spb.ru

Igor Kotenko

Visualization Model for Monitoring of Computer Networks Security based on the Analogue of Voronoi Diagrams

Keywords: visual analytics, visualization of security data, graphical models, computer networks, Voronoi diagrams

In this paper we propose an approach to the development of the computer network visualization system for security monitoring, which uses a conceptually new model of graphic visualization that is similar to the Voronoi diagrams. The proposed graphical model uses the size, color and opacity of the cell to display host parameters. The paper describes a technique for new graphical model construction and gives examples of its application along with traditional graph based and other models.

Introduction

Computer networks are rapidly growing today. Meanwhile, the more devices are in the network, the harder it is to ensure its security. This problem is met by operators of security systems of corporate level (e.g., security information and event management systems, SIEM systems), when the analyzed computer network is measured not only by hundreds of employees' workplaces and high order technical equipment, but also by smart doors, servers, various sensors of climate, security, etc.

To cope with the control of growing networks we need to apply systems for monitoring network security, which give us possibility to visualize the computer network and parameters of its state in a simple and efficient manner. But, as a rule, in such systems the network is presented with the application of rather traditional graphical models, for example, in the form of graphs or tables, that are difficult to understand in the case of large networks and display of a variety of parameters. In order to cope with this problem it is necessary to improve the efficiency of visualization means by complex use of various graphical models such as graphs, matrices, treemaps, parallel coordinates, etc. in the framework of the multiple view concept [START_REF] Wang | Guidelines for Using Multiple Views in Information Visualization[END_REF]. At the same time it is necessary to increase the efficiency of visualization of particular graphic models.

In the case of visualization of computer networks and their security, different techniques are developed that allow clustering of segments of the network (e.g., based on clustering of graph elements) or encapsulation of the state parameters [START_REF] Shi | Scalable network traffic visualization using compressed graphs[END_REF]. Yet another solution to the problem is to develop conceptually new graphical models, which are able to present information in a form that is new for the user and that allows to increase the efficiency of the user's work.

The novelty of this paper is to use a conceptually new graphic visualization model similar to the Voronoi diagrams, which allows to increase the effectiveness of visual analysis for the computer network security, for example, as one of functions of the SIEM system. It is expected that this model will be used in the developed visualization system in the framework of the supported multiple view concept.

The main contribution of the paper lies in the fact that it offers a new technique of visualization of network security, as well as reveals the theoretical and practical side of how this technique can be used in SIEM systems. The organization of the paper is as follows. Section 2 analyzes existing graphical models that can be used to visualize parameters of computer network security with description of their advantages and disadvantages. In section 3 we describe the developed conceptually new graphical model. Section 4 discusses the developed system for visualization of computer network security and provides examples of application of the proposed graphical model in the framework of this system. In section 5 and 6 the proposed graphical model is evaluated, as well as its comparison with other graphical models is performed. Section 7 discusses conclusions and future research directions.

Review of computer network visualization techniques

Within the developed computer network visualization system, graphs, matrices, and treemaps were mainly used to realize the multiple view concept. The listed graphical model (Fig. 1) have different advantages and disadvantages [START_REF] Tufte | Visual Explanations[END_REF], which can also be expressed in terms of informativity and ease of perception and use. Informativity can be represented as the detailization level, which is expressed in the compliteness of simultaneously displayed data. Ease of perception and use can be represented in the form of speed and simplicity of user interactions with displayed data. It is obvious that different graphical models have different ratio of informativity and ease of perception and use. Thus, often increase of the informativity due to the resulting congestion of the graphical model negatively affects the ease of perception and use, and vice versa. Some examples are an informative but difficult to read table, and an uninformative, but easy to understand semaphore shown in Fig. 2.

Fig. 2. Examples of a table and a semaphore

Let us consider from this point of view the graphic models [START_REF] Klyshinskij | Review of the methods of multidimensional data visualization[END_REF][START_REF] Marty | Applied Security Visualization[END_REF], which are most often used to visualize the security of computer networks. Matrices [START_REF] Kwan-Liu | Cyber Security Through Visualization[END_REF] (Fig. 1, right) are efficient in displaying the relations of elements of a small computer network that has complex topology and where each host has many connections. Security parameters, determined based on network traffic, can be set using the color and transparency of cells, located at intersections of rows and columns. However, the size of the cells depends on the dimension of the matrix, which is set by the number of network hosts. With the increase in the number of rows and columns of the matrix, the size of a particular cell goes to one pixel, the perception of the color tone and the more transparency of which is difficult. Inefficient use of space of matrices should be noted, when most of the cells remains blank, which has negative effect on the users' perception of parameters, especially when large computer networks are visualized. It is also worth considering that the matrices can visualize the links parameters, but not the parameters of the computer network hosts.

On the other hand, matrices efficiently display clusters of network, and they can be used to construct attack graphs [START_REF] Noel | Understanding Complex Network Attack Graphs through Clustered Adjacency Matrices[END_REF] (Fig. 3, left) and to analyze entire segments of the network, rather than individual hosts.

Matrices can also be expanded by displaying multidimensional data in 3D space. For example, the 3D analogue of a matrix is the dispersion chart [START_REF] Lau | The spinning cube of potential doom[END_REF], which is represented as a cube which axes are local IP addresses, global IP addresses and port numbers, and the color of the dot shows successful and unsuccessful attempts to establish TCP connections (Fig. 3, right). Thus, the presence of long lines or planes, consisting of the points on the dispersion chart, can inform about network scanning.

For efficient visualization of host parameters (for example, the data of the vulnerability scanners [START_REF] Kwan-Liu | Cyber Security Through Visualization[END_REF] (Fig. 4, left), assessment of criticality of assets, etc.) of a computer network, one can also use treemaps [START_REF] Harrison | Nessus Vulnerability Visualization for the Web[END_REF] (Fig. 1, in the center). Treemaps are efficient to display parameters, as they can handle the color, depth and size. It should be noted that in the presence of cells of large area, one can also use transparency as an additional parameter. However, treemaps are suitable to visualize purely hierarchical networks and are unable to visualize the interaction parameters. However, along with matrices, treemaps can be expanded, for example, for constructing attack graphs [START_REF] Williams | GARNET: A Graphical Attack Graph and Reachability Network Evaluation Tool[END_REF], when the steps of the intruder are displayed by directed edges (Fig. 4, right). Graphs [START_REF] Mcguffin | Simple Algorithms for Network Visualization: A Tutorial[END_REF], as the most common way to visualize computer networks, are efficient for topology visualization and can display parameters of hosts using the vertices of the graph, as well as interaction parameters using edges. However, like matrices, graphs inefficiently use space, leaving large empty areas. However, graphs are often used to display the network topology, when host type is used as vertices (Fig. 5).

For visualization of parameters vertices of graphs can be replaced by glyphs [START_REF] Novikova | Analytical Visualization Techniques for Security Information and Event Management[END_REF]. Glyphs can be represented in a pie chart in which the number of equal sized segments depends on the number of displayed parameters (Fig. 6). The parameters themselves can be expressed in the form of segments' color and their transparency. The glyphs can be augmented by the ring with the same number of segments that display the previous parameter value. Due to this, it is possible to produce a historical analysis.

The disadvantage of glyphs is that when they are used the graphical model becomes overloaded, thereby ease of perception deteriorates. It is worth saying that the graphs have many variations of how they are built and used. For example, graphs can be used for visual analysis of granting access to resources [START_REF] Montemayor | Information Visualisation for Rule-based Resource Access Control[END_REF] with the role access system that supports hierarchies or groups of users, to visualize patterns of network traffic [START_REF] Glatz | Visualizing Big Network Traffic Data using Frequent Pattern Mining and Hypergraphs[END_REF], for visualizing logs of financial transactions [START_REF] Novikova | Analytical Visualization Techniques for Security Information and Event Management[END_REF], for visualization of computer attacks [START_REF] Mansmann | Visual Support for Analyzing Network Traffic and Intrusion Detection Events using TreeMap and Graph Representations[END_REF][START_REF] Kotenko | Common Framework for Attack Modeling and Security Evaluation in SIEM Systems[END_REF], etc. At the analysis of existing graphic models the features were outlined that allow the user to efficiently analyze information. First, the user perceives better spatial mapping (plane and spacial figures), while the color and shape of figures are optional parameters. It is easy to demonstrate: the most important security parameters in treemaps (Fig. 4) are displayed with size of the planes; graphs (Fig. 5 and Fig. 6) operate with the sizes of the vertices; and at rendering in the form of a matrix (Fig. 3) the operator most often searches for and analyzes structures in the form of planes and lengthwise lines, consisting of individual cells. Second, it can be concluded that to visualize the links parameters it is better to use matrices, but they are unable to visualize hosts. For visualization of hosts parameters it is better to use treemaps, but they are unable to visualize links. In the case when it is necessary to visualize both hosts and links one needs to use graphs. Thus, the analysis of the advantages and drawbacks of existing graphical models showed that the simultaneous display of parameters of hosts and links is only possible with the use of graphs. However, if the vertices in the graphs are displayed in the form of planes, the operator will be able to analyze the information faster. Thus, analysis of relevant works showed that the creation of graphical models that will allow to effectively display both parameters (as in the treemaps) and the topology (as in graphs) is a perspective approach.

3

The proposed graphical model

For monitoring of computer network security a graphical model is proposed [START_REF] Kolomeec | Methodological Primitives for Phased Construction of Data Visualization Models[END_REF], which visually resembles Voronoi diagram [START_REF] Aurenhammer | Voronoi Diagrams and Delaunay Triangulations[END_REF], however, it is not the same from a mathematical point of view (Fig. 7). The main idea is to integrate capabilities of graphs to visualize the topology of the computer network and treemaps to visualize parameters. The solution appeared from the representation of network hosts in the form of cells, and links between hosts in the form of links between these cells. At the same time in the graphical model there are separators (dark grey in Fig. 7) which divide the cells that are next to each other, but have no links. For ease of understanding, we can present an analogy in the form of a maze: cells-polygons that represent hosts can be interpreted as the maze rooms; the links of the cells that represent the relationship between hosts can be interpreted as doors between the maze rooms; the separators that represent the lack of links can be interpreted as the maze walls. The algorithm for constructing the proposed graphical model is more complex than the algorithms for constructing graphs, matrices or treemaps, and consists of four steps: (1) building of the convex hull of a given planar graph; (2) implementation of the restricted Delaunay triangulation [START_REF] Aurenhammer | Voronoi Diagrams and Delaunay Triangulations[END_REF]; (3) formation of cells, based on triangulation; (4) selection of separators.

Let us consider the algorithm for constructing the proposed graphical model in more detail, on the example of the implemented software tool that provides a visual interface to display the security parameters of computer networks.

The proposed graphical model is implemented on the basis of the graph adjacency matrix. The first step builds a planar graph, which is supplemented by the convex hull. The convex hull is required for obtaining the convex figure, which is used to perform the next step. Graph which will be used in this example, and the result of the first step are shown in Fig. 8 and Fig. 9 respectively. In the next step for the resulting figure we should produce a restricted Delaunay triangulation [START_REF] Aurenhammer | Voronoi Diagrams and Delaunay Triangulations[END_REF]. It is worth noting the importance of implementing exactly limited triangulation, as it allows to triangulate the figures taking into account already existing relations and to avoid crossing of edges. The result of this step is shown in Fig. 10. The next step is to form the cells that will serve as the basis of the graphical model. For this we need to associate a subset of the triangles, obtained as the result of triangulation, with the corresponding vertex of the graph. For each vertex of the host (in Fig. 11 it is selected as a light grey circle) we find a subset of triangles (in Fig. 11 they are highlighted in gray) that includes this vertex. Next for triangles of this subset the weight centers are determined (in Fig. 12 they are shown as light gray points), union of which gives the desired polygon (in Fig. 12 it is highlighted in light gray edges). The resulting polygon corresponds to the host, based on which we defined the subset of triangles of the triangulation. The result of the step is shown in Fig. 13, where the edges of the figure resulting from the triangulation are black, and the edges of the desired cells are red (in Fig. 13 it is light grey). The next step is to outline the separators. Since each cell corresponds to a specific host, the edges of the cells-hosts, with which there is a link, can be designated a certain color, for example gray. All other edges will be separators and will have an appropriate color (e.g. red in Fig. 14). To obtain a figure of a certain shape (e.g. rectangle as in Fig. 14) we can also add points to the cells of the convex hull, or move these points to the required positions. In Fig. 14 the resulting shape is depicted with addition of new points, and Fig. 15 shows a figure with relocation of the common points of the polygons-hosts of the convex hull. It should be noted that from the mathematical point of view [START_REF] Aurenhammer | Voronoi Diagrams and Delaunay Triangulations[END_REF], the proposed model is not a Voronoi diagram, despite the visual similarity, since the Voronoi diagram has the different mathematical meaning, and existing algorithms for its construction do not allow to visualize the computer network topology when constructing the chart based on vertices of the graph.

Examples of application of the proposed graphical model

The proposed graphical model can be used to visualize the security parameters of a computer network or to analyze the behavior of the attacker. It may be required when informing the operator of SIEM system about the threats of a security breach or by visual analysis of computer network security. Let us consider examples of application of the proposed graphical model in the framework of implementation of the visualization system in more details.

Description of the visualisation system

To analyze the security of a computer network the visualization system is developed, which supports the display of the computer network using both the classical methods of visualization, such as graphs, graphs augmented with glyphs, treemaps and matrices and the graphical model, proposed in this paper. The example of dashboard of the developed visualization system is depicted in Fig. 16. The system includes the ability to manage data sources, aggregation and correlation of data collected from sources and tools for visual analytics of computer network security. Tabs to navigate to the relevant controls are at the top of the dashboard in Fig. 16.

In the left part of the dashboard in Fig. 16 there is enumeration of representations of computer networks formed on the basis of data from different sources.

The central part of the dashboard is the implementation of the multiple view concept, when the data is displayed in different views: as graph; as the graph augmented with glyphs; charts and diagrams; treemaps; matrices and as graphical model, formed as the analogue of Voronoi diagrams, presented in this paper.

In various usage scenarios, the user selects the graphical model, which is capable to visualize data the user needs at the moment most efficiently. In the following subsections the examples of scenarios of using the developed graphical model in comparison with the graphs will be presented.

Description of the source data

Data of a computer network (Fig. 17), which consists of 9 segments, will be used as source data for visualization.

Fig. 17. Segments of a computer network

The segment, which consists of external users who have remote connection to the computer network via the Internet, is represented in block 1. Block 2 displays the web server for remote connectivity, as well as hosts needed to operate the web server. Block 3 represents the security system located between the web server and the demilitarized zone (DMZ). The demilitarized zone is shown by block 4 and includes segments of the internal network of the company. Block 5 and block 9 are the computers of the internal network users. Block 6 displays the devices connected to the network through Wi-Fi connection. Block 7 is the server of storing and processing data. Block 8 corresponds to the virtualization server, together with the virtual machines placed on it.

Example 1. Visualization of the state of computer network security

To visualize the security of a computer network for each host the indicators of protection from attacks and possible damage in case of compromise of this host are calculated. In the proposed graphical model the security level against attacks can be represented in the form of a polygon color -hosts that have passed the threshold are in red color (in Fig. 18 -dark grey), and possible damage is shown as the size of the polygon (area of polygon). We shall also consider the example of visualization based on the graph, where the security is represented by vertex color, and the possible damageas the radius of the vertex. Let us consider two variants.

In the first case (Fig. 18) the vulnerable hosts are red (dark grey in Fig. 18) are strongly scattered. Almost every network segment has vulnerable hosts. Despite the fact that damage at their compromise a small, scatteredness gives greater variability of actions for the attacker due to the presence of many potential hosts from which attack may occur.

Fig. 19 shows the corresponding visualization based on the graph, where the value of the possible damage in case of compromise is outlined by radius of vertex, and the presence of vulnerability is marked in red color (in Fig. 19 light gray, as they are almost invisible, they are indicated by arrows). It is obvious that the proposed graphical model allows us to more quickly identify vulnerable hosts and produce visual analysis of the damage done when they are compromised.

In another case (Fig. 20), it is clear that virtualization platform is under vulnerability, and therefore all machines, located on it, as well as some computers of internal users are also under vulnerability. Corresponding visualization based on the graph is shown in Fig. 21. As in the previous case, visual analysis of potential damage and identifying the vulnerable segment are more efficient when using the proposed model.

Example 2. Visualization of the attack route

Visualization of the route of the attack (Fig. 22) may be noted as another example of using the proposed graphical model. The host from which an attacker carries out an attack is one of the computers of internal users. This host can be denoted in blue (dark gray in Fig. 22). All hosts, to which the actions of the attacker were recorded, can be also displayed in blue, however with different degree of transparency, which will depend on the intensity of the actions of the attacker. Thus it is possible to analyze which segment of the network is affected by the attack and if the attacker compromised the most important hosts or not. One can also estimate how close the attacker came to certain hosts for subsequent selection of protection strategies. In Fig. 22 it is seen that in such a scenario only some elements of the network will be affected (virtualization platform and servers for storing and processing data), however the potential damage from compromise is unreasonably great.

Evaluation of the proposed graphical model

In order to evaluate the proposed graphical model two groups of indicators were identified -performance indicators and functionality indicators. Functionality indicators define the set of scenarios that the proposed graphical model is able to visualize, such as network connectivity, network size, abilities to display topology, possibility to display the parameters of hosts and links between hosts. Performance indicators define the efficiency of user's perception of information and ease of working with data, for example: the efficiency of indicators perception, efficiency of finding a way to attack, efficiency of analysis of network segments.

The evaluation of the proposed graphical model was carried out in comparison with the visualization of a computer network as a graph.

The evaluation was performed using a survey of experts. All experts note that in a scenario, when the computer network is a non-planar graph, application of the proposed graphical model is impossible, which is the main drawback of the proposed graphical model. In some scenarios, when non-planar graphs are rarely used or not used at all, the proposed graphical model is by an order more efficient. It is also noted that the proposed graphical model is not efficient in visualization of small computer networks, and to visualize them it is more efficient to use matrices and graphs.

Thus, the graphical model is most appropriate to be used when rendering a medium to large scale computer networks. High computational complexity of the algorithm is also a disadvantage, however this can be shortened by using a client-server architecture, when only coordinates of the cells will be transferred to the thin client.

As a whole experts agree in opinion that the proposed graphical model has more options for visualizing metrics and network segments in comparison with the classical methods of visualization that are based on graphs, matrices and treemaps.

Thus, the proposed graphical model is an alternative to visualization of computer networks in the form of graphs, treemaps and matrices. The use of human spatial perception (the location of the cell-computers relative to each other) and the absence of necessity of edges provides a number of advantages at the cognitive level of perception of visualization. On the other hand, a disadvantage of the presented graphical model is that it can help to visualize exclusively planar graphs, which reduces the usage scope. However, in some scenarios, when graphs are not planar, are rarely used or not used at all, the proposed graphical model is by order more efficient. It is also worth noting that the use of the proposed graphical model is applicable for visualization of any object that can be represented as a planar graph.

It is supposed that the proposed graphical model can improve the efficiency of visual analytics using the graphical models to visualize the computer networks within the multiple view concept in SIEM systems.

Discussion

Based on the evaluation of the proposed graphical model a comparative 4) fine -graphical model supports this method of visualization, the efficiency of the perception is high and the user easily analyzes the information. Graphs are not restricted in display of size and color of vertices, size and color of connections, clustering, and can also display any topological types. Graphs have a limit on the transparency of the elements, as with the transparency of less than 30% the item will be hard to read. Graphs also can display hosts and connections by a limited set of geometric shapes to be visualized in 3D, provided that clustering of the vertices is done. The only thing that graphs do not support is display of nesting. Treemaps do well with display of hosts using the size and color of planes. It is the use of planes that allows the user to efficiently analyze information. With it, treemaps have no restrictions in the display of hosts with transparency, because even at minimal transparency the outlines of the plane are kept, and the user cannot miss it, as it would be in the case of graphs. However, treemaps cannot display the parameters of the links and cannot be displayed in 3D. But at the same time they do not have restrictions on clustering of hosts and perfectly display nesting, as it is the basis of the very concept of the treemaps. However, because of the nesting treemaps allow to display exclusively hierarchical networks.

Matrices cannot display hosts, but allow to visualize them when using the cells' colors and their transparency. Matrices cannot operate with cells' sizes. Transparency of the matrix cells have the same restriction as graphs (at least 30% transparency), but they have no restrictions in colors. Cells themselves can also be represented as a limited set of geometric shapes and can be displayed in 3D. Matrices cannot visualize nesting of the hosts, however, with the help of links they can efficiently display the network's clusters that will be represented in the form of planes (Fig. 1, right part). Matrices can display without limitations any topological types, however they are most efficient in finding clusters in complex non-planar networks.

The proposed graphical model uses the size, color and opacity of the cell for displaying of parameters of the host, as in the trees of maps, so their estimations coincide. To display the links' parameters it uses size, color and transparency of edges of cells. Edges are visualized as lines like in graphs, therefore their estimations are also identical, except for the transparency -in the proposed model, even at zero transparency the contour of lines is preserved. The proposed graphical model can be repre-sented in the form of a 3D polyhedron, where the edges will correspond to the cells. The graphical model also has the capability of clustering, as shown in Fig. 18, and through the use of planes, inside which one can place similar planes, it supports nesting. The proposed graphical model also supports the display of hierarchical and planar networks, but does not support non-planar network.

Thus, we can see that the proposed model in some cases (in the hierarchical and planar network topologies) can provide an alternative to graphs, treemaps or matrices, or be used as a supplement to graphical models within the multiple view concept, based on which the dashboard of security systems are built.

The development of the presented graphical model is been continued. At the moment the development of the algorithm of the polymorphism of cells is performed, in order to be able to change the size and shape of the cells without violating the topology. The development of the algorithm of display of proposed graphical model in 3D is carried out, this can be achieved by imposing points of the graphical model on the sphere, and then to draw cross-sections of the sphere at the points that correspond to the cells.

Another direction for future research is to analyze the possibility and efficiency of using the proposed graphical model for visualization of processes associated with information security, but which are not associated with computer networks. Any process that is currently visualized using planar graphs (some examples are displayed in section 2), can be visualized with the proposed graphical model. Thus, despite the fact that in the implemented system the proposed model is used to find vulnerabilities, risk assessment and other parameters of a computer network, the scope and possibility of application to ensure information security are much wider.

Conclusion

In this paper the analysis of existing visualization methods was performed and the new graphical model based on the analogue of Voronoi diagrams was presented. We demonstrated that the proposed graphical model of visualization of computer networks allows in some cases to display data more efficiently, compared to already existing graphical models. The developed visualization system was presented used to analyze the security of computer networks, and examples of visual analysis of the computer network state were provided. The estimation of the proposed model for visualizing the security parameters of computer network was done, and the comparison of its efficiency with graphs, matrices and treemaps was performed. Directions of future research based on the use of the proposed graphical model were presented, in particular the development of algorithms for polymorphism and nesting of cells, display in 3D and using the proposed model for visual analytics of processes and objects, which were previously represented as graphs.

Fig. 1 .

 1 Fig. 1. Examples of graphical models (left to right): graphs, treemaps and matrices

Fig. 3 .Fig. 4 .

 34 Fig. 3. Usage of a matrix to visualize attack graphs (left) and as the dispersion chart (right)

Fig. 5 .Fig. 6 .

 56 Fig. 5. A graph displaying host type

Fig. 7 .

 7 Fig. 7. The proposed graphical model

Fig. 8 . 9 .

 89 Fig. 8. Planar graph Fig. 9. The construction of the convex hull

Fig. 10 .

 10 Fig. 10. The graph and his triangulation

Fig. 11 .Fig. 12 .

 1112 Fig. 11. The determination of triangles required to build a cell

Fig. 13 .

 13 Fig. 13. The result of building the cells of the graphical model

Fig. 14 . 15 .

 1415 Fig. 14. Adding new points to the polygons Fig. 15. Relocation of points of the polygons Therefore, the graphical model allows to display hosts in the form of planes, and the links between the hosts -in the form of contact planes. The formal description of the algorithm to build the presented graphical model can be represented as the following pseudo code: // step 1 -making the convex hull of planar graph graph = getPlanarGraph(adjacencyMatrix) graph = addConvexHullToGraph(graph) //step 2 -making the constrained Delaunay triangulation listOfTriangles = makeСonstrainedDelaunayTriangulation(graph) listOfPolygons //step 3 -making cells based on triangulation for each dot from graph { listOfTrianglesThatContainsDot = getTriangles(dot, listOfTriangles); polygon for each triangle from listOfTrianglesThatContainsDot{ triangleCenter = getCenterOfTriangle(triangle) addPontToPolygon(polygon, triangleCenter) }} //step 4 -making separators for each polygon from listOfPolygons{ for each line from polygon{ connectedPolygon = getConnectedPolygon(line, polygon) connectionIsExist = isConnectionExist(polygon, connectedPolygon, graph) if connectionIsExist == true{ color = grey }else{ color = red } setLineColor(line, color) }} // the end of algorithm is drawing of cells draw(listOfPolygons)

Fig. 16 .

 16 Fig. 16. The dashboard of the developed visualization system

Fig. 18 .Fig. 19 .Fig. 20 .

 181920 Fig. 18. Scattered unprotected hosts are presented on the basis of the proposed model

Fig. 21 .

 21 Fig. 21. Two segments of the network with unprotected hosts, presented on the basis of the graph

Fig. 22 .

 22 Fig. 22. Visualization of the attack route

table (

 (

	tab. 1)
	was built with three other graphical models -graphs, treemaps and matrices. It should
	be noted that graphical models are considered in the minimum version, i.e. without
	additions and excluding the presentation of information in the form of text and signa-
	tures. The table shows the display capabilities and the number of simultaneously dis-
	played parameters of security of network hosts, network connections, possibility of
	extension of graphical models for visualization in 3D space, possibility of visualiza-
	tion of topological parameters and possibility of display of networks of different topo-
	logical types. The cells contain evaluations of the efficiency of parameters display
	perception in one way or another. Four estimations are outlined: (1) does not support -
	graphical model does not support this method of visualization; (2) supports -graph-
	ical model supports this method of visualization, but with restrictions; (3) good -
	graphical model supports this method of visualization; (

Table 1 .

 1 Comparison of the possibilities of graphs, treemaps, matrices, and the proposed model

			Graphs	Treemaps	Matrices	The proposed model
	Hosts pa-rameters display	size color transparency	good good supports	fine fine good	does not support does not support does not support	fine fine good
		size	good	does not support	does not support	good
	Display of parameters of links	color transparency	good supports	does not support does not support	good supports	good good
		shape	supports	does not support	supports	supports
	Possibilities of extension for display in 3D	supports	does not support	supports	supports
	Display of topology parameters	hosts clasterization incapsulation of hosts (nestness)	good does not support	good fine	fine does not support	good fine
	Display of different topological types	hierarchical planar non-planar	good good good	fine does not support does not support	good good fine	fine fine does not support

Acknowledgements.

This research is being supported by the Ministry of Education and Science of The Russian Federation (contract 14.604.21.0137, unique contract identifier RFMEFI60414X0137) in SPIIRAS.