\documentclass[CRMATH,Unicode,biblatex,XML]{cedram} \TopicFR{Analyse numérique} \TopicEN{Numerical analysis} \addbibresource{CRMATH_Hoch_20230672.bib} \usepackage{amssymb} \usepackage{mathtools} \graphicspath{{Figures/}} \newcommand{\divrm}{\mathrm{div}} \newcommand\n{\nabla} \newcommand{\dd}{\mathrm{d}} \newcommand\RR{\mathbb{R}} \newcommand{\density}{\mathrm{density}} \newcommand{\energy}{\mathrm{energy}} \newcommand{\sumrm}{\mathrm{sum}} \newcommand{\defrm}{\mathrm{def}} \renewcommand\div{\n\cdot} \renewcommand\O{{\mathcal O}} \renewcommand\S{{\mathcal S}} \newcommand\F{\mathcal{F}} \newcommand\U{\mathcal{U}} \newcommand\V{{\mathcal V}} \newcommand\mytextcolor[2]{#2} \newenvironment{hypotheses}{\begin{enonce}{Hypothesis}}{\end{enonce}} \let\tfrac\frac %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \makeatletter \let\save@mathaccent\mathaccent \newcommand*\if@single[3]{% \setbox0\hbox{${\mathaccent"0362{#1}}^H$}% \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}% \ifdim\ht0=\ht2 #3\else #2\fi } %The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath: \newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna} %If there's a superscript following the bar, then no negative kern may follow the bar; %an additional {} makes sure that the superscript is high enough in this case: \newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}} %Use a separate algorithm for single symbols: \newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}} \newcommand*\wide@bar@[3]{% \begingroup \def\mathaccent##1##2{% %Enable nesting of accents: \let\mathaccent\save@mathaccent %If there's more than a single symbol, use the first character instead (see below): \if#32 \let\macc@nucleus\first@char \fi %Determine the italic correction: \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}% \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}% \dimen@\wd\tw@ \advance\dimen@-\wd\z@ %Now \dimen@ is the italic correction of the symbol. \divide\dimen@ 3 \@tempdima\wd\tw@ \advance\@tempdima-\scriptspace %Now \@tempdima is the width of the symbol. \divide\@tempdima 10 \advance\dimen@-\@tempdima %Now \dimen@ = (italic correction / 3) - (Breite / 10) \ifdim\dimen@>\z@ \dimen@0pt\fi %The bar will be shortened in the case \dimen@<0 ! \rel@kern{0.6}\kern-\dimen@ \if#31 \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}% \advance\dimen@0.4\dimexpr\macc@kerna %Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows: \let\final@kern#2% \ifdim\dimen@<\z@ \let\final@kern1\fi \if\final@kern1 \kern-\dimen@\fi \else \overline{\rel@kern{-0.6}\kern\dimen@#1}% \fi }% \macc@depth\@ne \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}% \macc@set@skewchar\relax \let\mathaccentV\macc@nested@a %The following initialises \macc@kerna and calls \mathaccent: \if#31 \macc@nested@a\relax111{#1}% \else %If the argument consists of more than one symbol, and if the first token is %a letter, use that letter for the computations: \def\gobble@till@marker##1\endmarker{}% \futurelet\first@char\gobble@till@marker#1\endmarker \ifcat\noexpand\first@char A\else \def\first@char{}% \fi \macc@nested@a\relax111{\first@char}% \fi \endgroup } \makeatother \let\oldbar\bar \renewcommand*{\bar}[1]{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}} \let\oldtilde\tilde \renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}} %\let\tilde\widetilde %\let\hat\widehat \let\oldhat\hat \renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}} %\let\tilde\widetilde \renewcommand*{\to}{\mathchoice{\longrightarrow}{\rightarrow}{\rightarrow}{\rightarrow}} \let\oldmapsto\mapsto \renewcommand*{\mapsto}{\mathchoice{\longmapsto}{\oldmapsto}{\oldmapsto}{\oldmapsto}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} \hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red} \newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}} \newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}} \newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}} \newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% \title{An induced limitation in the reconstruction step for Euler equations of compressible gas dynamics in arbitrary dimension} % \alttitle{ Une limitation induite pour l'étape de reconstruction dans les équations d'Euler en dimension quelconque} \author{\firstname{Philippe} \lastname{Hoch}} % \address{CEA-DAM, DIF, 91297, Arpajon Cedex, France} \email{philippe.hoch@cea.fr} \subjclass{00X99} \keywords{\kwd{Compressible Euler system} \kwd{arbitrary order reconstruction} \kwd{induced admissible limitation}} \altkeywords{\kwd{Équations d'Euler compressibles} \kwd{reconstruction d'ordre arbitraire} \kwd{limitation induite admissible}} \begin{abstract} % anglais We are interested in the limitation process for the reconstruction of quantities related to Euler's equations of compressible gas dynamics for a general pressure law of type $P(\rho,\epsilon)$ (density, specific internal energy). For example, for perfect gas laws, we recall the constraints $\rho>0$ and $\epsilon>0$, and that the velocity $\mathbf{U}$ is \emph{a priori} not bounded in the continuous problem. Nevertheless it is in $L^2(\Omega,\rho)$ as a consequence of the relation on the energies $E=\epsilon + \frac{1}{2} |\mathbf{U}|^2$ in $L^1(\Omega,\rho)$ (due to global conservation of total energy $\rho E$). We show a similar result for conservative reconstruction in any space dimension and for an arbitrary reconstruction order. The use of the Leibniz formula on the specific variables $\epsilon$, $\mathbf{U}$ and $\frac{1}{2}|\mathbf{U}|^2$ allows to obtain also such a discrete \emph{induced} control of reconstructed velocity thanks to control of reconstructed density and energies. We build a direct limitation on the weight variable $\rho$ and also especially on the specific variable $\epsilon$. In particular, the latter makes it possible to limit, in an \emph{induced} way, the velocity $\mathbf{U}$. The limited reconstruction of the conservative variables is deduced from the assembly of these different limitation processes. We illustrate in dimension $d=1$ and $d=2$ on some test cases, our reconstructions of orders 2 and 3. \end{abstract} \begin{altabstract} On s'intéresse au processus de limitation pour la reconstruction des quantitées liées aux équations d'Euler de la dynamique des gaz compressibles pour une loi de pression générale de type $P(\rho,\epsilon)$ (densité, énergie interne massique). Par exemple, pour la loi des gaz parfait, les contraintes sont $\rho>0$ et $\epsilon>0$, et la vitesse $\mathbf{U}$ n'est \emph{a priori} pas bornée dans le problème continu. Elle est néanmoins dans $L^2(\Omega,\rho)$ comme conséquence de la relation sur les énergies $E=\epsilon + \frac{1}{2} |\mathbf{U}|^2$ dans $L^1(\Omega,\rho)$ (par conservation globale de l'énergie totale $\rho E$). On montre un principe similaire dans le cadre d'une reconstruction conservative en dimension d'espace quelconque et pour un ordre de reconstruction lui aussi arbitraire. L'utilisation de la formule de Leibniz sur les variables massiques $\epsilon$, $\mathbf{U}$ et $\frac{1}{2}|\mathbf{U}|^2$ permet en effet d'obtenir en discret un contrôle \emph{induit} de la vitesse reconstruite grâce au contrôle des reconstructions de la densité et des énergies. Nous construisons une limitation directe sur la variable poids $\rho$ et aussi surtout sur la variable massique $\epsilon$. En particulier, cette dernière permet de limiter, de manière \emph{induite}, la vitesse $\mathbf{U}$. La reconstruction limitée des variables conservatives se déduit de l'assemblage de ces différents processus de limitation. Nous illustrons sur des cas tests, en dimension $d=1$ et $d=2$, les reconstructions d'ordres 2 et 3 ainsi obtenues. \end{altabstract} \begin{document} \maketitle \section{Introduction/Framework} \mytextcolor{red}{ We are interested in a novel limitation process for finite volume schemes. We propose to limit some variable in order to fulfill some principles (physical or invariant domain validity) in such a way that it induce the limitation of some other variables. Here, we apply this strategy in the case of compressible Euler system of gas dynamics. In this case, we focus on the construction of the velocity limitation induced by the limitation of specific internal energy (and density)}. These equations model the conservation of mass, momentum and total energy, and write: %\newpage \begin{equation} \label{eq:cons_laws} \begin{array}{llllllllllllllll} \partial_t\boldsymbol{\U}(t,\mathbf{x}) +\div(\boldsymbol \F(\boldsymbol \U(t,\mathbf{x}))) = \mathbf{0} \\ % \bm{\U}(0,\mathbf{x})= \bm{\U}_i, \qquad \mbox{if}\:\: \mathbf{x} \in \Omega_i.% \noindent \end{array} \end{equation} where $\mathbf{x}=(x_1,\ldots ,x_d)$ and: \begin{equation} \begin{gathered} \boldsymbol \U = (\rho, \rho \boldsymbol U, \rho E), \qquad \boldsymbol \F (\boldsymbol \U) = (\rho \boldsymbol U, \rho \boldsymbol U \otimes \boldsymbol U + P I_d, \boldsymbol U ( \rho E + P) ) \\ \boldsymbol U \:\:\:\:\: {\rm the \:\:\: velocity}, \qquad %:\:\:\: % {\rm and} \\ E \:\:\:\:\: {\rm the \:\:\: specific \:\:\: total \:\:\: energy\:\:\: with} \\ \label{eq:rel_tot_nrj}%eq:rel_alg_nrjtot} E = \epsilon + \frac{1}{2}|\boldsymbol U|^2, \:\:\:\: {\rm and} \:\:\:\: \epsilon \:\:\:\: {\rm the\:\:\: specific \:\:\: internal \:\:\: energy}. \end{gathered} \end{equation} \mytextcolor{blue}{ In the following, vectors are noted in bold. Let $\mathbf{a}$ and $\mathbf{b}$ be two vectors in $\RR^d$, their dot product is denoted by $\left(\mathbf{a},\mathbf{b}\right)$ and tensorial product by $\mathbf{a} \otimes \mathbf{b}$. For our purpose the pressure $P$ is considered as $P(\rho, \epsilon)$, where $\rho \in I_{\density}$ and $\epsilon \in I_{\energy}$ ($I_Q$ is an open interval of $\RR$, for exemple for perfect gas EOS~\eqref{eq:EOS_GP}: $I_{\density}=I_{\energy}=]0,+\infty[$):} \begin{equation} \label{eq:EOS_GP} P = (\gamma-1) \rho \epsilon, \qquad (\gamma>1 : {\rm adiabatic \:\: constant}). \end{equation} \mytextcolor{blue}{ Note however that this work is valid for general pressure laws $P=P(\rho,\epsilon)$ (as far as~\eqref{eq:cons_laws} is hyperbolic). Consider a two point finite volume discretization of~\eqref{eq:cons_laws}, let's define some notations. % \begin{notation} \begin{itemize} \item $\Omega_j$ is a cell, $|\Omega_j|$ his volume. $\partial \Omega_j$ denote the cell boundary, and $f$ is a $d-1$ dimensional hypersurface ($f$ is an edge for $d=2$, a face for $d=3$, etc, $f \subset \partial \Omega_j$), $|f|$ the associated $d-1$ measure, and the cell index $k$ will denote the neighboring cell sharing $f$ with $j$. \item $r$ is a node, $k_i(r)$ design the set of cells containing the node r. \item in dimension $d=2$, $r+1/2$ will also denotes an edge (see also Figuree~\ref{Fig: V_j}). \end{itemize} % \end{notation} The so called finite volume unknown is the cell average $\boldsymbol \U_{j}(t)\coloneqq \frac{1}{|\Omega_j|} \int_{\Omega_j} \boldsymbol{\U}(t,\mathbf{x}) \dd x$, we refer e.g.~to~\cite{Toro,GoRav}: \begin{align} \label{first_order_transport_direct_euler_std_cont} \frac{\dd }{\dd t} \boldsymbol \U_{j}(t) + \frac{1}{|\Omega_j|}\sum_{f \in \partial \Omega_j} \int_{f} \F (\boldsymbol \U(t,\mathbf{x})).\mathbf{N} \dd s &= 0, \quad\mbox{let}\:\: \V_{j}\:\: \mbox{an approximation of} \: \U_j, \\ \label{first_order_transport_direct_euler_std} \frac{\dd }{\dd t} \boldsymbol \V_{j}(t) + \frac{1}{|\Omega_j|}\sum_{f \in \partial \Omega_j} |f| G(\boldsymbol \V_{j},\boldsymbol \V_k, \tilde{\mathbf{N}}_f) &= 0. \end{align} } $\tilde{\mathbf{N}}_f$ is an exact or averaged unit normal to $f$. % (in case of $f$ is not an hyper-plan of $\mathbf{R}^d$). The function $G$ in~\eqref{first_order_transport_direct_euler_std} is a numerical flux that is locally conservative ($G(a, b, \mathbf{N}) + G( b, a, -\mathbf{N})=0$) and consistent ($ G( a, a, \mathbf{N})= \F (a).\mathbf{N}$), for example: Roe, Rusanov, HLL, Relaxation, VFFC... \\ \mytextcolor{blue}{ We consider a high-order spatial reconstruction, represented by a (high degree polynomial) function $\mathcal{R}_{j}(\mathbf{x},\boldsymbol \V)$. We also need to deal with a high-order quadrature flux formula $(\omega_l^f,\mathbf{x}^f_l)_{l=1}^{s}$ (e.g. Gauss Legendre), so that we consider: \begin{equation} \label{high_order_transport_direct_euler_std} \frac{\dd }{\dd t} \boldsymbol \V_{j}(t) + \frac{1}{|\Omega_j|} \sum_{f \in \partial \Omega_j} \sum_{l} \omega_l^f |f| G\left(\mathcal{R}_{j}\left(\mathbf{x}^f_{l},\boldsymbol \V\right),\mathcal{R}_{k}\left(\mathbf{x}^f_{l},\boldsymbol \V\right), \tilde{\mathbf{N}}_f\left(\mathbf{x}^f_{l}\right)\right) = 0. \end{equation} In case of dimension $d=2$ (then from~\eqref{first_order_transport_direct_euler_noeud} to~\eqref{high_order_transport_direct_euler_composite}), let us also consider the extension of numerical fluxes to nodes then to composite finite volume schemes (see~\cite{PH_NodalSolv:hal-03585115}). In the latter, the numerical fluxes are defined at all co-dimension object (from $1$ to $d$) on quadrature points of $\partial \Omega_j$.}%For example, in two dimensions, These are either localised at nodes $G_j^r$, or at some interior edge points $G_j^{r+1/2}$. We recall that these nodal and composite numerical fluxes $G_j^r, G_j^{r+1/2}$ are locally conservative and consistent (and that the edge part $r+1/2$ may be any of the classical numerical edge flux of~\eqref{high_order_transport_direct_euler_std}). The associated local normals at these points are defined by vectors $\mathbf{C}_j^{r}=\frac{1}{2}(\boldsymbol{x}_{r+1}-\boldsymbol{x}_{r-1})^{\perp}$, $\mathbf{C}_j^{r+1/2}=(\boldsymbol{x}_{r+1}-\boldsymbol{x}_r)^{\perp}$ and for \mytextcolor{blue}{edge invariant} quadrature formula: \begin{equation} \label{first_order_transport_direct_euler_noeud} \hspace{-.3cm} \frac{\dd}{\dd t} \boldsymbol \V_{j}(t) + \frac{1}{|\Omega_j|} \left( (1-\theta) \sum_{r \in \partial \Omega_j} G_j^{r}(\boldsymbol \V_{k_1(r)},\ldots ,\boldsymbol \V_{k_m(r)}) \:\: \mathbf{C}_j^{r} + \theta \hspace{-.2cm} \sum_{r+1/2 \in \partial \Omega_j} G_j^{r+1/2}(\boldsymbol \V_{j},\boldsymbol \V_{k})\:\: \mathbf{C}_j^{r+1/2} \right) = 0. \end{equation} Here again, we consider extension to higher order in both quadrature and polynomial reconstruction of~\eqref{first_order_transport_direct_euler_noeud} \begin{hypotheses}{\label{hypothesis_order_quad}} \mytextcolor{blue}{ We assume that there exists $\{ q_i, \theta_i \}_{i=1}^s$ some interior quadrature parameters and weights ($\theta_i \ge 0$) such that:} \begin{equation} \label{eq:2Dcomp_souscvx} 1-\sum _{i=1}^s\theta_i \ge 0. \end{equation} Let $q_0=0$ and $q_{s+1}=1$ be the extremity parameters, and the associated extremity weights: \begin{equation} \label{eq:2Dcomp_node_ext} \theta_0=\theta_{s+1}=\frac{1}{2}\left(1-\sum _{i=1}^s \theta_i\right). \end{equation} \end{hypotheses} Let us now consider, the arbitrary order composite (in terms of geometric objects of non-unique co-dimension fluxes type) %(in terms of non-unique co-dimension geometrical fluxes) scheme~\cite{PH_NodalSolv:hal-03585115} (here, in dimension $d=2$): \begin{multline} \label{high_order_transport_direct_euler_composite} \frac{\dd }{\dd t} \boldsymbol \V_{j}(t) + \frac{1}{|\Omega_j|} \left( (1 -\sum _{i=1}^s \theta_i) \sum _{r \in \partial \Omega_j} G_j^{r}(\mathcal{R}_{k_1(r)}(\mathbf{x}_{r},\boldsymbol \V),\ldots ,\mathcal{R}_{k_m(r)}(\mathbf{x}_{r},\boldsymbol \V))\:\: \mathbf{C}_j^{r} + \right. \\ \hspace{-.3cm}\left. \sum _{i=1}^s \theta_i \sum _{e \in \partial \Omega_j} G_j^{r+1/2}(\mathcal{R}_{j}(\mathbf{x}^e(q_i),\boldsymbol \V),\mathcal{R}_{k}(\mathbf{x}^e(q_i),\boldsymbol \V))\:\: \mathbf{C}_j^{r+1/2} \right) = 0. \end{multline} \begin{itemize} \item Taking $s=1$ in~\eqref{eq:2Dcomp_souscvx}~\eqref{eq:2Dcomp_node_ext} %\eqref{high_order_transport_direct_euler_composite} gives $q_1=\frac{1}{2}$, $\theta_1=\theta \in [0,1]$ so recovering~\eqref{first_order_transport_direct_euler_noeud} (we may reach the third/fourth order with $\theta=\frac{2}{3}$). \item \mytextcolor{blue}{Taking $(\theta_i,q_i)=(\omega_i,q_i), \forall i=1,\ldots ,s$ (Gauss Legendre) in~\eqref{eq:2Dcomp_souscvx}~\eqref{eq:2Dcomp_node_ext} so that $\sum _{i=1}^s \omega_i =1$ and then by~\eqref{eq:2Dcomp_node_ext} $ \theta_0=\theta_{s+1}=0$ %\eqref{high_order_transport_direct_euler_composite} so recovering arbitrary high-order pure edges schemes~\eqref{high_order_transport_direct_euler_std} (see e.g.~\cite{DumbserMunz}).}%, ($\sum _{i=1}^s \omega_i =1$). \end{itemize} %\input{introduction_limitation_processes.tex} \mytextcolor{red}{ A major numerical problem is the preservation of the set of admissible states. We recall that for usual equations of state, these are nothing but the positivity of density and temperature which, in many cases, corresponds to positivity of specific internal energy $\epsilon$. } \mytextcolor{red}{ One of the main difficulty is due to the fact that such invariant domain are expressed using the definition of the specific internal energy (or temperature), which is a nonlinear function of conservative variables (density, momentum, and total energy). Nevertheless, these latter are of prime importance to approximate correct weak (and entropic) solutions. Therefore, a limitation process of the specific energy leads to complex modifications of the conservative variables. In the literature, the question of the ``good'' variables to limit is not really addressed and this is the purpose of our contribution. %Our paper deals with a control of the reconstruction step under stability criteria. Up to our knowledge, for $d$-dimensional Euler gas-dynamic equations, no previous works propose a direct algebraic procedure based on non-linear reconstructions (non polynomial) of specific quantities (using Leibniz formula) of arbitrary order.} %Moreover, we also show that this process produce the same limitation coefficient betweenkinetic energy and internal energy inside the reconstruction of total energy. \mytextcolor{blue}{ The litterature on limitation processes for hyperbolic conservation laws is very rich and is still an active research domain.} \mytextcolor{red}{For finite volume methods it may be either included in the reconstruction (of averages unknowns $\V_j$) process or directly applied to the reconstructed flux $G$ (with blending technics).} \mytextcolor{blue}{Here, we adopt the philosophy of the former strategy: limitation of high-order reconstruction. It can be divided into two classes: the a-priori and a-posteriori processes.} \mytextcolor{blue}{ The latter includes APITALI~\cite{ale2d_strategy,HL_APITALI_Vecteur}, MOOD~\cite{MOOD} processes. These methods act by systematicaly %some criteria checking process of testing the numerical solution under some criteria with an a-posteriori limitation. The process first selects ``bad'' cells, then the main action is to decrease for this set some high-order parts (in compact or in a hierarchical way). Due to the existence of a low order (generaly first order) scheme fullfiling by construction the given criteria, the iterative loop always converge by choosing appropriate decreasing sequences of degree's reconstruction (see~\cite{ACPHNS}). % The iterative loop always converge on the total mesh (with update of bad set). this is due to the existence of a (generaly) low order scheme fullfiling by construction the given criteria. } \mytextcolor{blue}{The first class includes all the developpment based on the seminal constructive work of Harten based on \emph{a priori} TVD limitation criteria in dimension $d=1$. Unfortunately, due to a rather negative result of~\cite{GoodLevq} implying that a TVD scheme is at most first order in dimension $d \ge 2$, many efforts and contributions have been done on relaxed stability constraints. Among all the works, one may cite UNO/ENO/WENO methods~\cite{HO,HARTEN19973,RAbgrall,Weno}, maximum principle when the continuous PDE satisfies such a property~\cite{ShuLB,VdKas,QuadStab}, and invariant domain preservation~\cite{GMPT,ShuMaireVilar,ZhangShu12}.} \mytextcolor{blue}{The paper is organised as follow, we present an arbitrary high-order reconstruction for volumic (density), but also and more important of the specific (primary) quantities ($\mathbf{U}$ and E). The latter is based on arbitrary order multi-variate Leibniz formula. The second section is devoted to a non-linear reconstruction of internal energy, we show that this reconstruction depends on both high-order velocity terms and high-order internal energy terms.} \mytextcolor{red}{The high-order terms is a made of three contributions, each of them possessing different properties of homogeneity. %: the velocity is 2-homogeneous and the internal energy is linear. By imposing equality between the different contributions, we obtain, as a by-product, a natural dependence between limiters, thus inducing those of velocity reconstruction.} %Each of such high-order terms possesses such properties: velocity is 2-homogeneous and internal energy is linear. %1-homogeneous. %This is then obvious to make a link %Making a same level of contribution between these two scales, we obtain an induce control of high-order velocity terms from those of high-order internal energy own data. %. More specificaly, we show that this reconstruction depends on both high-order velocity tearms and high-order internal energy itself. % with respect to both high-order velocity terms (of order two) and high-order terms from his own data (of order one). %Using the relation linking specific total energy %and velocity reconstructions. In order to express high-order kinetic energy in terms of high-order velocity we re-use Leibniz. %We finaly obtain a non-linear reconstruction of internal energy that posses some homogeneous properties with respect to both %high-order velocity terms (of order two) and high-order terms from his own data (of order one). %A natural way of making a same scaling is then to adopt a quadratic dependency between this two terms. %This directly induce a link between high-order velocity and high-order internal energy. \mytextcolor{blue}{The third section is devoted to the construction of the so called induced limitation process: a scalar limiter for internal energy is proposed, for which a limiter is automaticaly deduced (and not specifically designed!) for the velocity. In the last section, we asses this induced limitation strategy on some discontinuous solution in dimension $d=1$ and $d=2$.} %\subsection{Nonlinear \emph{a priori} strategies for the reconstructions step} \subsection{Obtaining arbitrary order for volumic and specific variables} %\section{Obtaining arbitrary order for volumic and specific variables} We deal with arbitrary order reconstruction for a pure volumic field Q (e.g. density). If $n+1$ is the order (e.g. same order as the quadrature formula in~\eqref{high_order_transport_direct_euler_std} or in hypothesis~\ref{hypothesis_order_quad} for~\eqref{high_order_transport_direct_euler_composite}), we write: \begin{equation} \label{eq:reconst_prod_dens_var_volumic} \mytextcolor{red}{ P(\mathbf{x},Q)= \bar{Q}_j + (\nabla Q)_j (\mathbf{x}-\mathbf{x}_j) + \frac{1}{2} ((\nabla^2 Q)_j (\mathbf{x}-\mathbf{x}_j),(\mathbf{x}-\mathbf{x}_j)) +\cdots + \frac{1}{n!} (\nabla^n Q)_j \overbrace{(\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j)}^{n\text{ times}}.} \end{equation} Here $\mathbf{\nabla} = ^t(\partial_{x_1},\ldots ,\partial_{x_d})$ denotes the gradient operator, $(\nabla^n Q)_j$ is a n-multilinear form acting on $\RR^{n \times d}$, and $\mathbf{x}_j=\frac{1}{|\Omega_j|} \int _{\Omega_j} \mathbf{x} \, \dd \mathbf{x} $ is the centro\"id of cell $\Omega_j$ (or a point inside $\Omega_j$). Each $(\nabla^l Q)_j$ must be an approximation of $(\nabla^l Q) (\mathbf{x}_j)$, and $\bar{Q}_j$ is given by \begin{equation} \label{eq:cte_highorder_varQ} \bar{Q}_j = \begin{cases} Q_j - \sum _{l=1}^{n} c^l_j (Q), &\mbox{if conservation is mandatory,} \\ Q_j, &\mbox{otherwise}, \end{cases} \end{equation} with \:\: $ Q_j=\frac{1}{|\Omega_j|} \int _{\Omega_j} Q \dd \mathbf{x} $ \:\: and \begin{equation} \label{eq:cte_correct_highorder} \forall l=1,\ldots ,n \qquad c^l_j(Q) \coloneqq \frac{1}{l!} \frac{1}{|\Omega_j|}\int_{\Omega_j} (\nabla^l Q)_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j) \,\dd \mathbf{x}. \end{equation} \mytextcolor{blue}{ We recall that high-order derivatives in~\eqref{eq:reconst_prod_dens_var_volumic} may be represented by the tensorial form using canonical basis $e_{i} \in \RR^d$ = $(0,\ldots 0,\underbrace{1}_{i^{th}},0\ldots,0) $ and with the Einstein sum convention reads as: \begin{align} \label{eq:diff_tensor_form} (\nabla^n Q) &= (\partial_{x_{i_1}} \ldots \partial_{x_{i_n}} Q) \:\: (e_{i_1} \otimes\cdots \otimes e_{i_n}) \intertext{or} %\\ \label{eq:diff_tensor_form_a} (\nabla^n Q) (\mathbf{a},\ldots ,\mathbf{a}) &=\sum_{i_1=1}^d\ldots \sum_{i_n=1}^d (\partial_{x_{i_1}} \cdots \partial_{x_{i_n}} Q) \:\: a_{i_1}\ldots a_{i_n} \qquad \mbox{with}\:\: \mathbf{a}=\mathbf{x}-\mathbf{x}_j. \end{align} Notice that $ (\partial_{x_{i_1}}\ldots \partial_{x_{i_n}} Q) \prod _{j=1}^n a_{i_j}$ can be written as $ \prod _{i=1}^d(\partial_{x_{i}} ^{m_i} Q)\prod _{i=1}^d a_{i}^{m_i}$ with $m_i= \# \{j ; i_j=i \}$. Using the multi-index notation, let $\mathbf{m}$ be a list of $d$ natural numbers $\mathbf{m}=(m_1,m_2,\ldots ,m_d)$, $|\mathbf{m}|=\sum _{i=1}^d m_i$. Using% the convention that $\mathbf{m}$ being a list of $d$ natural numbers $\mathbf{m}=(m_1,m_2,\ldots ,m_d)$: \begin{equation} \label{eq:multi_index_notation} \partial^{\mathbf{m}} \:\:\:\: \mbox{corresponds to} \:\:\:\: \overbrace{\partial_{x_1}\ldots\partial_{x_1} }^{m_1}\overbrace{\partial_{x_2}\ldots\partial_{x_2} }^{m_2}\: \cdots \: \overbrace{\partial_{x_d}\ldots\partial_{x_d} }^{m_d}, \end{equation} we can rewrite~\eqref{eq:diff_tensor_form_a} \begin{align} \label{eq:diff_tensor_form_a_final} (\nabla^n Q) (\mathbf{a},\ldots ,\mathbf{a}) &= \sum _{|m|=n} \binom{n}{\mathbf{m}} (\partial^{\mathbf{m}} Q)\prod _{i=1}^{d} a_i^{m_i}, \qquad \mbox{with} \:\:\binom{n}{\mathbf{m}} = \frac{n!}{m_1!\ldots m_d!} = \frac{n!}{\prod _{i=1}^d m_i!}. % &=\sum _{|m|=n} \left( \inom{n}{\mathbf{m}} \prod _{i=1}^{d} a_i^{m_i} \partial^{\mathbf{m}} Q. \end{align} } \begin{remark} \mytextcolor{blue}{ We also may use other multivariate calculus by using the following notation: \begin{equation} \label{eq:diff_multi_variate_dotprod} \forall n\ge1 \:\:\:\:\: (\nabla^n Q) (\mathbf{a},\ldots ,\mathbf{a}) = [\mathbf{a}, \nabla]^n Q = \left( \sum _{i=1}^d (a_i \partial_{x_i}) \right)^n Q, %\qquad (\mbox{where} \:\: \mathbf{a}\coloneqq \mathbf{x}-\mathbf{x}_j). \: \:\: \mbox{where} \:\:\mytextcolor{red}{ \partial_{x_k} a_l = \delta_{kl} \:\: \mbox{(Kronecker symbol).}} \end{equation} Here the right hand side of~\eqref{eq:diff_multi_variate_dotprod} represents a non symetric linear application ($[A,B] \ne [B,A]$) of the vector direction $\mathbf{a}$ with the vector gradient $\nabla$ to the power $n$, the resulting differential operator being applied to $Q$ (see the two examples below). Using the multinomial Newton formula, the operator $[\mathbf{a}, \nabla]^n$~\eqref{eq:diff_multi_variate_dotprod} writes} \begin{equation} \label{eq:diff_multi_variate_dotprod_mindices} %[\mathbf{a}, \nabla]^n \left( \sum _{i=1}^d (a_i \partial_{x_i}) \right)^n = \sum _{|\mathbf{m}|=n} \binom{n}{\mathbf{m}} \prod _{i=1}^{d} (a_i \partial_{x_i})^{m_i},% \:\:\:\: \mbox{where} \binom{n}{\mathbf{m}} = \frac{n!}{m_1!...m_d!} = \frac{n!}{\prod _{i=1}^d m_i!}, \end{equation} and we finally get for $n\ge1$ (see also~\eqref{eq:diff_tensor_form_a_final}): \begin{equation} \label{eq:diff_multi_variate_dotprod_mindices_final} (\nabla^n Q) (\mathbf{a},\ldots ,\mathbf{a}) =\sum _{|\mathbf{m}|=n} \binom{n}{\mathbf{m}} \prod _{i=1}^{d} a_i^{m_i} \partial^{\mathbf{m}} Q . \end{equation} \end{remark} \begin{exemple} For example, if $n=2$, and $d=2$, the right hand side of~\eqref{eq:diff_tensor_form_a_final} (or~\eqref{eq:diff_multi_variate_dotprod_mindices_final}) applied to $Q$~gives \begin{equation*} %(a_1 \partial_{x_1} + a_2 \partial_{x_2})^2 Q (\nabla^2 Q) (\mathbf{a},\mathbf{a}) = a_1^2 \partial^2_{x_1} Q + 2 a_1 a_2 \partial_{x_1}\partial_{x_2} Q + a_2^2 \partial^2_{x_2} Q . \end{equation*} which is nothing but the tensorial form (with Hessian) in~\eqref{eq:reconst_prod_dens_var_volumic} (with $\mathbf{a}=\mathbf{x}-\mathbf{x}_j$). \end{exemple} \begin{exemple} Taking $n=3$, and $d=2$, the right hand side of~\eqref{eq:diff_tensor_form_a_final} (or~\eqref{eq:diff_multi_variate_dotprod_mindices_final}) applied to Q gives: \begin{equation*} %(a_1 \partial_{x_1} + a_2 \partial_{x_2})^3 Q (\nabla^3 Q) (\mathbf{a},\mathbf{a},\mathbf{a}) = a_1^3 \partial^3_{x_1} Q + 3 a_1^2 a_2 \partial^2_{x_1}\partial_{x_2} Q + 3 a_1 a_2^2 \partial_{x_1} \partial^2_{x_2} Q + a_2^3 \partial^3_{x_2} Q. \end{equation*} \end{exemple} %\noindent We now focus on nonlinear strategies for the reconstruction of specific variables using Leibniz formula. \subsubsection{Using Leibniz rule for specific quantities}\label{subsection:Leibniz} For any specific quantity $S$, whenever \mytextcolor{red}{$P(\mathbf{x},\rho)>0$}, we define its nonlinear reconstruction by the \emph{rational function}:% $R^{n}$: \begin{equation} \label{eq:massic_qtt_n} \mytextcolor{red}{ R(\mathbf{x},S) \coloneqq \frac{P(\mathbf{x},\rho S)}{P(\mathbf{x},\rho)}. } \end{equation} Where the reconstruction of \mytextcolor{red}{$P(\mathbf{x},\rho S)$} is given by~\eqref{eq:reconst_prod_dens_var_volumic} with $Q=\rho S$. Instead of considering the volumic variable $\rho S$, the idea is to deal with the variations of each variable $\rho$ and $S$ independently (Leibniz formula) and then gather all the information (see the Subsections~\ref{subsection:HighOrderMassic} and Remark~\ref{subsection:GatheringForAssembling} below). Assume that the fields $\rho$ and $S$ are regular, we can use the {Leibniz rule} at any differentiation order of the product $\rho S$ so that: Leibniz multi-variate approach using the identity~\eqref{eq:diff_multi_variate_dotprod}\eqref{eq:diff_multi_variate_dotprod_mindices}\eqref{eq:diff_multi_variate_dotprod_mindices_final}, replacing $Q$ by the product $\rho S$: \begin{equation} \label{eq:diff_multi_variate_dotprod_mindices_Leibniz} % (\mathbf{a}, \nabla)^n (\rho S) =\left( \sum _{i=1}^d (a_i \partial_{x_i}) \right)^n (\rho S)= \sum _{|\mathbf{m}|=n} \left( \left( \begin{array}{llll} n \\ \mathbf{m} \end{array} \right) \prod _{i=1}^{d} (a_i \partial_{x_i})^{m_i} \right) (\rho S) (\nabla)^n (\rho S) (\mathbf{a},..,\mathbf{a}) = \sum _{|\mathbf{m}|=n} \left( \begin{array}{llll} n \\ \mathbf{m} \end{array} \right) \prod _{i=1}^{d} a_i^{m_i} \partial^\mathbf{m} (\rho S). \end{equation} %where we set $|\mathbf{m}|=\sum _{i=1}^d m_i$. %\noindent Noting $\mathbf{0_d}=(0,\ldots ,0)$, we can use the multi-index Leibniz formula: \begin{equation} \label{eq:diff_multi_variate_dotprod_mindices_Leibniz_generical_multi_ind} \partial^{\mathbf{m}} (\rho S) = \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \binom{\mathbf{m}}{\mathbf{k}} \partial^{\mathbf{k}} \rho\:\: \partial^{\mathbf{m}-\mathbf{k}} S \qquad \mbox{with} \:\:\: \binom{\mathbf{m}}{\mathbf{k}} = \prod _{i=1}^d \left( \begin{array}{llll} m_i \\ k_i \end{array} \right). \end{equation} In~\eqref{eq:diff_multi_variate_dotprod_mindices_Leibniz_generical_multi_ind}, we recall some conventions for multi-indices comparison operators: \begin{definition}\label{eq:def_tilde_m} Let $\mathbf{k}$ and $\mathbf{m}$ be two multi-indices, we have the classical component-wise definition: \begin{equation*} \mathbf{k} \le \mathbf{m} \qquad %\mbox{iff} \mbox{means that} \qquad \forall i \:\: k_i \le m_i.%, \\ % \mathbf{k} < \mathbf{m} \qquad \mbox{meaning that} \qquad \forall i \:\: k_i < m_i. \end{equation*} %We define an ``intermediate'' operator: \mytextcolor{blue}{ We define the $<$ operator \begin{equation} \label{eq:compar_operator_le_neq} \mathbf{k} < \mathbf{m} \qquad %\mbox{iff} \mbox{means that} \qquad \forall i \:\: k_i \le m_i \:\:\: \mbox{and} \:\:\:\: \exists i \:\: k_i < m_i. \end{equation} } \end{definition} Now, for~\eqref{eq:diff_multi_variate_dotprod_mindices_Leibniz_generical_multi_ind}, we finally end up with \begin{equation} \label{eq:diff_multi_variate_dotprod_mindices_identif_massic_final} \forall \mathbf{a} \in \RR^d, \qquad \nabla^n (\rho S) (\mathbf{a},\ldots ,\mathbf{a}) = \sum _{|\mathbf{m}|=n} \left( \begin{array}{llll} n \\ \mathbf{m} \end{array} \right) \prod _{i=1}^{d} a_i^{m_i} \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \left( \begin{array}{llll} \mathbf{m} \\\mathbf{k} \end{array} \right) \partial^{\mathbf{k}} \rho \:\: \partial^{\mathbf{m}-\mathbf{k}} S. \end{equation} For example if $n=2$ and $d=2$, then~\eqref{eq:diff_multi_variate_dotprod_mindices_identif_massic_final} writes: \begin{enumerate} \item $ m_1=2 \:\: m_2=0$ \[ \binom{2}{(2 \: 0)} a_1^2 %(x_1 -x_{1,j})^2 \left[ \binom{2}{0} \overbrace{\partial_{x_1}^0 \rho \: \partial_{x_1}^{2-0} S}^{k_1=0} + \binom{2}{1} \overbrace{\partial_{x_1}^1 \rho \: \partial_{x_1}^{2-1} S}^{k_1=1} + \binom{2}{2} \overbrace{\partial_{x_1}^2 \rho \: \partial_{x_1}^{2-2} S}^{k_1=2} \right] \] \item $m_1=1 \:\: m_2=1$ \begin{multline*} \binom{2}{(1 \: 1)} a_1 a_ 2 %\\ %(x_1 -x_{1,j}) (x_2 -x_{2,j}) \\ \left[ \binom{1}{0}\binom{1}{0}) \overbrace{\partial_{x_1}^0\partial_{x_2}^0 \rho \: \partial_{x_1}^{1-0}\partial_{x_2}^{1-0} S}^{k_1=0,k_2=0} + \binom{1}{0}\binom{1}{1} \overbrace{\partial_{x_1}^0\partial_{x_2}^1 \rho \: \partial_{x_1}^{1-0}\partial_{x_2}^{1-1} S}^{k_1=0,k_2=1} \right.\\ \left. \hspace{2cm} + \binom{1}{1}\binom{1}{0} \overbrace{\partial_{x_1}^1\partial_{x_2}^0 \rho \: \partial_{x_1}^{1-1}\partial_{x_2}^{1-0} S}^{k_1=1,k_2=0} + \binom{1}{1}\binom{1}{1} \overbrace{\partial_{x_1}^1\partial_{x_2}^1 \rho \: \partial_{x_1}^{1-1}\partial_{x_2}^{1-1} S}^{k_1=1,k_2=1} \right] \end{multline*} \item $ m_1=0 \:\: m_2=2 $ \[ \binom{2}{(0 \: 2)} a_2^2 %(x_2 -x_{2,j})^2 \left[ \binom{2}{0}) \overbrace{\partial_{x_2}^0 \rho \: \partial_{x_2}^{2-0} S}^{k_2=0} + \binom{2}{1} \overbrace{\partial_{x_2}^1 \rho \: \partial_{x_2}^{2-1} S}^{k_2=1} + \binom{2}{2} \overbrace{\partial_{x_2}^2 \rho \: \partial_{x_2}^{2-2} S}^{k_2=2} \right] \] after summation, this gives the more common tensorial form: \begin{equation} (\nabla^2 (\rho S)\mathbf{a},\mathbf{a}) = \left( (\rho \nabla^2 S + (\nabla \rho \otimes \nabla S) + (\nabla S \otimes \nabla \rho) + S \nabla^2 \rho ) \mathbf{a},\mathbf{a}\right). \end{equation} \end{enumerate} \subsubsection{Computing arbitrary high-order specific quantities}{\label{subsection:HighOrderMassic}} In view of formula ~\eqref{eq:diff_multi_variate_dotprod_mindices_identif_massic_final} (with $k_i=0$ terms), we need to compute at least \begin{equation} \mytextcolor{red}{ \label{eq:all_high_order_massic} % \forall l=1,\ldots ,n \:\:\: \forall \mathbf{m};\:\: |\mathbf{m}|=l \qquad \prod _{i=1}^{d} \partial_{x_i}^{m_i} S, \qquad \forall l=1,\ldots ,n \:\:\: \forall \mathbf{m};\:\: |\mathbf{m}|=l } \end{equation} \mytextcolor{red}{ which is nothing but all $ {\nabla^l S}, \:\:\: l\le n$. Hence, we choose a numerical method to compute high-order terms in each cell $\Omega_j$: $({\nabla S})_j, \ldots, ({\nabla^n S})_j $ (each must be an approximation of $(\nabla^{i} S) (\mathbf{x}_j)$). Note also that by symmetry $k_i=m_i$ we need to compute all the high-order term for $\rho$: ${\nabla^l \rho}, \:\:\: l\le n $ using ~\eqref{eq:diff_multi_variate_dotprod_mindices_final}.} % \subsubsection{Gathering sub product of independent high-order variables}{\label{subsection:GatheringForAssembling}} \begin{remark}{\label{subsection:GatheringForAssembling}} We use formula~\eqref{eq:diff_multi_variate_dotprod_mindices_final} for the volumic quantity $\rho$ to obtain $\nabla^n \rho$ ($n \ge 1$) (replacing $S$ by $\rho$ in~\eqref{eq:all_high_order_massic}). For the specific quantity $S$, the high-order derivatives~\eqref{eq:all_high_order_massic} are gathered with those of $\rho$ in~\eqref{eq:diff_multi_variate_dotprod_mindices_identif_massic_final} to deduce \mytextcolor{red}{$P(\mathbf{x},\rho S)$} in~\eqref{eq:reconst_prod_dens_var_volumic}. \end{remark} \section{Obtaining homogeneity properties for non-linear reconstruction on \texorpdfstring{$\mathbf{U}$}{U} and \texorpdfstring{$\epsilon$}{epsilon}} \subsection{Arbitrary order (limitation-free)} In case of arbitrary order reconstruction of conservative quantities, we have defined in all cell $j$ \begin{align} \label{eq:rec_ro_n} \mytextcolor{red}{P(\mathbf{x},\rho)}&= \bar{\rho}_j + \sum _{l=1}^n \frac{1}{l!} (\nabla^{l} \rho)_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j), \\ \label{eq:rec_roU_n} \mytextcolor{red}{\boldsymbol{P}(\mathbf{x},\rho \mathbf{U})}&= \bar{\rho \mathbf{U}}_j + \sum _{l=1}^n \frac{1}{l!} (\nabla^{l} (\rho \mathbf{U}))_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j),\\ \label{eq:rec_roE_n} \mytextcolor{red}{P(\mathbf{x},\rho E)}&= \bar{\rho E}_j + \sum _{l=1}^n \frac{1}{l!} (\nabla^{l} (\rho E))_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j). \end{align} Using ~\eqref{eq:diff_multi_variate_dotprod_mindices_identif_massic_final}, ~\eqref{eq:rec_roU_n} and~\eqref{eq:rec_roE_n} may be written as \begin{align} \label{eq:rec_roU_n_dev} \boldsymbol{P}(\mathbf{x},\rho \mathbf{U})&= \bar{\rho \mathbf{U}}_j + \sum _{l=1}^n \frac{1}{l!} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} \mathbf{U})_j \\ \label{eq:rec_roE_n_dev} \mytextcolor{red}{P(\mathbf{x},\rho E)}&= \bar{\rho E}_j + \sum _{l=1}^n \frac{1}{l!} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} E)_j \end{align} \begin{notation} The Leibniz rule in~\eqref{eq:rec_roU_n_dev} or~\eqref{eq:rec_roE_n_dev} for the volumic variable $\rho S$ (S = $\mathbf{U}$ or E) shows a separate dependancy on both $\rho$ and the specific variable $S$. Hence, we will use the following notation \begin{equation} \mytextcolor{red}{P(\mathbf{x},\rho S)}= P(\mathbf{x},[\rho]_0^n,[S]_0^n) \ %\\ % \begin{equation} \label{eq:notation_intervall_global_high_order} \begin{minipage}[t]{70mm} where $[\,.\,]_i^{k}$ denotes partial derivatives of global order between $i$ and $k$. \end{minipage} % \end{equation} \end{equation} \end{notation} We may write now the primary specific variables $\mathbf{U}$ and $E$ as rational functions. As usual, isolating the constant value in cell from variable function, we have \begin{definition} Let the reconstructed density be given by~\eqref{eq:rec_ro_n} and assume that it does not vanish in $\Omega_j$. Then the rational reconstruction of $\mathbf{U}$ and $E$ given by~\eqref{eq:massic_qtt_n} rewrites as \begin{equation} \hspace{-.4cm} \label{eq:rec_U_n_ini1} \left\{ \begin{array}{lllllllllll} \mytextcolor{red}{\boldsymbol{R}(\mathbf{x},\mathbf{U})}&\coloneqq &\mytextcolor{red}{\frac{\boldsymbol{P}(\mathbf{x},\rho \mathbf{U})}{P(\mathbf{x},\rho) }}, \\ \mytextcolor{red}{\boldsymbol{R}(\mathbf{x},\mathbf{U})} &=& \mathbf{U}_j + \mytextcolor{red}{\boldsymbol{HR}(\mathbf{x},\mathbf{U}), \qquad \left( \boldsymbol{HR}(\mathbf{x},\mathbf{U})= \frac{\boldsymbol{P}(\mathbf{x},\rho \mathbf{U})-\mathbf{U}_j P(\mathbf{x},\rho) }{P(\mathbf{x},\rho)}=\frac{\boldsymbol{HP}(\mathbf{x},\rho \mathbf{U})}{P(\mathbf{x},\rho)}\right)}, \end{array} \right. \end{equation} and \begin{equation} \hspace{-.4cm} \label{eq:rec_E_n_ini1} \left\{ \begin{array}{lllllllllll} \mytextcolor{red}{R(\mathbf{x},E)}&\coloneqq &\mytextcolor{red}{ \frac{P(\mathbf{x},\rho E)}{P(\mathbf{x},\rho) }}, \\ \mytextcolor{red}{R(\mathbf{x},E)} &=& E_j + \mytextcolor{red}{HR(\mathbf{x},E), \qquad \left( HR(\mathbf{x},E)= \frac{P(\mathbf{x},\rho E)- E_j P(\mathbf{x},\rho) }{P(\mathbf{x},\rho)}=\frac{HP(\mathbf{x},\rho E)}{P(\mathbf{x},\rho)}\right)}, \end{array} \right. \end{equation} where each of the \mytextcolor{red}{$HR(\mathbf{x}, S)$} terms in~\eqref{eq:rec_U_n_ini1} and~\eqref{eq:rec_E_n_ini1} ($S=E$ or $S=\mathbf{U}$) may be considered as a rational high-order correction with respect to the constant term. \end{definition} \begin{definition}{\label{def_internal_nrj_reconstruction}Specific internal energy reconstruction}\\ % Now, we give the reconstruction associated The specific internal energy reconstruction is defined as: %of the specific (secondary) internal energy $\epsilon$ is defined by using the reconstruction of the specific (primary) variables $\mathbf{U}$~\eqref{eq:rec_U_n_ini1} and $E$~\eqref{eq:rec_E_n_ini1}: \begin{equation} \label{eq:rel_reconstruction_epsEU2} \mytextcolor{red}{R(\mathbf{x},\epsilon)}\coloneqq \mytextcolor{red}{R(\mathbf{x},E) - \frac{1}{2} |\boldsymbol{R}(\mathbf{x},\mathbf{U})|^2.} \end{equation} It is then obtained by using the reconstruction of the -primary- specific variables $\mathbf{U}$~\eqref{eq:rec_U_n_ini1} and $E$~\eqref{eq:rec_E_n_ini1}, and it may be therefore considered as a -secondary- specific variable. \end{definition} %\newpage \begin{lemma} The high-order part \mytextcolor{red}{$HR(\mathbf{x},\cdot\, )$} in~\eqref{eq:rec_U_n_ini1} and~\eqref{eq:rec_E_n_ini1} has the following properties: \begin{enumerate} \item it reads as a rational function (see notation in~\eqref{eq:notation_intervall_global_high_order}): \begin{equation} \label{eq:rec_S_form_prop} \mytextcolor{red}{HR(\mathbf{x},S)\coloneqq \frac{HP(\mathbf{x},[\rho]_0^{n-1},[S]_1^{n})}{P(\mathbf{x},\rho)}} \end{equation} % \begin{equation} % \label{eq:notation_intervall_global_high_order} %\mbox{where} \:\: [q]_i^{k} \:\: \mbox{denotes partial derivatives of global order between $i$ and $k$}. % \end{equation} \item the numerator in~\eqref{eq:rec_S_form_prop} is linear with respect to all $S$ derivatives: \begin{equation} \label{eq:homogeneity_1} \mytextcolor{red}{\forall \lambda \in \RR \qquad HP(\mathbf{x},[\rho]_0^{n-1}, [\lambda S]_1^{n}) = \lambda HP(\mathbf{x},[\rho]_0^{n-1},[S]_1^{n}).} \end{equation} \end{enumerate} \end{lemma} \begin{proof} First, using~\eqref{eq:cte_correct_highorder}, \begin{align*} c^l_j(\rho S) &= \frac{1}{l!} \frac{1}{|\Omega_j|}\int_{\Omega_j} (\nabla^l (\rho S))_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j) \dd \mathbf{x}, \\ &= \frac{1}{l!}\frac{1}{|\Omega_j|}\int_{\Omega_j} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} S)_j\, \dd \mathbf{x} \end{align*} We recall that: \begin{equation} \label{eq:dec_leibniz_multi_prod} \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho) \:\: (\partial^{\mathbf{m}-\mathbf{k}} S) = \sum _{0 \le k_1 \le m_1}\cdots \sum _{0 \le k_d \le m_d} \prod _{i} \binom{m_i}{k_i} \left(\prod _{i} \partial_{x_i}^{k_i}\right) \rho \left(\prod _{i} \partial_{x_i}^{m_i-k_i}\right) S. \end{equation} In order to clarify manipulation of multi-index sums, we use Definition~\ref{eq:def_tilde_m}~\eqref{eq:compar_operator_le_neq}. Now, isolating the constant term for $S$ ($k_i=m_i$, $\forall i$) in~\eqref{eq:dec_leibniz_multi_prod}, we have: \begin{equation} \label{eq:isolating_leibniz_cteS} \sum _{\mathbf{0_d} \le \mathbf{k} \le \mathbf{m}} \binom{\mathbf{m}}{\mathbf{k}} \partial^{\mathbf{k}} \rho \:\: \partial^{\mathbf{m}-\mathbf{k}} S = S \partial^{\mathbf{m}} \rho + \mytextcolor{blue}{\sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}}\binom{\mathbf{m}}{\mathbf{k}} \partial^{\mathbf{k}} \rho \:\: \partial^{\mathbf{m}-\mathbf{k}} S, \end{equation} hence \begin{equation*} c^l_j(\rho S)= S_j c^l_j(\rho) + \delta^l_j(\rho S), \end{equation*} where we have noted \begin{equation} \label{eq:CtR} \delta^l_j(\rho S)\coloneqq \frac{1}{l!} \frac{1}{|\Omega_j|}\int_{\Omega_j} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \mytextcolor{blue}{ \sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}}\binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} S)_j\, \dd \mathbf{x}. \end{equation} We thus have \begin{equation} \label{eq:rho_S_isol} \bar{\rho S}_j = S_j \bar{\rho}_j - \sum _{l=1}^n \delta_j^{l} (\rho S). \end{equation} Since \begin{equation*} \mytextcolor{red}{P(\mathbf{x},\rho S)} = \bar{\rho S}_j + \sum _{l=1}^n \frac{1}{l!} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \mytextcolor{blue}{\sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} S)_j, \end{equation*} we can isolate the constant term for S (see~\eqref{eq:isolating_leibniz_cteS}). Thus, using~\eqref{eq:rho_S_isol}, we have \begin{align} \nonumber \mytextcolor{red}{P(\mathbf{x},\rho S)}{}={}& S_j \left(\bar{\rho}_j + \sum _{l=1}^n \frac{1}{l!} (\nabla^{l} \rho)_j (\mathbf{x}-\mathbf{x}_j,..,\mathbf{x}-\mathbf{x}_j) \right) - \sum _{l=1}^n \delta_j^{l} (\rho S) \\ \nonumber & \hspace*{25mm}+ \sum _{l=1}^n \frac{1}{l!} \sum _{|m|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \mytextcolor{blue}{\sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} S)_j \\ \label{eq:rho_S_isol2} \mytextcolor{red}{P(\mathbf{x},\rho S)} {}\coloneqq{}& \mytextcolor{red}{ S_j P(\mathbf{x},\rho) + HP(\mathbf{x},[\rho]_0^{n-1},[S]_1^{n})}, \end{align} where (see also~\eqref{eq:CtR}): \begin{equation}\label{eq:Hrest_massic_var_leibniz} HP(\mathbf{x},[\rho]_0^{n-1},[S]_1^{n}) \coloneqq - \sum _{l=1}^n \delta_j^{l} (\rho S) + \sum _{l=1}^n \frac{1}{l!} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \mytextcolor{blue}{\sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} S)_j \qedhere \end{equation} \end{proof} \mytextcolor{red}{ Notice that by keeping explicit dependency of different energies $\epsilon$ and $\frac{1}{2} |\mathbf{U}|^2$ in total specific energy $E$, with $E_j=\epsilon_j+\frac{1}{2} |\mathbf{U}_j|^2$, we have:} \begin{align} \nonumber \mytextcolor{red}{HR(\mathbf{x},E)}{}={}&\mytextcolor{red}{ \frac{P(\mathbf{x},\rho \epsilon)+P(\mathbf{x},\rho \frac{1}{2} | \mathbf{U}|^2) - (\epsilon_j + \frac{1}{2} | \mathbf{U}_j|^2) P(\mathbf{x},\rho) }{P(\mathbf{x},\rho)}}\\ \nonumber {}={}&\mytextcolor{red}{\frac{\big(P(\mathbf{x},\rho \epsilon) - \epsilon_j P(\mathbf{x},\rho)\big) + \big(P(\mathbf{x},\rho \frac{1}{2} | \mathbf{U}|^2) - \frac{1}{2} | \mathbf{U}_j|^2 P(\mathbf{x},\rho)\big) }{P(\mathbf{x},\rho)}}\\ \label{eq:RelE_eps_cin} \mytextcolor{red}{HR(\mathbf{x},E)} {}\coloneqq{}& \mytextcolor{red}{HR^{(E)}(\mathbf{x},\epsilon)+ HR^{(E)}\left(\mathbf{x},\frac{1}{2}|\mathbf{U}|^2\right)}. \end{align} \mytextcolor{red}{ this may define another way of computing the reconstructed specific total energy: \begin{align} \label{eq:tot_spec_E_developped} R^{\sumrm}(\mathbf{x},E){}\coloneqq {}&R^{(E)}(\mathbf{x},\epsilon) + R^{(E)}\left(\mathbf{x},\frac{1}{2} | \mathbf{U}|^2\right),\\ \nonumber {}={}& \epsilon_j + HR^{(E)}(\mathbf{x},\epsilon) + \frac{1}{2} |\mathbf{U}_j|^2 + HR^{(E)}\left(\mathbf{x},\frac{1}{2}|\mathbf{U}|^2\right)= E_j + HR^{(E)}(\mathbf{x},\epsilon) + HR^{(E)}\left(\mathbf{x},\frac{1}{2}|\mathbf{U}|^2\right) \end{align} As we will see thereafter, the core of the approach is to express the high-order part of the kinetic energy in relation~\eqref{eq:RelE_eps_cin} as a function of high order terms of velocity already computed in \mytextcolor{red}{$\boldsymbol{R}(\mathbf{x},\mathbf{U})$}. Indeed, the variable $\epsilon$ needs a specific attention: its definition makes use of the kinetic energy, so that the interplay between $\epsilon$ and $\mathbf{U}$ is not trivial. } \subsection{Computing high-order kinetic energy terms with multi-variate Leibniz formula} \mytextcolor{red}{ We recall the two specific total energy reconstruction~\eqref{eq:rel_reconstruction_epsEU2} and~\eqref{eq:tot_spec_E_developped}: \begin{equation} \label{eq:equal_two_rel_totE} \begin{cases} R^{\defrm}(\mathbf{x},E)=R(\mathbf{x},\epsilon) + \frac{1}{2} |R(\mathbf{x},\mathbf{U})|^2, \\ R^{\sumrm}(\mathbf{x},E)=R^{(E)}(\mathbf{x},\epsilon) + R^{(E)}(\mathbf{x},\frac{1}{2} | \mathbf{U}|^2). \end{cases} \end{equation} Note that imposing $R^{\defrm}(\mathbf{x},E)=R^{\sumrm}(\mathbf{x},E)$ in~\eqref{eq:equal_two_rel_totE} gives an equation for $R(\mathbf{x},\epsilon)$. Hence, owing to~\eqref{eq:rec_E_n_ini1},~\eqref{eq:rec_U_n_ini1}, the relation~\eqref{eq:rel_reconstruction_epsEU2} read \begin{align*} \mytextcolor{red}{R(\mathbf{x},\epsilon)}&=\mytextcolor{red}{ E_j + HR(\mathbf{x},E) -\frac{1}{2} |\mathbf{U}_j+\boldsymbol{HR}(\mathbf{x},\mathbf{U})|^2}, \\ &= \mytextcolor{red}{E_j + HR(\mathbf{x},E) -\frac{1}{2} \left( |\mathbf{U}_j|^2 + 2 (\boldsymbol{HR}(\mathbf{x},\mathbf{U}),\mathbf{U}_j) + |\boldsymbol{HR}(\mathbf{x},\mathbf{U})|^2 \right)}, \\ &= \mytextcolor{red}{E_j - \frac{1}{2} |\mathbf{U}_j|^2 + HR(\mathbf{x},E) - (\boldsymbol{HR}(\mathbf{x},\mathbf{U}),\mathbf{U}_j) -\frac{1}{2} |\boldsymbol{HR}(\mathbf{x},\mathbf{U})|^2}. \end{align*} } Using~\eqref{eq:RelE_eps_cin}, we obtain: \begin{align} \label{eq:rec_epsilon_two_components_homogenous} \mytextcolor{red}{R(\mathbf{x},\epsilon)} =\mytextcolor{red}{\epsilon_j + HR^{(E)}(\mathbf{x},\epsilon) + \underbrace{\left(HR^{(E)}\left(\mathbf{x},\frac{1}{2}|\mathbf{U}|^2\right) - (\boldsymbol{HR}(\mathbf{x},\mathbf{U}),\mathbf{U}_j) \right)}_{\coloneqq D(\mathbf{x},\mathbf{U})} - \frac{1}{2} |\boldsymbol{HR}(\mathbf{x},\mathbf{U})|^2}. \end{align} \begin{proposition} The scalar valued rational function \mytextcolor{red}{$D(\mathbf{x},\mathbf{U})$} in~\eqref{eq:rec_epsilon_two_components_homogenous} has the following properties: \begin{enumerate} \item it writes (see notation in~\eqref{eq:notation_intervall_global_high_order}) \begin{equation} \label{eq:rec_DeltaUkin_form_prop} \mytextcolor{red}{D(\mathbf{x},\mathbf{U})\coloneqq \frac{D(\mathbf{x},[\rho]_0^{n-1},[\mathbf{U}]_1^{n-1})}{P(\mathbf{x},\rho)}} \end{equation} \item the numerator of~\eqref{eq:rec_DeltaUkin_form_prop} is a second-order homogeneous polynomial w.r.t the last variable: \begin{equation} \label{eq:homogeneity_2} \mytextcolor{red}{ \mbox{{\em order 2 homogeneous:}}\:\: \forall \lambda \in \RR \: \qquad D(\mathbf{x},[\rho]_0^{n-1},[\lambda \mathbf{U}]_1^{n-1}) = \lambda^2 D(\mathbf{x},[\rho]_0^{n-1},[\mathbf{U}]_1^{n-1})} \end{equation} \end{enumerate} \end{proposition} \begin{proof} By~\eqref{eq:RelE_eps_cin} and~\eqref{eq:rec_S_form_prop},~\eqref{eq:rec_U_n_ini1} (same denominator), \mytextcolor{red}{$D(\mathbf{x},\mathbf{U})$} reads as a rational function \begin{align} \mytextcolor{red}{D(\mathbf{x},\mathbf{U})} &=\mytextcolor{red}{ \frac{HP(\mathbf{x},\rho\frac{1}{2}|\mathbf{U}|^2) - (\boldsymbol{HP}(\mathbf{x},\rho\mathbf{U}),\mathbf{U}_j)}{P(\mathbf{x},\rho)}} \\ &=\mytextcolor{red}{\frac{HP(\mathbf{x},[\rho]_0^{n-1},[\frac{1}{2} |\mathbf{U}|^2]_1^{n}) - \left(\boldsymbol{HP}(\mathbf{x},[\rho]_0^{n-1},[\mathbf{U}]_1^{n}),\mathbf{U}_j\right)}{P(\mathbf{x},\rho)}} \label{eq:DjU} \end{align} Using~\eqref{eq:Hrest_massic_var_leibniz} for the kinetic energy ($S=\frac{1}{2}|\mathbf{U}|^2$), the first term of the numerator in~\eqref{eq:DjU} reads as \begin{multline} \label{eq:Hrest_massic_var_leibniz_ke} \mytextcolor{red}{ HP\left(\mathbf{x},[\rho]_0^{n-1}, \left[ \frac{1}{2} |\mathbf{U}|^2 \right] _1^{n}\right)} \\ = \left(- \sum _{l=1}^n \delta_j^{l} \left(\rho \frac{1}{2} |\mathbf{U}|^2\right) + \sum _{l=1}^n \frac{1}{l!} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \mytextcolor{blue}{\sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} \left(\frac{1}{2} |\mathbf{U}|^2\right)_j\right), \end{multline} and the second term reads as: \begin{multline} \label{eq:Hrest_massic_var_leibniz_velocity_dot_Cte} (\mytextcolor{red}{\boldsymbol{HP}(\mathbf{x},[\rho]_0^{n-1}, [\mathbf{U} ] _1^{n})}, \mathbf{U}_j) \\ = \left(- \sum _{l=1}^n \delta_j^{l} (\rho \mathbf{U}) + \sum _{l=1}^n \frac{1}{l!} \sum _{|\mathbf{m}|=l} \binom{l}{\mathbf{m}} \prod _{i=1}^{d} (x_i-x_{i,j})^{m_i} \mytextcolor{blue}{\sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}}} \binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \rho)_j \:\: (\partial^{\mathbf{m}-\mathbf{k}} \mathbf{U})_j , \mathbf{U}_j \right). \end{multline} Now, subtracting~\eqref{eq:Hrest_massic_var_leibniz_velocity_dot_Cte} to~\eqref{eq:Hrest_massic_var_leibniz_ke}, we must give an expression of \begin{equation} \label{eq:diff_analysis} \left(\partial^{\mathbf{m}-\mathbf{k}} \left(\frac{1}{2} |\mathbf{U}|^2\right)\right)_j - ((\partial^{\mathbf{m}-\mathbf{k}} \mathbf{U})_j , \mathbf{U}_j) \end{equation} The second term in~\eqref{eq:diff_analysis} writes (omitting cell index $j$): \begin{equation} \label{eq:diff_analysis_sec} (\partial^{\mathbf{m}-\mathbf{k}} \mathbf{U} , \mathbf{U})= \sum _{s=1}^d u_s \partial^{\mathbf{m}-\mathbf{k}} u_s. \end{equation} and the first reads as: \begin{align} \nonumber \partial^{\mathbf{m}-\mathbf{k}} (\frac{1}{2} |\mathbf{U}|^2) &= \frac{1}{2} \sum _{s=1}^d \partial^{\mathbf{m}-\mathbf{k}} (u_s u_s), \\ &= \frac{1}{2} \sum _{s=1}^d \:\ \sum _{\mathbf{0_d} \le \mathbf{t} \le \mathbf{m}-\mathbf{k}} \binom{ \mathbf{m}-\mathbf{k}}{\mathbf{t}} \partial^{\mathbf{t}} u_s \: \partial^{\mathbf{m}-\mathbf{k}-\mathbf{t}} u_s. \label{eq:diff_analysis_pre} \end{align} Isolating the constant term with respect to $u_s$, ie considering the two extremal terms in the sum over $\mathbf{t}$, that is: \begin{equation*} \forall i \:\: t_i=0 \:\: \mbox{ and } \:\: \forall i \:\: t_i=m_i-k_i, \end{equation*} and using the previous Definition~\ref{eq:def_tilde_m} (see~\eqref{eq:compar_operator_le_neq}) also for the two bounds $\mathbf{t}=\mathbf{0}_d$ and $\mathbf{t}=\mathbf{m}-\mathbf{k}$, %the sum ``lower bound'', we have: \begin{equation} \label{eq:isolt_two_external_terms} \sum _{\mathbf{0_d} \le \mathbf{t} \le \mathbf{m}-\mathbf{k}} \binom{\mathbf{m}-\mathbf{k}}{\mathbf{t}} \partial^{\mathbf{t}} u_s \:\: \partial^{\mathbf{m}-\mathbf{k}-\mathbf{t}} u_s= (u_s) \: \partial^{\mathbf{m}-\mathbf{k}} u_s + \partial^{\mathbf{m}-\mathbf{k}} u_s \: (u_s) + \mytextcolor{blue}{\sum _{\mathbf{0_d} < \mathbf{t} < \mathbf{m}-\mathbf{k}}} \binom{\mathbf{m}-\mathbf{k}}{\mathbf{t}} \partial^{\mathbf{t}} u_s \:\: \partial^{\mathbf{m}-\mathbf{k}-\mathbf{t}} u_s \end{equation} so that using~\eqref{eq:isolt_two_external_terms},~\eqref{eq:diff_analysis_pre} and~\eqref{eq:diff_analysis_sec} in~\eqref{eq:diff_analysis}, we finally end up with~\eqref{eq:rec_DeltaUkin_form_prop}: \begin{equation} \label{eq:expressionform_homogenous2} \mytextcolor{red}{D(\mathbf{x},[\rho]_0^{n-1},[\mathbf{U}]_1^{n-1}) = %-D_j+ HP\left(\mathbf{x},[\rho]_0^{n-1},\left[\frac{1}{2}\big|\mathbf{U}\big|^2\right]_1^{n-1} \right)}, \end{equation} with \begin{multline} \hspace*{-3mm} HP\left(\mathbf{x},[\rho]_0^{n-1},\left[\frac{1}{2}\big|\mathbf{U}\big|^2\right]_1^{n-1} \right) = - \sum _{l=1}^n \left( \delta_j^{l} \left(\rho \left(\frac{1}{2} |\mathbf{U}|^2\right)\right) - (\delta_j^{l} (\rho \mathbf{U}),\mathbf{U}_j) \right) \\ + \Mk \sum _{l\mk =\mk 1}^n \frac{1}{l!} \sum _{|\mathbf{m}|\mk =\mk l} \mk \binom{l}{\mathbf{m}}\Mk \prod _{i\Mk =\Mk 1}^{d}\mk (x_i\Mk -\Mk x_{i,j})^{m_i} \hspace{-.15cm} \sum _{\mathbf{0_d} \le \mathbf{k} < \mathbf{m}} \Mk\binom{\mathbf{m}}{\mathbf{k}} (\partial^{\mathbf{k}} \mk \rho)_j \Mk \left(\frac{1}{2} \sum _{s=1}^d \, \sum _{\mathbf{0_d} < \mathbf{t} < \mathbf{m}-\mathbf{k}} \binom{\mathbf{m}-\mathbf{k}}{\mathbf{t}} (\partial^{\mathbf{t}} u_s)_j (\partial^{\mathbf{m}-\mathbf{k}-\mathbf{t}} u_s)_j\Mk \right)\Mk . \label{eq:finalform_homogenous2} \end{multline} In the last internal sum in~\eqref{eq:finalform_homogenous2}, due to the strict inequality bounds \mytextcolor{blue}{$\mathbf{0_d} < \mathbf{t} < \mathbf{m}-\mathbf{k}$} %( only two cases are possible %leading to the two events (see~\eqref{eq:compar_operator_le_neq}), we have: \begin{itemize} \item Lower bound: $\exists i^-$ such that $t_{i^-}>0$ and \item Upper bound: $\exists i^+$ such that $(m_{i^+}-k_{i^+}-t_{i^+})>0$, \end{itemize} %so that $\forall \boldsymbol{\alpha}; \mathbf{0_d} < \boldsymbol{\alpha}$, meaning that $\exists i; \alpha_i>0$ considering now $u_s \to \lambda u_s$, is nothing but multiplying $\partial^{\mathbf{t}} u_s$ and $\partial^{\mathbf{m}-\mathbf{k}-\mathbf{t}} u_s $ by $\lambda$, which gives the expected 2-homogeneity property. \end{proof} Hence, we have obtain in~\eqref{eq:expressionform_homogenous2} and~\eqref{eq:finalform_homogenous2} high-order terms for kinetic energy as a result of those of velocity (already computed!), note however that the highest order is $n-1$ (and not $n$). \subsection{A limited reconstruction of internal energy} With the previous results~\eqref{eq:rec_epsilon_two_components_homogenous}~\eqref{eq:expressionform_homogenous2}, we can establish the following unlimited reconstruction of specific internal energy (here, let just suppose that $P(\mathbf{x},\rho) \ne 0$): \begin{align} \label{eq:rec_eps_unlim} R(\mathbf{x},\epsilon) = \epsilon_j + \frac{HP(\mathbf{x},[\rho]_0^{n-1},[\mathbf{\epsilon}]_1^n)} {P(\mathbf{x},\rho)} +\frac{HP\left(\mathbf{x},[\rho]_0^{n-1},\left[\frac{1}{2}|\mathbf{U}|^2\right]_1^{n-1}\right)} {P(\mathbf{x},\rho)} -\frac{1}{2} \Bigg|\frac{\boldsymbol{HP}(\mathbf{x},[\rho]_0^{n-1},[\mathbf{U}]_1^n)}{P(\mathbf{x},\rho)} \Bigg|^2 \end{align} The high-order part of specific internal energy reconstruction is build with three high-order contributions. The first comes from $\epsilon$ data itself, %~\eqref{eq:rec_eps_lim1}, the second term %in~\eqref{eq:rec_eps_lim2} is made from kinetic energy data all derivatives are obtain by those coming from velocity thanks to Leibniz formula (see~\eqref{eq:expressionform_homogenous2}~\eqref{eq:finalform_homogenous2}). Finally, the third term %in~\eqref{eq:rec_eps_lim3} comes from the high-order velocity reconstruction and the algebraic relation~\eqref{eq:rel_reconstruction_epsEU2}. \\ Let $\lambda_{\rho}$, $\lambda_{\epsilon}$, $\lambda_{\mathbf{U}}$ be function parameters (called limiters) for $\rho$, $\epsilon$ and $\mathbf{U}$. We begin by introducing the limited version of the reconstructed volumic density (at least verifying positivity at quadrature points, cf~\eqref{eq:limi_hierarch_rho} or~\eqref{eq:limi_globally_rho}): \begin{equation} \label{eq:rec_rho_lim} \mytextcolor{red}{ P^{\lambda}(\mathbf{x},\rho) \:\:\:\: \mbox{and \: noting} \:\: [\rho^{\lambda}]_i^{k} \:\: \mbox{the resulting limited version of}\:\: [\rho]_i^{k} } \end{equation} \mytextcolor{red}{ Now, collecting all the properties linked to high-order part of all the specific variable of the previous subsections, we can rewrite a function parameterized version of %reconstructed velocity $\mathbf{U}$ and specific internal energy $\epsilon$ of~\eqref{eq:rec_epsilon_two_components_homogenous} (see also~\eqref{eq:rec_U_n_ini1}\eqref{eq:rec_E_n_ini1}\eqref{eq:rel_reconstruction_epsEU2}\eqref{eq:rec_S_form_prop}): \begin%{equation} {subequations}\label{eq:rec_eps_lim} \begin{align} \label{eq:rec_eps_lim1} \hspace{-1cm} R^{\lambda}(\mathbf{x},\epsilon) = \epsilon_j &+ \lambda_{\epsilon}\frac{HP(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{\epsilon}]_1^n)} {P^{\lambda}(\mathbf{x},\rho)}%\label{eq:rec_eps_lim1} \\&+(\lambda_{\mathbf{U}})^2 \frac{HP\left(\mathbf{x},[\rho^{\lambda}]_0^{n-1},\left[\frac{1}{2}|\mathbf{U}|^2\right]_1^{n-1}\right)} {P^{\lambda}(\mathbf{x},\rho)}\label{eq:rec_eps_lim2} \\ &-\frac{(\lambda_{\mathbf{U}})^2}{2} \Bigg|\frac{\boldsymbol{HP}(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{U}]_1^n)}{P^{\lambda}(\mathbf{x},\rho)} \Bigg|^2 \label{eq:rec_eps_lim3} \end{align} \end{subequations} } \mytextcolor{red}{The high-order part in~\eqref{eq:rec_eps_lim1} %(specific internal energy) and~\eqref{eq:rec_eps_lim3} %(velocity) use linear dependancy~\eqref{eq:rec_S_form_prop}~\eqref{eq:homogeneity_1} and the one in ~\eqref{eq:rec_eps_lim2} %(specific kinetic energy) is derived from the 2-homogeneous property~\eqref{eq:rec_DeltaUkin_form_prop}~\eqref{eq:homogeneity_2}~\eqref{eq:expressionform_homogenous2}. } % (notice that in each case $n \ge 1$). \section{Construction of induced limitation for velocity reconstruction}\label{section:induced_limitation} In view of physical constraints for the continuous system $\rho \in I_{\density}$, $\epsilon \in I_{\energy}$, we may decide to apply a \emph{direct} limitation only for $\rho$, $\epsilon$ for arbitrary order reconstruction (not only second or third order like in~\cite{PH_NodalSolv:hal-03585115}). %We give thereafter some assumptions: %\begin{itemize} %\item Assume that any conservative quantity $\boldsymbol \U = (\rho, \rho \mathbf{U}, \rho E)$ is approximated by polynomials of degree $n$ : % $P(\mathbf{x},\boldsymbol \U)$.\\ %\item Assume there exists two scalar valued polynomials of degree $n$, noted $\tilde{P}$ (containing high order terms of: internal energy $\epsilon$ and velocity $\mathbf{U}$) and $\hat{P}$ (containing high order terms of kinetic energy $\frac{1}{2}|\mathbf{U}|^2$) such that: %\begin{equation*} % \begin{array}{lllllllllll} % \tilde{P}(\mathbf{x},[\rho]_0^{n-1}, [\lambda f_i]_1^n)=\lambda \tilde{P}(\mathbf{x},[\rho]_0^{n-1},[f_i]_1^n), \\ % \hat{P}(\mathbf{x},[\rho]_0^{n-1}, [\lambda\mathbf{U}]_1^{n-1})= \lambda^2 \hat{P}(\mathbf{x},[\rho]_0^{n-1}, [\mathbf{U}]_1^{n-1}).\\ % \end{array} %\end{equation*} %\end{itemize} %\subsection{A limited reconstruction of internal energy and consequence} %Now, we can write a limited version of reconstructed specific internal energy $\epsilon$ of~\eqref{eq:rec_epsilon_two_components_homogenous}: % \begin{equation} % \label{eq:rec_eps_lim} % \hspace{-1cm} % R^{\lambda}(\mathbf{x},\epsilon) = \epsilon_j + \lambda_{\epsilon}\frac{\tilde{P}(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{\epsilon}]_1^n)} % {P^{\lambda}(\mathbf{x},\rho)} % + (\lambda_{\mathbf{U}})^2 \frac{\hat{P}(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{U}]_1^{n-1})} % {P^{\lambda}(\mathbf{x},\rho)} % - \frac{(\lambda_{\mathbf{U}})^2}{2} % \Bigg|\frac{\tilde{\mathbf{P}}(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{U}]_1^n)}{P^{\lambda}(\mathbf{x},\rho)} \Bigg|^2 % \end{equation} \mytextcolor{red}{We could limit independently $\lambda_{\epsilon}$ and $\lambda_{\mathbf{U}}$ (both beeing non-negative scalars), but it is interesting to notice that choosing in~\eqref{eq:rec_eps_lim1}{\rm --}\eqref{eq:rec_eps_lim3},} \begin{equation} \label{link_epsU_lim} % \lambda_{\epsilon} \ge 0 \lambda_{\mathbf{U}}=\sqrt{\lambda_{\epsilon}} % \qquad \mbox{(with \:} \lambda_{\epsilon} \ge 0 \mbox{)} \end{equation} \mytextcolor{red}{ we obtain an induced limiter for the velocity (built from that of internal specific energy).} %In view of~\eqref{eq:rec_eps_lim}-\eqref{eq:rec_eps_lim3} and~\eqref{link_epsU_lim}, we establish the following result. % \begin{lemma} \mytextcolor{blue}{ Using~\eqref{eq:rec_eps_lim1}{\rm --}\eqref{eq:rec_eps_lim3} and~\eqref{link_epsU_lim}, we obtain a \emph{direct} limited reconstruction for the internal energy $\epsilon$:} \begin{multline} \label{eq:lim_eps} R^{\lambda}(\mathbf{x},\epsilon) \\ = \Mk \epsilon_j \Mk + \Mk \lambda_{\epsilon} \left( \frac{HP(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{\epsilon}]_1^n)} {P^{\lambda}(\mathbf{x},\rho)} \Mk +\Mk \frac{HP(\mathbf{x},[\rho^{\lambda}]_0^{n-1},\left[\frac{1}{2}|\mathbf{U}|^2\right]_1^{n-1})} {P^{\lambda}(\mathbf{x},\rho)} \Mk -\Mk \frac{1}{2}\Mk \left|\frac{\boldsymbol{HP}(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{U}]_1^n)}{P^{\lambda}(\mathbf{x},\rho)} \right|^2 \right) \end{multline} and an \emph{induced} limited velocity reconstruction for $\mathbf{U}$ using~\eqref{eq:rec_U_n_ini1}~\eqref{eq:homogeneity_1} and~\eqref{link_epsU_lim}: \begin{equation} \label{eq:lim_velocity} \mytextcolor{red}{ \boldsymbol{R}^{\lambda}(\mathbf{x},\mathbf U)}=\mytextcolor{red}{ \mathbf{U}_j + \lambda_{\mathbf{U}} \frac{\boldsymbol{HP}(\mathbf{x},[\rho^{\lambda}]_0^{n-1},[\mathbf{U}]_1^n)}{ P^{\lambda}(\mathbf{x},\rho)}}, \qquad \mbox{with} \:\: \lambda_{\mathbf{U}}=\sqrt{\lambda_{\epsilon}}. \end{equation} \mytextcolor{red}{ Now, we have obtained what we claimed at the begining: the velocity reconstruction~\eqref{eq:lim_velocity} is limited thanks to the limitation of density and the specific internal energy reconstruction~\eqref{eq:lim_eps} which are here, the variables of interest. Note also, that in view of~\eqref{eq:lim_eps} and~\eqref{link_epsU_lim}, the reconstructed total specific energy $E$ have a same limitation factor from his specific internal and kinetic energies. This result seems to be optimal in terms of limiter for $E$ which balances the limitation of the two sub energies $\epsilon$ and $\frac{1}{2} |\mathbf{U}|^2$.} %\begin{proof} %For~\eqref{eq:lim_eps}, apply the relation~\eqref{link_epsU_lim} in~\eqref{eq:rec_eps_lim}. For~\eqref{eq:lim_velocity}, just use~\eqref{eq:rec_U_n_ini1}. %\end{proof} %\end{lemma} % \begin{remarque} \: \paragraph{Consequences} \begin{itemize} \item Each of~\eqref{eq:lim_eps} and~\eqref{eq:lim_velocity} may be written as \mytextcolor{red}{$R^{\lambda}(\mathbf{x},Q)=Q_j+ \lambda_j \hat{HR}(\mathbf{x},Q )$}, where \mytextcolor{red}{$\hat{HR}(\mathbf{x},Q)$} is the ``high-order'' extension (polynomial for volumic variables and rational fraction for specific quantities). Taking $\lambda_j\equiv 1$ gives the unlimited version (and for conservative variables~\eqref{eq:rec_roU_n_dev}~\eqref{eq:rec_roE_n_dev}) while $\lambda_j \equiv 0$ reduces to first order $\mytextcolor{red}{R^{\lambda}(\mathbf{x},Q)}=Q_j$ (and then $\mytextcolor{red}{\boldsymbol{P}(\mathbf{x},\rho \mathbf{U})}= \rho_j \mathbf{U}_j$ and $\mytextcolor{red}{P(\mathbf{x},\rho E)}=\rho_j E_j$). \item Classical scalar limiters can be used, Dukowicz, Barth--Jespersen~\cite{BarJesper89}, MLP~\cite{MLPlimiting} or min-mod. They all give $\mytextcolor{red}{R^{\lambda}(\mathbf{x}_l,Q)} \ge 0$ for all quadrature points $\mathbf{x}_l$ of $\partial \Omega_j$. Some of them also impose that \begin{equation} \label{eq:TVB} \min _{j' \in V(j)} Q _{j'} \le \mytextcolor{red}{R^{\lambda}(\mathbf{x}_l,Q)} \le \max _{j' \in V(j)} Q _{j'} \end{equation} where $V(j)$ is the set of all neighbooring cells around $j$ (those having at least a common node with $j$), see Figure~\ref{Fig: V_j}. \begin{figure}[htb]\centering \vspace*{-2mm} \includegraphics[height=7cm]{Voisinage.pdf} \vspace*{-1mm} \caption{\label{Fig: V_j}Neighborood of generic cell $j$: $V_j$ a cell sharing a node or an edge with cell $j$, quadrature boundary points $\{x_r\}_{r \in \partial \Omega_j}$,$\{x_{r+1/2}\}_{r+1/2 \in \partial \Omega_j}$.} \vspace*{-2mm} \end{figure} \item If $\lambda_{\epsilon}=0$ then $\lambda_{\mathbf{U}}=0$ (first order taking also $\lambda_{\rho}=0$). If $\lambda_{\epsilon}=1$ then $\lambda_{\mathbf{U}}=1$ (unlimited scheme taking also $\lambda_{\rho}=1$). Now, in case of $\lambda_{\epsilon} \in ]0,1[$ then $\lambda_{\mathbf{U}} > \lambda_{\epsilon}$ (see Figure~\ref{fig:sqrt_eps}). \goodbreak \item we emphasize that in our limitation process, we need to define \emph{only one limiter} for both the internal energy $\epsilon$ in a direct way and all the component of the velocity $\mathbf{U}$ through~\eqref{link_epsU_lim}. \begin{figure}[htb]%!htb] \centering \includegraphics[scale=.45]{fig_compar_sqrteps.pdf} \caption{\label{fig:sqrt_eps} Behavior of the induced limiter of velocity $\lambda_{\mathbf{U}}$ with respect to limiter of internal specific energy $\lambda_{\epsilon}$ on [0,1] } \end{figure} \end{itemize} Note that the reconstruction of density appearing at the denominator of the velocity~\eqref{eq:lim_velocity} and of the specific internal energy~\eqref{eq:lim_eps} is itself limited by any of the two different cases: \begin{itemize} \item a hierarchical strategy $\lambda_{\rho}^l$ for $l^{th}$ order term $\nabla^l \rho$: \begin{equation} \label{eq:limi_hierarch_rho} \hspace{-1cm} \mytextcolor{red}{P^{\lambda}(\mathbf{x},\rho)}= \bar{\rho}_j^{\lambda} + \sum _{l=1}^n \frac{1}{l!} \lambda^l_{\rho}(\nabla^{l} \rho)_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j), \end{equation} with \begin{equation} \label{eq:cte_highorder_lim_hi} \bar{\rho}_j^{\lambda} = \rho_j - \sum _{k=1}^{n} c^{k,\lambda}_j(\rho) \end{equation} and: \begin{equation} \label{eq:cte_correct_highorder_lim_hi} c^{k,\lambda}_j(Q) \coloneqq \frac{1}{k!} \frac{1}{|\Omega_j|}\int_{\Omega_j}\lambda^k_{Q} (\nabla^k Q)_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j) \dd \mathbf{x}. \end{equation} \item or the same function $\lambda_{\rho}$ is used for all high-order terms (in the same way as~\eqref{eq:lim_eps}~\eqref{eq:lim_velocity}): \begin{equation} \label{eq:limi_globally_rho} P^{\lambda}(\mathbf{x},\rho)= \bar{\rho}_j^{\lambda} + \lambda_{\rho} \sum _{l=1}^n \frac{1}{l!} (\nabla^{l} \rho)_j (\mathbf{x}-\mathbf{x}_j,\ldots ,\mathbf{x}-\mathbf{x}_j), \end{equation} with (see~\eqref{eq:cte_correct_highorder}): \begin{equation} \label{eq:cte_highorder_lim_all} \bar{\rho}_j^{\lambda} = \rho_j - \lambda_{\rho} \sum _{k=1}^{n} c^{k}_j(\rho) \end{equation} \end{itemize} Here, $\lambda_\rho$ in~\eqref{eq:limi_globally_rho} or each $\lambda_{\rho}^i$ in~\eqref{eq:limi_hierarch_rho} are designed to obtain at least the positivity criterion \\$\mytextcolor{red}{P^{\lambda}(\mathbf{x}_l,\rho)} > 0$ for all $\mathbf{x}_l$ quadrature points of $\partial \Omega_j$ but also~\eqref{eq:TVB} (with $Q=\rho$). %Note that we could also use the construction of positive polynoms of arbitrary order. %\end{remarque} \eject \begin{remarque} We recover exactly the second and third-order cases of~\cite{PH_NodalSolv:hal-03585115}. In~\eqref{eq:lim_eps}~\eqref{eq:lim_velocity}, we have the following expression for %$\hat{P}$ \mytextcolor{red}{$HP\big(\mathbf{x},[\rho]_0^{n-1},[\frac{1}{2}\big|\mathbf{U}\big|^2]_1^{n-1} \big)$} % polynomials (order two homogeneous) \hspace{-1cm} \begin{itemize} \item For n=1: $\mytextcolor{red}{HP\big(\mathbf{x}\big)}\equiv 0$ \item For n=2: \begin{multline*} HP\big(\mathbf{x},\rho_j, (\nabla \mathbf{U})_j\big) \\ =\frac{1}{2} \,^t(\mathbf{x}-\mathbf{x}_j) (\rho_j (^t\nabla \mathbf{U})_j (\nabla \mathbf{U})_j ) (\mathbf{x}-\mathbf{x}_j) - \frac{1}{2}\frac{1}{|\Omega_j|} \int _{\Omega_j} \,^t(\mathbf{x}-\mathbf{x}_j) (\rho_j (^t\nabla \mathbf{U})_j (\nabla \mathbf{U})_j ) (\mathbf{x}-\mathbf{x}_j) \dd \mathbf{x} \end{multline*} \end{itemize} which are indeed 2-homogeneous polynomials w.r.t high-order velocity terms. \end{remarque} \subsection{Choices of cell by cell limiter \texorpdfstring{$\lambda_{\rho}, \lambda_{\epsilon}$}{lambda rho, lambda epsilon}} In view of~\eqref{high_order_transport_direct_euler_composite}, we need to evaluate high-order fluxes wherever we consider Riemann problems. Many limiter for the density~\eqref{eq:limi_globally_rho} and for the internal energy~\eqref{eq:rec_eps_lim} (that induce velocity limiters~\eqref{link_epsU_lim}) can be investigated. Here we use a minmod/Barth--Jespersen like limiter constructed in the spirit of~\eqref{eq:limi_globally_rho}. For each cell we define a scalar cell limiter $\lambda^{\rho}_j$: \begin{itemize} \item $\mytextcolor{red}{P^{\lambda}(\mathbf{x},\rho)}$, the high-order limited reconstruction of the density: \begin{equation} \label{eq:def_limiter_density_all} \mytextcolor{red}{P^{\lambda}(\mathbf{x},\rho)} = \rho_j + \lambda^{\rho}_j (\mytextcolor{red}{P(\mathbf{x},\rho)} -\rho_j). \end{equation} Now, we explain the definition of $\lambda^{\rho}_j$: \begin{enumerate} \item\label{Sect3.1_1} Evaluation of (unlimited) $\mytextcolor{red}{P(\mathbf{x},\rho)}$ at $\mathbf{x}_l$ coordinates of quadrature points (on boundary cell $\Omega_j$ nodes and edges) where theRiemann problem is posed, see~\eqref{high_order_transport_direct_euler_std},~\eqref{high_order_transport_direct_euler_composite}: the associated values \mytextcolor{red}{$P(\mathbf{x}_{l},\rho) $} are denoted by $P_{j,l}$. We also define: \begin{equation} \label{eq:MinMaxEvalRecDofCell} \begin{cases} M^{\rho,\partial \Omega_j}_j\coloneqq \max _{\mathbf{x}_l \in \partial \Omega_j} P_{j,l} \\ m^{\rho,\partial \Omega_j}_j\coloneqq \min _{\mathbf{x}_l \in \partial \Omega_j} P_{j,l} \end{cases} \end{equation} \item\label{Sect3.1_2} Evaluation of (unlimited) \mytextcolor{red}{$P(\mathbf{x},\rho)$} at $\mathbf{x}_{j'}$ centro\"id of neighboring cell $\Omega_{j'}$, $j'$ belongs to the neighborhood $V_j$ of $j$ ($j \not \in V_j$) see Figure~\ref{Fig: V_j}, the associated values \mytextcolor{red}{$P(\mathbf{x}_{j'},\rho) $} are denoted $P_{j,j'}$, we also define: \begin{equation} \label{eq:MinMaxEvalRecCentNeighbCell} \begin{cases} M^{\rho,V_j}_j\coloneqq \max _{j' \in V_j} P_{j,j'}, \\ m^{\rho,V_j}_j\coloneqq \min _{j' \in V_j} P_{j,j'}. \end{cases} \end{equation} \item\label{Sect3.1_3} Computation of extremal values of previous reconstructed values: \begin{equation} \label{eq:MinMaxAllVal} \begin{cases} M\rho_j \coloneqq \max\left(M^{\rho,\partial \Omega_j}_j,M^{\rho,V_j}_j\right), \\ m\rho_j \coloneqq \min\left(m^{\rho,\partial \Omega_j}_j,m^{\rho,V_j}_j\right). \end{cases} \end{equation} \item\label{Sect3.1_4} Computation of extremal bounds of solution in neighborhood $V_j$: \begin{equation} \label{eq:MinMaxValNeighbCell} \begin{cases} \rho_j^{\max}\coloneqq \max(\rho_j, \max _{j' \in V_j} \rho_{j'}), \\ \rho_j^{\min}\coloneqq \min(\rho_j, \min _{j' \in V_j} \rho_{j'}). \end{cases} \end{equation} \item\label{Sect3.1_5} Evaluation of min and max gap to obtain spatial local bound preservation: %slope: \begin{equation} \label{eq:SlopeMinMax} \begin{cases} s^{\rho,\max}_j\coloneqq \begin{cases} 1, &\text{if } |M\rho_j -\rho_j| \le \delta, \qquad (\delta \ll 1)\\ \frac{\rho_j^{\max}-\rho_j}{M\rho_j -\rho_j}, & \text{else}. \end{cases} \\ s^{\rho,\min}_j\coloneqq \begin{cases} 1, &\text{if } |m\rho_j -\rho_j| \le \delta, \\ \frac{\rho_j^{\min}-\rho_j}{m\rho_j -\rho_j}, &\text{else}. \end{cases} \end{cases} \end{equation} \end{enumerate} Finally, we take: \begin{equation} \label{eq:lim-density-unique-cell} \lambda^{\rho}_j = \max(0,\min(1, \min(s^{\rho,\min}_j,s^{\rho,\max}_j))). \end{equation} \item $\mytextcolor{red}{R^{\lambda}(\mathbf{x},\epsilon)}$, the high-order limited reconstruction for internal energy (see~\eqref{eq:lim_eps}), we adopt the same approach replacing $\mytextcolor{red}{P^{\lambda}(\mathbf{x},\rho)}$ (resp. $\rho_j$) by $\mytextcolor{red}{R^{\lambda}(\mathbf{x},\epsilon)}$ (resp. $\epsilon_j$). We obtain \begin{equation} \label{eq:lim-epsilon-unique-cell} \lambda^{\epsilon}_j = \max(0,\min(1, \min(s^{\epsilon,\min}_j,s^{\epsilon,\max}_j))), \end{equation} and deduce from~\eqref{eq:lim_eps} \begin{equation*} \mytextcolor{red}{R^{\lambda}(\mathbf{x},\epsilon)} = \epsilon_j + \lambda^{\epsilon}_j \mytextcolor{red}{\hat{R}(\mathbf{x},\epsilon)} \end{equation*} \item $\mytextcolor{red}{\boldsymbol{R}^{\lambda}(\mathbf{x},\mathbf{U})}$, the high-order induced limited reconstruction for velocity ~\eqref{eq:lim_velocity}: \begin{equation*} \mytextcolor{red}{\boldsymbol{R}^{\lambda}(\mathbf{x},\mathbf{U})} = \mathbf{U}_j + \lambda^{\mathbf{U}}_j \mytextcolor{red}{\hat{\boldsymbol{R}}(\mathbf{x},\mathbf{U})} \overbrace{=}^{\eqref{link_epsU_lim}} \mathbf{U}_j + \sqrt{\lambda^{\epsilon}_j} \mytextcolor{red}{\hat{\boldsymbol{R}}(\mathbf{x},\mathbf{U})} \end{equation*} \end{itemize} \begin{remarque}\label{eq:rem_induced_limitation} \:~ \begin{itemize} \item We emphasize here that other \emph{a priori} limitation strategies are possible, and will be investigated in future works. \item In case of the weaker positivity condition on a target quantity $f = \rho$ (density) or $f = \epsilon$ (specific internal energy), we just set $s^{f,\max}_j=1$ and $f_j^{\min}=0$ (or $f_j^{\min}=\eta$ meaning that $\rho \ge \eta$ (or $\epsilon \ge \eta$) for example $\eta=10^{-15}$) in~\eqref{eq:SlopeMinMax}. Nevertheless, we have noticed that in some strong shock test cases, imposing only positivity may lead to oscillations.% the only positivity enforcement may drive to more oscillations. \item Step~\eqref{Sect3.1_2} has been added to~\eqref{eq:MinMaxAllVal} and gives better results in some test cases. \mytextcolor{red}{ \item In order to obtain both $R^{\lambda}(\mathbf{x},\rho)>\rho_{\min}$ and $R^{\lambda}(\mathbf{x},\epsilon)>\epsilon_{\min}$ at each quadrature point in~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}, we need of course that the constant term (giving by finite volume unknowns at each time step) satisfies: \begin{align} \label{eq:rho_eps_validate} \rho_j>\rho_{\min} \quad\mbox{and} \quad \epsilon_j>\epsilon_{\min} \end{align} This is essentially linked to qualitative property of the underlying numerical flux of having enough numerical viscosity and also that the underlying first order scheme (with first order reconstruction) fullfills constraints~\eqref{eq:rho_eps_validate}. } \end{itemize} \end{remarque} %\vspace{-.3cm} \section{Numerical results}\label{section:num_results} \mytextcolor{blue}{ This section is devoted to numerical test cases involving high-order (second and third-order here) reconstruction with the limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim_eps}~\eqref{eq:lim_velocity}. The high-order time integration is obtained with RK3-TVD~\cite{shuosher}. For the tests we present here, we take $\delta=10^{-15}$ in~\eqref{eq:SlopeMinMax}. We use a composite extension (node and edges) of Roe (with entropy correction) and Rusanov numerical fluxes~\cite{GoRav,Toro} (see~\cite{PH_NodalSolv:hal-03585115} and~\eqref{high_order_transport_direct_euler_composite} with $\theta=\frac{2}{3}$ or $\theta=\frac{\pi}{4}$). For sanity check, we test the overall reconstruction limiter-free $\lambda^{\rho} \equiv 1$ in~\eqref{eq:limi_hierarch_rho} (or~\eqref{eq:limi_globally_rho}), and $\lambda^{\epsilon} \equiv 1$ in~\eqref{eq:lim_eps} (and then~\eqref{eq:lim_velocity}) on a regular solution then on a series of classical non regular problems. } \mytextcolor{red}{ As highlighted by the last point of remark~\eqref{eq:rem_induced_limitation}, in order to obtain physical valid constraints~\eqref{eq:rho_eps_validate} (if not fullfill by an high-order scheme), we apply an APITALI process (see~\cite[Section~4.4]{PH_NodalSolv:hal-03585115} and~\cite{ale2d_strategy,HL_APITALI_Vecteur}). For the approximation of discontinuous solutions, the a-posteriori control permits to obtain such constraints. Especially, in the next sub-section~\ref{section:effect_apriori}, in order to make comparisons, some computations are done without \emph{a priori} limiters ($\lambda^{\rho} \equiv 1$~\eqref{eq:limi_globally_rho} and $\lambda^{\epsilon} \equiv 1$~\eqref{eq:lim_eps} and then~\eqref{eq:lim_velocity}). In this case, we observe the following qualitative behavior: \begin{enumerate} \item If we use \emph{a posteriori} strategy by just imposing positivity~\eqref{eq:rho_eps_validate}, the numerical solution at final time exhibits some oscillations (see Figure~\ref{Fig:1D_ComparAprioriLimUnlim:Sod} and Figure~\ref{Fig:1D_ComparAprioriLimUnlim:Noh}). \item If we do not use \emph{a posteriori} strategy, the code may crash in very few time steps by non physical states. \end{enumerate} } %\vspace{-.3cm} \subsection{Third-order checking on regular Taylor--Green Vortex solution} This standard test case is used to assess the expected (high-) order of a numerical method. The domain is $[-5,5] \times [-5,5]$, and the analytic regular solution is given by: \begin{equation} \rho(t,\mathbf{x})=(T_{\infty}+\delta T)^{1/(\gamma-1)},\:\: \mathbf{U}(t,\mathbf{x})=\mathbf{U}_{\infty} + \delta \mathbf{U},\:\: p(t,\mathbf{x})=\rho^{\gamma}, \end{equation} with \begin{equation} \begin{aligned} %\mbox{with}\qquad \delta \mathbf{U}(t,\mathbf{x})&=\frac{\beta}{2\pi} \exp(1-r^2) (-x_2,x_1)\\ \delta T(t,\mathbf{x})&=-\frac{(\gamma-1)\beta^2}{8 \gamma \pi^2} \exp(1-r^2)\\ \end{aligned} \end{equation} and $T_{\infty}=1$, $\mathbf{U}_{\infty}=\mathbf{0}$ (stationnary), %$T_{\infty}=1$, $\beta=5$, $\gamma=1.4$, $\mathbf{x}=(x_1,x_2)$ and $r^2=|\mathbf{x}|^2$, $t_f=1$. \mytextcolor{blue} { We note (see Figure~\ref{fig:TL_comp}) that the strategy based on gathering the high-order terms in reconstructed fields: velocity $\mathbf{U}$ and internal energy $\epsilon$ (by the way of Leibniz formula) gives indeed a third order convergence (without limitation). We also compare results with the value $\theta=2/3$ (Simpson formula) to the value $\theta=\pi/4$ (the conical degenerate case, see~\cite{PH_NodalSolv:hal-03585115}). Although only the former gives a third/fourth order quadrature formula for the fluxes, we observe for both cases approximately the same error and order convergence rate. } \begin{figure}[!htb]\centering \includegraphics[scale=.3]%[height=5.5cm,width=11cm] {cv_vo_teta23.pdf} %\\ \includegraphics[scale=.3]%[height=5.5cm,width=11cm] {cv_comp_pi4_23.pdf} \caption{\label{fig:TL_comp}Convergence of third order composite Rusanov limiter-free with $\theta=2/3$ and $\theta=\pi/4$.} \end{figure} \subsection{1D test cases with non-smooth solution} \subsubsection{Sod Shock Tube} \mytextcolor{blue} { The standard Riemann problem of Sod shock tube is solved to asses the good behavior of limitation process, conservation property (discontinuity location and level of intermediate states for weak solutions) of the schemes. It initialy consists of a perfect gas ($\gamma=1.4$) at rest with left state %$(\rho,\mathbf{U},P)_{Left} = (1,\mathbf{0},1)$ and a right state $(\rho,\mathbf{U},P)_{Right} = (.125,\mathbf{0},0.1)$. $(\rho,u,P)_{Left} = (1,0,1)$ and a right state $(\rho,u,P)_{Right} = (.125,0,0.1)$. } \newpage \mytextcolor{blue} { For the Sod shock tube (cf Figure~\ref{Fig:1D_sod_o3}), we observe a good behavior of the induced limitation for the third order reconstruction, especially for the velocity field for which the limitation is linked with those of internal energy (cf~\eqref{eq:lim_velocity}). In Figure~\ref{Fig:1D_sod_cvo123}, for two different numerical composite fluxes, we observe the benefits made by increasing the order (one, two and three) of the reconstruction along with the induced limitation. } \begin{figure}[!htb] \begin{center} \centerline{ \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_density_o3_rusanov_two.pdf} \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_density_o3_roe_two.pdf}} \vspace*{5mm} \centerline{ \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_internalmassicenergy_o3_rusanov_two.pdf} \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_internalmassicenergy_o3_roe_two.pdf}} \vspace*{5mm} \centerline{ \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_velocity_o3_rusanov_two.pdf} \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_velocity_o3_roe_two.pdf} } \vspace{-.1cm} \caption{\label{Fig:1D_sod_o3}\mytextcolor{blue}{Mesh convergence of third-order reconstruction and induced limitation on Sod shock tube (Rusanov (Left) / Roe (Right)): density (Top), internal energy (Center), velocity (Below)}.} \end{center} \vspace{-1cm} \end{figure} \newpage \vspace{-.2cm} \begin{figure}[!htb] \centerline{ \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_density_o1o2o3_rusanov_two_100.pdf} \includegraphics[scale=.31]%[height=5cm,width=7.5cm] {sod_cv_density_o1o2o3_roe_two_100.pdf}} \vspace{-.2cm} \caption{\label{Fig:1D_sod_cvo123}\hspace{-.1cm}\mytextcolor{blue}{Comparison on 100 cells (for density) of order's 1, 2 and 3 with reconstruction and induced limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} on Sod shock tube (Rusanov (Left) / Roe (Right))}.}% \vspace{.2cm} \end{figure} \subsubsection{Osher--Shu test case} %\vspace{-.3cm} This test case~\cite{shuosher2} requires the scheme to be both high-order to reproduce the smooth physical high frequencies but also to have good limitation process to resolve a discontinuity. The initial conditions contains a moving Mach 3 shock wave that interacts with periodic perturbations in density: \begin{equation} (\rho,u,P)= \begin{cases} (3.857143, 2.629369, 10.333333),& -5 \le x_1 \le -4, \\ ( 1+ 0.2 \sin(5 x_1), 0, 1), &-4 < x_1 \le 5. \end{cases} \end{equation} \mytextcolor{blue}{ For the Osher--Shu test case (Figure~\ref{Fig:1D_OsherShu}), we note that the convergence rate of the Rusanov scheme is slower than the convergence rate of Roe scheme. %do converge less faster than the Roe scheme. Indeed for a coarse mesh (100 cells), even for third order reconstruction, the composite Rusanov scheme damps all high frequencies of the solution ($x_1 \in [0,2]$). They are almost flattened (compared to composite Roe flux). Obviously, for finer mesh (800 cells), this difference between both schemes is less significant. } \begin{figure}[!htb] \centerline{ \includegraphics[scale=.31]{dens-shuosher-o3-rusanov-two.pdf} \includegraphics[scale=.31]{dens-shuosher-o3-roeflux-two.pdf}} \vspace{-.2cm} \caption{\label{Fig:1D_OsherShu}\mytextcolor{blue}{Mesh convergence of third-order schemes with induced limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim_eps}~\eqref{eq:lim_velocity} on Osher--Shu test case (Rusanov (Left) / Roe (Right))}} \end{figure} \newpage \subsubsection{Blast waves test case} \mytextcolor{blue}{ The blast waves test case proposed by Colella/Woodward~\cite{BlastW} assesses the robustness of the overall numerical machinery. The solution results from the interaction of two reflecting shocks. A robust limitation strategy is mandatory (if not then negative internal energy occurs quasi immediately). } \begin{equation} (\rho,u,P)= \begin{cases} (1,0,1000), &0 \le x_1 \le 0.1,\\ (1,0,10^{-2}), &0.1 < x_1 < 0.9,\\ (1,0,100), &0.9 \le x_1 \le 1. \end{cases} \end{equation} \mytextcolor{blue}{ For the blast waves test case (Figure~\ref{Fig:1D_BlastWaves}), we observe less smeared convergence of Roe scheme with respect to the Rusanov but the numerical results do not exhibit Gibbs phenomenum. However, the resolution of contact discontinuity at $x_1 \simeq 0.6$ could be enhanced by using an HLLC composite flux or by other ad-hoc technics. } \begin{figure}[!htb]\centering % \centerline{ \includegraphics[scale=.31]%[scale=.4] {blastw_densite_o3_rusanov_two_b_ref.pdf} \includegraphics[scale=.31]%[scale=.4] {blastw_densite_o3_roeflux_two_b_ref.pdf} %} \caption{\label{Fig:1D_BlastWaves}\mytextcolor{blue}{Mesh convergence of third-order schemes with induced limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} for the Colella/Woodward's blast waves problem (Rusanov (Above) / Roe (Below))}.} \end{figure} %\newpage \subsubsection{Effect of \emph{a priori} induced limitation of section 3}\label{section:effect_apriori} \mytextcolor{red}{ In this paragraph, we test the behavior with and without \emph{a priori} induced limitation strategy in~\eqref{eq:limi_globally_rho}~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} (see Figure~\ref{Fig:1D_ComparAprioriLimUnlim:Sod}). } \paragraph{Sod} \: \newpage \begin{figure}[!htb] \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_densite_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_densite_avec_sans_lim_200.pdf} } % \caption{\label{Fig:1D_ComparAprioriLimUnlim:dens} Density} %\end{figure} %\begin{figure}[!htb] \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_vitesse_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_vitesse_avec_sans_lim_200.pdf} } %\caption{\label{Fig:1D_ComparAprioriLimUnlim:velocity} Velocity} %\end{figure} %\begin{figure} \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_internalnrj_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_internalnrj_avec_sans_lim_200.pdf}} \vspace{-.2cm} \caption{\label{Fig:1D_ComparAprioriLimUnlim:Sod}\mytextcolor{red}{ Sod shock tube: From Top to Bottom: Density, Velocity, Internal Energy. Comparison (third order Rusanov composite scheme) with and without a-priori induced limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} on 100 cells (left) and 200 cells (right).}} \end{figure} \mytextcolor{red}{ \paragraph{Noh} The Noh test case need of using special care for high-order methods. This is mainly because of the internal energy is nearly zero. More precisely, for this shock tube problem, the data are $\gamma=\frac{5}{3}$, $t_f=0.6$ and the initial condition is uniform: \begin{equation} (\rho,u,P)= (1,-1,10^{-13}) \end{equation} On the left side, a symmetry is imposed, and an inflow on the right side (see Figure~\ref{Fig:1D_ComparAprioriLimUnlim:Noh}). } \begin{figure}[!htb] \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_noh_densite_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_noh_densite_avec_sans_lim_200.pdf} } \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_noh_vitesse_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_noh_vitesse_avec_sans_lim_200.pdf} } \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_noh_internalnrj_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_noh_internalnrj_avec_sans_lim_200.pdf} } \centerline{ \includegraphics[scale=.3]%[scale=.4] {compar_noh_pressure_avec_sans_lim_100.pdf} \includegraphics[scale=.3]%[scale=.4] {compar_noh_pressure_avec_sans_lim_200.pdf} } % \vspace{-.2cm} \caption{\label{Fig:1D_ComparAprioriLimUnlim:Noh}\mytextcolor{red}{Noh test: From Top to Bottom: Density, Velocity, Internal Energy, Pressure. Comparison (third order Rusanov composite scheme) with and without a-priori induced limitation~\eqref{eq:limi_globally_rho}\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} on 100 cells (left) and 200 cells (right).}} \vspace{-2cm} \end{figure} \newpage \mytextcolor{red}{ Results show the good capability of the induced limitation described in previous section~\ref{section:induced_limitation} do deal with discontinuities (see e.g. velocity or pressure in Noh test case Figure~\ref{Fig:1D_ComparAprioriLimUnlim:Noh}). We recall that no direct limitation of the velocity reconstruction is used, it is only deduce from the limitation of an ad-hoc internal energy reconstruction and~\eqref{link_epsU_lim}. An \emph{a posteriori} limitation is mandatory (resp. useless) in order to reach final time in case of the induce \emph{a priori} limitation strategy is not used (resp. used). } \vspace{-.4cm} \subsection{2D test cases with non-smooth solution} \subsubsection{Lax--Liu (case 12)} \mytextcolor{blue}{ In two dimensional case, Lax--Liu proposed some 2D-Riemann problems with initial data resulting in different Shock/Contact/Rarefaction configurations (see also~\cite{LiskaWend} for some schemes comparison). For the case 12 of~\cite{LaxLiu}, the data on each quadrants are (see Figure~\ref{Fig:2D_Lax-Liu} Top Left): } \begin{equation} \label{eq:cauchy_laxliu12} \begin{aligned} &(\rho,\mathbf{U},P)_{2} = (1,(0.7276,0),1),& &(\rho,\mathbf{U},P)_{1} = (.5313,(0,0),0.4), \\ & (\rho,\mathbf{U},P)_{3} = (.8,(0,0),1),& &(\rho,\mathbf{U},P)_{4} = (1,(0,0.7276),1). \end{aligned} \end{equation} \begin{figure}[!htb] \vspace*{-6mm} \centerline{\includegraphics[scale=.4]{config_4zones.pdf} \hspace{1.4cm} \includegraphics[scale=.94]{sol12.xcf.pdf}} % \vspace{5mm} \centerline{ \includegraphics[scale=.35,clip=true,trim = 0 150 0 150]%[scale=.44] {Riemann2D_12_Rusanov_two.pdf}\hspace{-.4cm} \includegraphics[scale=.35,clip=true,trim = 0 150 0 150]%[scale=.44] {Riemann2D_12_Roeflux.pdf}} \caption{\label{Fig:2D_Lax-Liu}\mytextcolor{blue}{Density for third-order schemes with induced limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} for two dimensional Riemann problem (configuration 12 of~\cite{LaxLiu} Top Left and~\eqref{eq:cauchy_laxliu12} with a reference solution on Top Right see in~\cite{LiskaWend}) (Rusanov (Left) / Roe (Right)) with 400x400 cells.}} \end{figure} We observe in Figure~\ref{Fig:2D_Lax-Liu} a good agreement with respect to manufactured reference solution. Here, we emphasize that for the velocity treatment, we do not need to use an ad-hoc procedure for vectorial limitation such as VIP~\cite{VIP_LF,HL_APITALI_Vecteur}. \subsubsection{Double Mach Reflection} \mytextcolor{blue}{ For the test case proposed in~\cite{BlastW}, we compare results (see Figure~\ref{Fig:2D_DbleMach}) of our proposed limitation strategy with two numerical composite fluxes Rusanov and Roe. } \begin{figure}[!htb] \centerline{ \hspace{.2cm} \includegraphics[scale=.35,clip=true,trim=0 150 0 150]%[scale=.44] {dble_mach_densite_rusanov2_tri.pdf}\hspace{-.4cm} \includegraphics[scale=.35,clip=true,trim=0 150 0 150]%[scale=.44] {dble_mach_densite_roeflux_tri.pdf} } \caption{\label{Fig:2D_DbleMach}\mytextcolor{blue}{Density for third-order schemes with induced limitation~\eqref{eq:limi_globally_rho}~\eqref{eq:lim-density-unique-cell}~\eqref{eq:lim_eps}~\eqref{eq:lim-epsilon-unique-cell}~\eqref{eq:lim_velocity} for the double Mach reflection problem with composite schemes (Rusanov (Left) / Roe (Right)) on 40201 unstructured triangles.}} \end{figure} \mytextcolor{blue}{ At final time $t_f=2$, we note that upper shock at $x_1=2$ is well resolved (without oscillation and at the right position). Here, we do not observe the ``jet'' near the oblic reflection axis that may occur with standard edge flux schemes (carbuncle phenomenum). } %\newpage \section{Conclusion} In this paper, we propose an \emph{a priori} limitation process in the reconstruction step for the d-dimensional compressible Euler equations of gas dynamics. Taking an arbitrary order polynomial reconstruction of conservative quantities of this system, we have obtained sufficient stability condition on density $\rho>0$ and on the internal specific energy $\epsilon>0$. The main idea is to use the multi-variate (arbitrary order) Leibniz formula for primal specific quantities: the velocity $\mathbf{U}$, the specific total energy $E$. The high-order kinetic energy part is deduced (again by Leibniz formula) from high-order terms coming from those of velocity. Finally, using an explicit form of the non linear rational reconstruction $R(\mathbf{x},\epsilon)$, it is possible to obtain a direct limitation (for $\epsilon$) that induces a limitation for the velocity reconstruction $\boldsymbol R(\mathbf{x},\mathbf{U})$. \section*{Thanks} %\thanks{ The author want to thanks the referee for valuable remarks and comments, helping the better understanding of the document. And also thank to S.~Del~Pino and X.~Blanc for their carefull reading of the paper. %} \section*{Declaration of Interest} The author do not work for, consult, own shares in or receive funding from any company or organization that would benefit from this article, and have declared no affiliation other than their research organisations. \printbibliography \end{document}