Duval, L L Bouard

Keywords: maillages volumiques, multi-échelle, ondelettes, compression, simulation, ingénierie des réservoirs volume mesh, multiscale, wavelets, compression, simulation, reservoir engineering Résumé

Titre Ranement en résolution et précision pour la compression et la simulation de maillages volumiques en géosciences En bref, La recherche est entrée dans une ère marquée par l'utilisation intensive des données. La quantité croissante d'informations scientiques pose un dé dans plusieurs domaines d'application. C'est le cas notamment en simulation, où l'utilisation d'énormes volumes de données créent des goulets d'étranglement lors de calculs haute performance. Plusieurs méthodes de réduction des données sont actuellement à l'essai pour tenter de résoudre ce problème. Dans ce travail, nous nous concentrons sur la modélisation géologique et le worow de simulation en ingénierie réservoir. En géosciences, les modèles sont composés d'informations hétérogènes, notamment la géométrie du maillage et les propriétés pétrophysiques. Ces derniers peuvent contenir jusqu'à plusieurs millions de cellules. Des techniques d'upscaling/upgridding sont couramment utilisées pour réduire le temps de simulation. Pourtant, elles sont souvent ad hoc et ne répondent pas entièrement à tous les besoins de manipulation des données : visualisation, stockage, et génération d'une donnée à résolution et précision adaptées pour la simulation.

Nous proposons une méthodologie complète permettant un ranement en résolution et en précision, basée sur HexaShrink, un outil de décomposition multi-échelle basé sur les ondelettes. Nous évaluons ses capacités de compression sans perte puis avec perte en le combinant avec des codeurs entropiques génériques et évolués (type "zerotree") , ainsi que sa pertinence visuelle en appliquant la méthode sur une collection de maillages représentatifs.

Nous testons aussi de manière approfondie l'impact du ranement de la résolution et de la précision sur la simulation. Tous nos tests de simulation ont été réalisés sur Lundi, un modèle représentatif de diérents environnements géologiques, généré spéciquement pour ce travail. Les résultats se comparent positivement aux codeurs de référence SZ et ZFP, en dénissant des métriques d'évaluation objectives corrélées aux évaluations subjectives des résultats de simulation réservoir.

tenant au cadre universitaire, professionnel et personnel.

Je tiens à exprimer mes plus sincères remerciements à mes encadrants. Marc ANTONINI, mon directeur de thèse au laboratoire I3S à Sophia Antipolis ainsi que Frédéric PAYAN dont le soutien dans les moments les plus diciles m'a permis de poursuivre mes travaux. Je remercie tout particulièrement Laurent DUVAL, mon promoteur de thèse IFPEN, dont le nombre d'heures consacré à ce projet et à mon bien être général a très fortement été revu à la hausse. Vous avez, à n'en pas douter, démontré des compétences dans un encadrement à distance désormais plus qu'utiles. Merci de m'avoir fait découvrir le domaine de la compression et les disciplines associées, j'espère dans ma future carrière pouvoir partager l'essence et les bienfaits de ces traitements.

Un grand merci à Christophe PREUX pour son expertise dans le domaine de l'ingénierie réservoir et à Sébastien MAGAND pour son étude du big data et de son usage en simulation venue consolider la partie contexte de ce projet.

Thanks to my reviewers and examiners for reading and reporting on the manuscript, and for the enriching exchanges during the defense. C'est avec émotion et beaucoup de regrets, en l'absence de pot de départ en raison de la crise sanitaire, que je quitte toute l'équipe IFPEN de la Direction Sciences et Technologies du Numérique de Solaize. Je tenais à vous remercier pour votre accueil chaleureux et l'atmosphère quotidienne bienveillante, notamment Laurent PIGEON pour sa patience dans la résolution de mes problèmes informatiques. Merci au bureau des doctorants pour leur étude sur la théorisation des humeurs durant ces trois années de thèse partagées et vécues aux côtés de Johan l'ancien, David, le tic de tic et tac, et les jeunes Bassel et Adam dans l'ordre de succession. J'espère tous vous revoir très prochainement.

Je tiens aussi à remercier l'équipe I3S pour avoir rendu mes déplacements à l'Université Côte d'Azur agréables, notamment : Cyprien, Melpo, Arnaud, Danny, Eva. . . i ii Pour nir, je remercie ma famille qui a su m'entourer et m'encourager : mes parents, sans aucun doute mes plus grands supporters, mon frère, ma s÷ur, leurs conjoints respectifs pour l'exemple qu'ils m'ont toujours donné, ma grand-mère, pour sa force de caractère dont j'espère avoir hérité, Philippe et Marie-Do pour m'avoir plus que choyée lors de mes séjours à Antibes, Sylvain pour m'avoir accompagnée sur cette n de thèse, ma belle-famille pour m'avoir accueillie comme leur lle/s÷ur notamment lors du premier connement, mes ami(e)s de géologie et de médecine, du collège/lycée et plus particulièrement ma Marie. Enn merci à toi mon Gautier, une étape de plus pendant laquelle tu as su me supporter, écoutant le récit quotidien des réussites et des échecs, sachant trouver les mots justes pour me rassurer et rendre le moment plus doux. On a pas ni de grandir ensemble, j'ai hâte de la suite.

Résumé

Les sciences nous aident à appréhender le monde qui nous entoure : qu'il s'agisse du monde vivant, d'objets naturels à petites ou grandes échelles, ou de phénomènes physiques proches ou lointains.

Cette connaissance du monde est nourrie par des avancées technologiques et intensivement documentée par une donnée générée en grande quantité. Face à ces volumes de données croissants les capacités d'exploitation deviennent insusantes.

Nous nous intéressons dans cette étude à la gestion des données de simulation pour les géosciences grâce à la compression, deux notions qui seront dénies et contextualisées par la suite.

Dans cette introduction nous diérencierons trois types de données scientiques courantes : la donnée expérimentale, de simulation, d'apprentissage. La donnée expérimentale est une donnée produite par des outils de mesure dont la précision n'a de cesse d'augmenter (microscope, satellite, méthodes de tomographie). Ils permettent aujourd'hui d'observer et mesurer des objets et des évènements autrefois méconnus ou inaccessibles.

Ces dernières années les progrès du numérique ont démocratisé des outils qui viennent compléter la vision naturaliste fournie par les données expérimentales. Notamment, les outils de simulation approximent des phénomènes physiques dont la structure de base et le modèle mathématique sont connus. Les résultats de la simulation peuvent être visualisés et analysés pour améliorer la compréhension d'un objet : sa formation et son devenir (courants océaniques, phénomènes météorologiques, cosmologiques). Plus récemment les méthodes d'apprentissage ont été popularisées. Ce processus se base sur l'exploitation de bases de données existantes, et ce pour dénir des liens profonds et comprendre les rouages de phénomènes complexes.

Qu'importe leur type, le domaine, d'énormes quantités de données ont été générées par des scientiques en soif de connaissance : planiant des campagnes d'acquisition colossales, et construisant des centres de calculs aux performances exponentielles. Cette frénésie a été soutenue par la découverte de matériaux innovants (intégrant les terres rares dans leurs compositions), qui ont permis la miniaturisation des outils informatiques, augmenté les moyens, et facilité cette quête du savoir. Les iii iv découvertes et connaissances profondes qui en découlent constituent actuellement le socle et l'avenir de nos sociétés technologiques.

Pourtant certains chercheurs tirent la sonnette d'alarme : les données scientiques sont générées en volumes exponentiels, des volumes qui ne semblent pas se restreindre aux vues des futurs projets.

Or la manipulation de ces grands volumes est problématique à diérentes étapes du processus. En eet elle peut altérer des puissances de calcul (simulation, apprentissage), ralentir la transmission de la donnée de son lieux de production/calcul à l'outil d'analyse du scientique (ordinateur local), sans parler des besoins de stockages croissants sous-jacents. Finalement, ces volumes de données sont visualisés/analysés pour qu'en émanent des résultats parfois diciles à cibler, qui ont tendance à se noyer dans ces volumes.

Ainsi diérentes étapes du processus sont dites pénalisantes et appelées goulets d'étranglement par la communauté scientique. De plus la gestion des volumes de données entraine un coût signicatif de moyens, qu'ils soient matériels, énergétiques et économiques. Ce problème tend à limiter le nombre et l'ampleur des projets en recherche. Pire encore, il contraint actuellement des scientiques à détruire une partie de la donnée sans même l'avoir analysée.

Pour prévenir et anticiper ce déluge imminent, diérentes solutions constituent des axes de recherche majeurs. Ils visent à réduire la quantité de donnée en suivant cinq concepts distincts dont on souligne ici l'imbrication. Premièrement le compressed sensing manipule des bases de données naturelles qu'il échantillonne (spatialement/temporellement ou de manière aléatoire), ceci an de la résumer par un sous-ensemble. Ce sous-échantillonnage contrairement aux autres méthodes contient in ne une donnée intelligible, semblable à la donnée initiale. Les prochaines notions quant à elles se basent sur la transformation algorithmique de la donnée brute. De cette phase émerge des coecients représentatifs de la donnée initiale. Il s'agit d'une phase dite de concentration. S'en suit une condensation, cette étape ne conserve que la partie émergente de la donnée transformée pour en extraire moyennes et écarts-types etc.. La compression quand à elle intègre la donnée transformée, et génère des chiers binaires, choisissant ou non de tronquer la précision binaire (quantication).

Une dernière étape, dite de compaction (encodage), assimile les phases de concentration et de compression, structurant et réorganisant l'information binaire pour réduire le volume de données.

De ces notions imbriquées découle l'intuition qu'une donnée peut être hiérarchisée, et qu'une partie restreinte mais signiante de la donnée pourrait sure aux besoins scientiques. Le but ultime serait alors de pouvoir estimer pour une donnée la quantité minimale d'information utile : à l'image de l'entropie de Shannon -notion fondamentale de la théorie de l'information.

A l'échelle de ce projet de thèse, il serait présomptueux et démesuré de vouloir satisfaire l'ensemble de la communauté scientique. Cependant nous sommes convaincus que ce travail peut servir différents domaines et inspirer de nouvelles recherches. Ainsi ce projet de thèse s'est concentré sur une donnée de simulation utilisée en géosciences. Nous avons construit une méthodologie adaptée, motivée par des besoins métiers et une connaissance de la donnée, basée sur des notions et outils de compression existants.

Nous présentons dans un premier temps la donnée sur laquelle l'ensemble de notre étude repose. Il s'agit d'un maillage volumique : une structure numérique permettant de modéliser des formations Résumé vii choisissant un encodeur évolué. Les encodeurs génériques testés jusqu'à présent n'ont tenu compte ni du volume de la donnée, ni de la décomposition multi-échelle opérée par HexaShrink.

A l'instar de JPEG2k nous avons complété la décomposition en ondelette d'un encodeur progressif type zerotree. Cet outil exploite l'organisation multi-échelle héritée de la transformation. Se dessine en eet des structures, qui ont pour origine des coecients dans l'approximation (ou sous bandes de haut niveau), là où se concentrent les fortes valeurs et s'étendent dans plusieurs sous bandes de bas niveaux vers des coecients de plus faibles valeurs. Semblables à des arbres, ces structures peuvent être remplies de zeros à la lumière de la profondeur binaire d'où le terme zerotree. Des plans de bits les plus signicatifs aux plans de bits les moins signicatifs, la donnée est progressivement encodée, et les arbres de zéros remplacés par un seul symbole. De cette manière le budget binaire est économisé et la compaction de la donnée optimisée.

Le ux binaire généré par l'algorithme zerotree est donc ordonné selon sa signiance. De cette manière la transmission même partielle d'une partie de la donnée encodée permet par décodage et transformation inverse de reconstruire une version approchée de la donnée initiale, à précision numérique inférieure. Le reste de l'information ajoute du détail lors de la reconstruction, son utilisation complète garantit une reconstruction de la donnée à l'identique.

Aux résultats obtenus avec des encodeurs génériques (LZMA, bzip2, gzip) sur des données pétrophysiques [valeurs ottantes], nous comparons ceux générés grâce au zerotree. Bien qu'intéressants les taux de compression en sans perte sont limités comme le relaye la communauté scientique.

Par la suite, nous décidons de changer la précision de nos propriétés pétrophysiques. Nous poursuivons dès lors l'hypothèse initiale qu'une partie de l'information peut être susante pour des besoins scientiques en jouant sur les notions de précision numérique ou résolution spatiale.

En modiant la précision d'une donnée on modie sa qualité. Ce changement peut être évalué de manière objective ou subjective. L'évaluation objective se base sur un calcul mathématique, mesurant la diérence entre la donnée originale et la donnée dégradée. L'évaluation subjective quant à elle se base sur l'opinion d'un panel d'observateurs non-qualiés. Selon leur appréciation (visuelle, auditive..), une note subjective est attribuée à la dégradation.

Contrairement à la donnée multimédia, la donnée scientique peut servir de support à la simulation (notre cas). Par conséquent son altération n'est pas seulement visuelle, mais son impact doit aussi être évalué dans un workow de simulation et estimé par des professionnels. En complément des métriques objectives classiques, nous développons de nouvelles méthodes adaptées an d'évaluer avec justesse la qualité d'une donnée scientique. La confrontation des résultats obtenus avec diérentes métriques a pour but de valider ou rejeter certaines méthodes d'évaluation, plébiscitées par des travaux sur la compression de données scientiques.

A propos de la littérature, diérents outils de réduction et compression ont déjà été intégrés avec succès à des workows de simulation dans diérents domaines (climatologie, cosmologie). La réduction et la compression avec perte permettent de considérablement réduire la quantité d'information binaire, et ainsi soulager diérents points clés de la simulation HPC. L'application de ces méthodes reste encore marginale, car son impact est dicile à évaluer. Cependant, le sujet gagne en popularité au vue de récentes publications, du partage de benchmark et de la diusion de codes open source. De plus de nombreuses conférences eurissent à ce sujet. viii Pour poursuivre notre étude nous avons créé un workow de simulation (écoulement en milieu poreux) pour évaluer dans un premier temps la méthode multi-échelle HexaShrink, et étudier dans un second temps l'impact d'une précision ranable.

Nous construisons un cas d'étude comportant un réservoir faillé, baptisé lundi, sur lequel nous paramétrons une simulation d'écoulement, classiquement réalisée en ingénierie réservoir pour simuler l'exploitation d'hydrocarbures. Étant amené à répéter un nombre de fois conséquent la même simulation, les dimensions de lundi (128 × 128 × 32 cellules) sont restreintes, pour ne pas être contraint par les temps de simulation excessifs. Aussi les trois failles que comporte lundi permettront de tester HexaShrink, et d'évaluer ses performances d'upgridding. Concernant les propriétés pétrophysiques deux environnements de dépôt ont conditionné la modélisation de quatre jeux distincts, chacun composé de deux propriétés volumiques : porosité -perméabilité. Les deux environnements s'inspirent de la donnée SPE10 [START_REF] Christie | Tenth SPE comparative solution project: A comparison of upscaling techniques[END_REF], un modèle réservoir utilisé pour évaluer l'upscaling.

L'environnement uviatile comporte des objets saillants, appelés chenaux. Or l'upscaling de cette propriété est périlleux, car la modication des chenaux peut considérablement modier l'écoulement des uides dans le réservoir. Les trois autres environnements sont moins contraignants, la distribution spatiale des propriétés y est plus diuse. lundi est modélisé à l'image d'un quarter ve spot model : un cas d'école en ingénierie réservoir.

La forme générale du maillage est marquée par une pente, dont on peut évaluer le pendage le long d'une de ses diagonales en surface. Pour simuler l'exploitation d'un champ pétrolier, on place au coin le plus bas un puit injecteur opposé à un puit producteur, situé sur le coin de plus haute altitude.

Initialement le reservoir contient deux phases verticalement séquencées : une phase huileuse surplombant une phase aqueuse. L'huile située dans la partie haute du reservoir est extraite par l'injection d'eau dans la partie basse par le puit injecteur. Sous pression, la phase aqueuse pousse l'huile jusqu'au puit producteur où elle est extraite.

Les paramètres de simulation (pression, débit de production) restent xes pour chaque expérience an d'assurer une continuité, et pouvoir comparer les résultats des simulations.

Les résultats sont des paramètres physiques mesurés aux puits qui renseignent l'ingénieur réservoir de l'avancement de la production. Dans cette étude, nous nous concentrons sur la courbe water cut". Elle évalue l'évolution de la saturation en eau des volumes extraits durant l'exploitation du champs. Un remaniement des données d'entrée peut modier le prol de la courbe. Son analyse sert à valider des méthodes d'upscaling. En eet, le changement opéré sur le maillage doit peu ou ne pas impacter les résultats de simulation obtenus avec un maillage haute résolution.

Les diérentes versions de lundi, les plus basses résolutions générées avec HexaShrink ainsi que la donnée à précision ranable sont évaluées selon leur courbe water cut". Dans la littérature, s'agissant des études d'upscaling, les courbes sont tracées pour permettre à un oeil expert de se faire une opinion quant à la méthode utilisée. Ici nous proposons d'automatiser cette phase de validation en dénissant des enveloppes d'acceptabilité autour de la courbe water cut" référence. De l'enveloppe la plus proche identique" à une enveloppe plus éloignée, nous dénissons cinq critères subjectifs pour évaluer les résultats de simulation.

Quatre niveaux de décomposition HexaShrink sont appliqués au maillage lundi. Lors du pro-Résumé ix cessus de reconstruction, quatre résolutions intermédiaires sont générées et testées par simulation.

Comme détaillé avec la méthodologie HexaShrink, le maillage à la i ème résolution, noté res.-i est issu de la donnée d'approximation du i ème niveau de décomposition. Ce maillage contient (2 3 ) i fois moins de cellules que la donnée initiale. Nous utilisons les propriétés pétrophysiques de l'environnement dius, et observons l'impact de l'upgridding/upscaling d'HexaShrink sur les résultats de simulation. L'impact pour res.-1 est modéré, l'évaluation subjective donne des résultats acceptables, mais ne le sont plus à partir de la res.-3. Le temps de calcul de simulation diminue considérablement avec le nombre de cellules, passant de 2 heures pour la haute résolution à un calcul quasi instantané pour la res.-3. Ainsi nous démontrons que la décomposition multi-échelle HexaShrink obtient des résultats acceptables et diminue considérablement le temps de calcul de la simulation. D'avantages de tests pourraient conforter son utilisation pour manipuler des grilles peu complexes en ingénierie réservoir.

Dans un second temps nous étudions l'impact de la compression avec perte dans notre workow de simulation. Nous traitons dans cette expérience la haute résolution (128 × 128 × 32) de la perméabilité de lundi, à laquelle on applique un compandor, avant de procéder à la transformation en ondelettes, et générer la donnée à diérents niveaux de précision (zerotree, quantication). L'utilisation d'un compandor est rééchie et se justie par des pratiques courantes de modélisation basées sur les lois de la physique.

Pour visualiser nos résultats on diminue la précision numérique de nos données dont on rapporte l'évaluation objective en fonction de la quantité binaire de donnée compressée utilisée. La courbe est tracée pour diérentes métriques objectives et complété par l'évaluation subjective des résultats de simulation. Ainsi on vérie que l'évaluation objective de la donnée concorde avec la qualité des résultats de simulation obtenus avec. Nous concluons sur le fait que certaines métriques, notamment absolues, ne sont pas corrélées aux résultats de la simulation.

Nous démontrons que l'utilisation d'un compandor améliore la qualité objective/subjective et les performances d'un outil de compression. De plus notre approche est comparée à d'autres outils (SZ [START_REF] Cappello | Use cases of lossy compression for oating-point data in scientic data sets[END_REF], ZFP [START_REF] Lindstrom | Fixed-rate compressed oating-point arrays[END_REF]), connus et reconnus, optimisés pour la compression de données scientiques en HPC. Leurs résultats sont comparables à ceux obtenus avec notre approche.

Aussi, nous avons étudié l'anisotropie des données en géologie, ce pour améliorer notre méthode de compression. Les roches décrites par les propriétés pétrophysiques sont en eet façonnées par des évènements naturels anisotropes (sédimentation soumise à des courants orientés, pression lors de l'enfouissement du matériel, autres évènements tectoniques). Le considérant nous avons testé une méthode d'ondelettes en paquet, orientée selon l'anisotropie de la propriété. Cette décomposition particulière nécessite l'adaptation de l'algorithme zerotree, étudiée dans ses travaux par Christophe et al. (2008). Cette prise en compte de l'anisotropie améliore davantage la compression de cette propriété.

Dans un dernier temps à l'image de nouveaux outils (TTHRESH (Ballester-Ripoll et al., 2019))

qui étudient les seuillages optimaux des méthodes de compression (SZ, ZFP), nous confrontons l'ensemble de nos résultats objectifs/subjectifs, avec l'évolution de paramètres de compression mesurés Cette approche appliquée aux géosciences peut être adoptée dans d'autres domaines scientiques, prenant toujours en compte les besoins métier, et la connaissance de la donnée manipulée.

Contents

Remerciements i

Résumé iii

Publications and communications xi

Contents xii

Acronyms xvii 1 General introduction 1 Technological advances have raised our knowledge and structured the scientic process into a workow. Its dierent steps are fed by data, notably numerical data used for simulation, the topic of this study. However, whatever the stage of the workow, the data quantity increases exponentially. This is demonstrated by the endless deployment of data centers and the use of ever-increasing computing resources. Today's scientic projects process a thousand times larger data quantities than fteen years ago. The energetic, economic and ecological impact of their use encourages scientists to improve data processing by reducing data quantity. Several methods exist, in particular compression, widely developed for processing multimedia data. Some of these methods improve the gain by adjusting the precision of the data. This change is evaluated by quality metrics.

2

Chapter 1. General introduction

Researchers seek to understand the surrounding world using technological tools. Innovations from measuring instruments and intensive development of computational resources allow us to better apprehend our environment. Knowledge in diverse scientic elds leaped forward: from materials sciences to cosmology, for climate analysis or geosciences. Objects and phenomena at variable dimensions, dicult to study in the past, are now accessible through digital technologies. Prior to their development, sciences were documented by empirical descriptions and theoretical models. Numerical tools has considerably broadened the eld of possibilities by rst enriching measurement tools with digital components creating ever more sensitive sensors. Then, the rise of computational resources made it possible to simulate complex natural phenomena such as multi-phase ows, used to study large objects as hurricanes, the ow of liquid phases in the subsoil or on smaller scale air ow around a sail.

This constant improvement of tools and methods considerably increased the volume of experimental data. The growing number of exascale 1 projects in various scientic elds demonstrates our insatiable thirst of knowledge. While progress and discoveries seem innite, they are limited by concrete aspects as handling of huge data quantities. From public side, the actual cost of what is called data" is largely ignored and even growing awareness among scientists, high precision in quest for realism is the baseline from modeling to simulation, and leads to exponential rise of data volume and computing resources. The increase in data may give hope for an increase in the volume of knowledge, however, it could drown out the information.

Regarding simulation, the topic of our study, beyond diculties in processing huge data quantities, it is foreseen that calculation supports (as supercomputers) are reaching some limits. While computational resources are improving rapidly for years, storage and bandwidth capacities are actually more limited. It consequently limits several steps of the simulation workow and slows down the overall computation performances. Later in this section, we contextualize simulation within numerical environment and communicate the real costs to highlight the emergency of the situation.

To provide partial answers to the problem of data quantity in science, we focus on compression, popularized in the 80's in the multimedia context. We adapted and evaluated methods to propose a data representation with tunable precision for needs in reservoir simulation.

Data feeds scientic workow

To address scientic issues, huge data quantities are daily generated and analyzed. Among them, there are dierent types with distinct structure and binary format. Whatever the type, we nd them at distinct steps of scientic workows. Their increasing consumption is illustrated by the worldwide development of ever more data centers. 1 Adjective refers to the computing resources, among 10 18 oating operations per second delivered by supercomputers.

Data types

For scientic data, the widely used standard/norm is IEEE-754, used for calculations at oatingpoint precision. Behind this standard is the number of bits required to write each value. For simple oating-point precision, a value is coded on 32 bits, as illustrated with π example on Figure 1.1.

Among the 32 bits, one bit is dedicated to the sign while eight others are used for exponent and 23 remaining correspond to the mantissa. They allow to write a number with seven digits after the decimal point. Other extended standards exist, they multiply the binary size by two or four.

They are respectively called double or quadruple precision. This use increases the dynamic range and bring respectively a decimal precision of 16 and 34 digits. It is notably used in simulation for ever more accurate calculation needs, but is currently questioned by scientists. Because for certain applications a xed precision may not be required and moreover would increase data volumes. For several years concurrent formats of IEEE-754, as POSIX, adjust the length of mantissa to change the precision of the data. However habits die hard, because change in standard would involve to redesign the data chain from visualisation to storage such as the data workow. Throughout classical workow, three types of data can naturally be considered: experimental data, simulated data and learning data", distinguishable because of their production. Experimental data is the result of a measure or observation made on real objects or phenomena. Considering geological eld, seismic and well data are precious and unique because respectively recovered during huge and expensive acquisition campaigns and drilling operations. By contrast simulated and learning data can be generated on request by code execution but are nonetheless expensive as we will see later. Simulated data yields realistic representation of natural phenomena as ow motion description. While learning data is composed by data about real events and used to train a machine for decision-making tasks, or generated by training during the process learning example. Three data types are highly linked and feed entire workow. Experimental data are used in a rst step for modeling and then inform simulation, while learning processes can enrich and correct the original model.

From the three types of data we can use in a workow, we focus on simulation and dierentiate the inputs, the pending data (used for calculation) and the outputs. For simplication, we interest rst in simulation inputs. The objective is evaluate the impact of compression on simulation by Chapter 1. General introduction using adapted metrics, or in other words to study the simulation integrity at renable precision.

Returning to a more general point of view, given the quantities of data, it is impossible to ignore the digital infrastructure and overall data management procedures.

Digital infrastructures for data handling

Digital data is handled in the physical world, stored on the disk, loaded into memory, transmitted via optical ber and analyzed by using CPU-GPU resources. The hardware layer is of great importance in studying the needs for data processing. The most representative infrastructures doing the job are the data centers, whose number has increased considerably in recent years. According to datacentermap.com (updated in 2020), there are around 4500 data centers worldwide in 122 countries, including more than 1800 in the United States. As for France, there are 156 data centers, 35 % of which are in the Ile-de-France (Paris) region.

A data center hosts a set of numerical infrastructures for storage, transfer and computing of data.

Combined to supercomputers one can perform High Performance Computing (HPC) to manage massive data quantities and simulate objects or complex phenomena. Handling and processing are now only evaluated in terms of oating-point operations per second (ops) only. For instance, the most powerful supercomputer named Fugaku (Japan) is able to realize more than 400 billiard of oating-point operations per second (400 petaops). To a lesser extent, the supercomputer used for simulation during our study has 440 teraops capacity. It is located at IFP Energies nouvelles -Solaize and photographed on Figure 1.2. The website top500.org references performances of worldwide largest supercomputers used for exascale projects 2 .

The concern comes from the increasing number of data centers growing in size and resources with each new generation. Consumed power and its impact is becoming a rst order matter.

Infobesity: source variety & growing production

The objective is to extract scientic insight from computed results in the large HPC facilities. According to the famous statement of mathematician Richard Hamming: The purpose of computing is insight, not numbers. All these resources and their increase leads to fears of infobesity. By denition, it consists in overload data at the risk losing the data essence and information.

Moreover having previously spoken about ops, we must be aware that the entire power is not exclusively dedicated to the pure mathematical power, but distributed between various handling steps. First of all the memory context is transferred from the HPC through the network, to be saved in a reliable storage [START_REF] Naksinehaboon | Reliabilityaware approach: An incremental checkpoint/restart model in HPC environments[END_REF], waiting for its analysis in a next step. The required time is mainly underestimated, and the storing cost as well. Obvious solution for speeding up the computing process is to increase supercomputer resources and op capacities. Nevertheless such approach saturates global digital infrastructure and reveals system limitations and bottlenecks.

2 https://www.exascaleproject.org/ registers and solutions american exascale projects. Additionally in simulation because the data can be regenerated on demand, the cost it represents and the notion of value disappears. The codes are therefore executed daily and data is produced massively, with great precision in a permanent quest for realism. Intermediate results, often voluminous, raise storage issues. Therefore some HPC simulation run are merely destroyed, because their retrieval is not feasible. This represents both an energetic waste, and a loss in simulation history and analysis.

It is crucial nowadays to consider the full data processing and report the actual price caused by the infobesity. Indeed, it seems inconceivable within the current trend to ignore the ecological cost of the so-called digitalization and its handling dependencies their energetic consumption (for cooling computers) and their CO 2 emissions for the user-cost.

Our concerns are the more or less hidden cost of all these operations. Let us also briey mention, from a material point of view, the tensions around the rare earths necessary to produce ever powerful computers and the around infrastructures. It contributes to the deterioration of the already complex geopolitical situations in producer countries. Cost of digital is therefore a growing concern for energy, economic as well as ecological and political reasons.

Scientic projects with oversized data

Here are two examples of experimental projects collecting huge data volumes analyzed by Big Data, from the micro to the macroscopic scale. We talk about genomic investigations and go next to the This is not only question of data size, but these large volumes result in multiple sub-problems throughout the workow.

Bottleneck & issues in data workow

In pushing the process and level of detail for simulation, scientists have pointed out some limits for data handling. While these issues were known and investigated for data visualization, storage and transmission, other steps have recently been identied as bottlenecks to HPC performance.

Energetic, economic and ecological cost

Let's start with some general gures successively considering BigData, HPC and data centers. Current forecast predicts the eld of BigData analytic will represent a world market of 90 billion dollars 4 in 2025, and experiences a very strong growth of 20 % in the next ve future years 5 . In regard HPC market, it is estimated to 60 billions dollars in 2025 6 with a growth of 7 %7 .

8

According to recent gures on global energy consumption, digital has been estimated at 1.9 % in 2013 and 3.3 % in 2020 with a forecast growth of 9 % per year [START_REF] Ahvar | Estimating energy consumption of cloud, fog and edge computing infrastructures[END_REF]. Only considering data centers, facilities have huge energy requirement as cooling post representing among half of their needs. It represents 1 % of electrical energy consumption in 2005, 1.8 % in 2012 and potentially will reach 5 % in 2030. This alarming observation is fueled by new technologies as machine learning (García-Martín et al., 2019). In the ood anticipation research is increasingly focusing on the benets of compression for saving space in data centers, and therefore energy [START_REF] Raïs | Leveraging energyecient non-lossy compression for data-intensive applications[END_REF].

Germany is actively working on this ecological and economical problematic [START_REF] Geveler | Future data centers for energy-ecient large scale numerical simulations[END_REF] in view of energetic transition of the country. Lossless compression Original data is perfectly reconstructed after compression/decompression.

It includes all the archiving techniques, generic such as the zip, and others, dedicated to data types.

It is used to store the data and recover it without any loss (documents, codes etc.). Compression gains are generally limited to a factor of 23.

Lossy compression This approach is widespread in multimedia (audio, pictures, video) because 1.3. Numerical diet for (scientific) data? 9 useful compression rates are mostly higher (820 fold in audio or image) and thus better adapted to the usage (cf. the above example with HD video). Lossy compression implies a variable data loss, which means the data can not be perfectly reconstructed, and the decompressed data will not be exactly similar to the initial data. Lossy compression involves a control of the data loss with respect to the compression ratio. Using the imperfectness of human perception (vision or audition in the multimedia context), data modications are allowed, at locations or in frequencies where they are not perceptible, or harmless for interpretation. The MP3 format, for instance, degrades frequencies barely audible to the human ear, or sounds masked by others. Concerning images, edges and textures are especially preserved, because the human eye is highly sensitive to their degradation.

Hence, compression gains are much higher, while respecting (hopefully sucient) visual or auditive quality. 9 Mathematical estimation of the minimal binary quantity required to transmit a sample of an information.

Data compression components

Chapter 1. General introduction

The second step is the quantization, that converts the transformed data into a nite number of symbols. By simplication of small coecients, sparsity of the data is further increased to reduce binary size of data and facilitate its handling. This use is for lossy compression only (MP3 example).

It signicantly improves compression performance while modifying precision of decompressed data.

Then, a coding step exploits the remaining redundancy in the quantized data, to further compact it into a binary le. It can be performed with quantities of entropy coders, the most popular being developed by [START_REF] Human | A method for the construction of minimum-redundancy codes[END_REF]. We could also cite dictionary approach, whose the best representatives belong to Lempel-Ziv family (LZMA for example, used in Subsection 3.4.1). Roughly, all these methods look for replacement of repeated symbols by shorter one to compact the information, and nally to save bit-budget.

To better compress data acquired at various sampling rates, and to allow ecient decompressions at dierent resolutions (think about screen sizes for seamless video viewing), recent standards have investigated data transformation using the concept of multiresolution. It permits to represent the data at dierent embedded (lower) resolutions (corresponding to low frequencies or approximations), complemented by several detail components (or high frequencies). Details and approximations are combined to reconstruct higher resolutions. Wavelet transforms are instances of such techniques, as applied in our study (cf. Section 3.2).

As said, multimedia already possess many compression standards, well-known to the general public: MP3 for music (developed in the MPEG-1/MPEG-2), JPEG for images, based on a discrete cosine transform [START_REF] Wallace | The JPEG still picture compression standard[END_REF], and all the MPEG-derived standards for video: H264, HEVC for the more recent ones [START_REF] Ohm | High eciency video coding: The next frontier in video compression[END_REF]. While some of these standards are old, possibly less ecient than more recent techniques, they remain de facto standard, because they are widespread. For instance, JPEG2000 [START_REF] Taubman | JPEG2000: Image Compression Fundamentals, Standards and Practice[END_REF] previously mentioned outperforms JPEG in quality, through its wavelet-based implementation, but did not nd its expected success in mainstream multimedia. In other words, maximum quality is not the only driver for acceptance. Usability is also an important criterion.

Classical performances evaluation

Beyond usability, scientist should be picky about other measurable criteria to compare methods and evaluate their performance. In this section, main metrics are listed according to the type of compression.

Lossless compression is assessed in term of execution speed and compression ratio (CR). Last parameter is related to the bits rate, which consists in the bit number used to code a sample of compressed data. Depending on the user's needs, one of the parameters may be more important.

If the need is storage, emphasis is given to compression performance because as previously saw in Subsection 1.2.3, space is precious and expensive. But when considering transfer and visualization, execution times become the main concern. For these applications compression and decompression should have a minimal impact on overall time. Some methods, even though based on same notions 1.3. Numerical diet for (scientific) data? 11 such as Lempel Ziv (dynamic dictionary) tools, obtain very dierent performances, due to the optimization of particular components. For example, LZAM is very competitive in term of compression performance compared to LZ4. But LZ4 is the faster of the two [START_REF] Szorc | Better compression with zstandard[END_REF].

Lossy compression removes information during quantization step to improve compression performance. It obviously degrades data obtained after decompression. This introduces a dilemma: the better the compression the heavier the alteration due to the compression (as illustrated on Lena, Figure 1.4, using jpeg compression standard). While data quantity must be minimal, a certain quality remains appreciable. By consequence, it appears that quality assessment is crucial with lossy compression, and choice of metrics could be more or less adapted to evaluate alteration and to guide suitable compression. There are two ways to assess the quality of decompressed data. First one is objective and based on mathematical principles, while the second is subjective and focuses on opinion of a panel of observers. Their role is to provide marks according to the degree of alteration, called mean opinion score (MOS). Degradation would not be necessary perceptible according to visible and audible spectrum. Such approach is heavy to set up, as human organization, contrary to objective evaluation, as computing approach.

Objective evaluation is easily computerized. It quanties dierence between original and degraded data. A majority of metrics arranged on Figure 1.5, derivate from Euclidean distance dened by L 2 , as normalized root mean square error (nRMSE) and variation of signal noise ratio (PSNR, SNR, Λ-SNR). Considering gray scale images for example (whose pixels value ranges between 0 and 255), these metrics operate pixel by pixel, computing dierence between original and degraded version. Such an evaluation is valuable and reliable for visual appreciation.

Chapter 1. General introduction

But we express certain reservations about its use for scientic data without prior testing. As depicted by Figure 1.5, we dierentiate relative from absolute approaches. Unlike the absolute error, the relative error is related to the original values. Hence, same error made on a sample would be minor if concerns a high value, and major if the original sample is small. As will be seen later in our experimental case (cf. Subsection 4.2.1), distinction is necessary and essential. Others metrics could also be tested to evaluate the degradation of the data structure (SSIM [Structural Similarity Index Measure], Fidelity, Q criteria). Finally, ultimate purpose would consist in correlating the subjective assessment with the objective scale. This in order to guarantee objective thresholds above which the decompressed data will be subjectively acceptable and guide the suitable compression.

Lr p = C c=1 P (c)-P (c) P (c) p 1 p 1 C Lr 1 = MRE Lr ∞ = max RE L p = C c=1 |P (c) -P (c)| p 1 p 1 C L 1 = MAE 1 C L 2 2 = MSE RMSE = √ MSE nRMSE =
In conclusion to this rst chapter, a general context has been established by addressing current issues about data production in scientic elds. The data volumes generated in mass are measured in binary quantity and results from precision-resolution often over estimated. Injected at dierent steps of the scientic workow, we distinguish experimental from simulation and learning data. Reect of its over consumption, the increasing number of data centers reects excess of modern time. To highlight its consequences, we refer to its economic, energetic and ecological cost. In parallel high technologies reveal limits in data management. Compression appears as a promising solution to reduce data quantity. Commonly used with multimedia data, and assimilated by default to standard formats invading our everyday life, the compression pipeline comprises dierent components whose nature inuences its eciency. Performance is notably evaluated by various parameters. Finally, 1.3. Numerical diet for (scientific) data? 13 compression could be conservative or at renable precision to increase its gain. In the latter case, it is required to assess the quality of decompressed data in order to evaluate the compression impact.

In this work, compression inspired us to manage scientic simulation data. We restricted the scope to geosciences. Natural context and interest objects, called reservoir, are introduced in the following chapter. Such areas modeled by meshes to execute ow simulation are intensively studied.

In geosciences, the meshes are complex composite data hardly manageable, and time-consuming for simulation because of potential large dimensions.

CHAPTER 2 Reservoir modeling context

We focus on the eld of geosciences and study deep geological layers. Of great interest to the oil and gas industry, they may contain rocks called reservoirs, whose porosity concentrates energy resources.

The protability of their exploitation notably relies on simulations on numerical models. They consist in meshes integrating structural discontinuities and lled with petrophysical properties. The latter are numerical data characterized by particular dynamics that, sometimes, can be preprocessed with compandor to increase the quality of the compressed data. Our knowledge of the data and of adjacent methods allows us to dene the basis of an innovative representation that preserves the geological features, and thus to limit the impact of compression on the simulation.

Chapter 2. Reservoir modeling context

There are many scientic elds using large data volumes for dedicated applications. Among them, earth sciences are undeniably intensive users, to study large earth objects and stand therefore for domains of choice in our study.

Geosciences are split into sub disciplines because Earth is composed by several envelopes and each discipline will focus on a particular one. From the deeper to the higher, we distinguish the Earth's core, mantle, crust, biosphere-hydrosphere and atmosphere, they are studied by diverse communities from the geology to the climatology. Although resorting to dissimilar objects and models, they may use comparable methodologies. Understanding of physical phenomena occurring there bases on local eld measures feeding a simulation model made of equations system. However geosciences could refer to all these natural objects, we only focus on what happen on Earth's crust, and will use for simplication the term geosciences for this limited area.

The part of geosciences that is interested in energy production, notably the oil & gas industry, remains of strong economical interest. Hence it keeps on investigating part of the subsoil and improving elements of the workow used by geologists, geophysicists and reservoir engineers, illustrated in Figure 2.1. To create a geomodel, they use experimental data measured on the eld with dierent methods. These tools provide data at dierent format and incertitude levels: from hard data extracted by exploration wells to topographic images recovered during seismic acquisition campaigns. (Right) generated geomodel is not directly usable for simulation and requires the generation of an additional grid at lower resolutions to reduce computation time. Upscaling is operated by reservoir engineers.

Then, by simulation are predicted at certain incertitude level the production parameters. Finally, the real parameters recovered during exploitation will enrich and correct the model in a feedback loop.

Terms of the steps of the workow will be further explained. Although specic to the geosciences, we estimate a generic scientic workow would be close to this one. Throughout this workow, a complex geological object can be modeled by a geological volume mesh (GVM), mix between a volume mesh (VM) and geological mesh (GM). It corresponds in heterogeneous data composed by diverse components to represent the 3D structure and the media properties. A reservoir mesh (RM) models a particular geological area (introduced in Subsection 2.1.1). RM could be used as support to simulate the ow of liquid or/and gas that happen within.

In a quest for realism, it can be extremely detailed and dimensions can be extremely large 1 . If so, their direct simulation would be less tractable because it would be too time consuming. To prevent it, solutions exist (namely upscaling, cf. discontinuities (faults network) may strongly impact the liquids ow and so the simulation. Therefore it seems essential to identify them and preserve them throughout the workow.

A primer in geosciences for reservoir

The Earth crust could be represented by succession of horizontal rock layers. This stacking structure could be formed in diverse geological contexts. Among them we focus on conditions that give rise to the formation of fossil resources. This geological story started 20 to 350 millions years ago, with progressive deposit of sediments in particular structures called sedimentary basins.

Nowadays exploited, these structures could extend on huge surfaces, up to several hundred kilometers across, and the reservoir technical attainable set between 600 meters to 8 kilometers deep.

They originally consist of depressions in the earth's crust (most of the time) covered by water (oceans, sees, uvial systems or lakes), where the transport of sediments of variable origin stops.

The majority of sediments are detrital and come from continental erosion, but can also be organic or linked to activity of living organism (of animal as plankton shell or vegetal origin). The nature 1 up to billion of cells Chapter 2. Reservoir modeling context of sediments determines the types of rocks resulting from the sedimentation process.

Each formation (in sedimentary basin and others) could be described according to its deposit age and its lithologic composition, while its physical state would be characterized by petrophysical properties: the rock density, its porosity, its permeability and its phase saturation (water, gas, oil present in the porous system). Such information represents properties data, part of GM components, detailed next in Subsection 2.1.2. Among existing properties, we subsequently focus on two of them: porosity and permeability.

Porosity, noted φ, is the percentage of the empty volume V e on the total rock volume V tot , as described by Equation 2.1,

φ = V e V tot .
(2.1)

While permeability, noted K, evaluates the ability of rock to be crossed by liquids under a gradient of pressure noted ∇P . The Equation 2.2 is an extension of Darcy law for multi-phase ow. It computes for a liquid phase (water, oil, gas) its velocity, v, in function of its relative permeability, kr, its density, ρ and its viscosity, µ, subject to a permeable media and the gravity, g. K is measured in Darcy (homogeneous to one surface, one Darcy is equal to 9.87 × 10 -13 m 2 ). Maturation process, generation and storage of hydrocharbon Coming back to the formation of hydrocarbon resources, we know that the raw material is rich in carbon. It is mostly issued from decomposition of living matter deposed on the basin oor. To be preserved from oxygen and decarbonisation, very specic conditions need to be set up around carbon matter. This requires fast covering of organic elements by impermeable layers of mud, creating anaerobic environment. This progressively allows to bury high quantity of carbon under successive sediment formations.

At depth, under extreme pressures and temperatures, organic sediments are changed into source rock, as shown on the right of the Figure 2.2. Then starts the formation of oil and gas during cooking process. These hydrocarbons will next migrate along faults across covered formation up to porous layer called reservoir rock, whose porosity is naturally occupied by water, but replaced by oil and gas during migration. In particular stratigraphic conditions 2 if impermeable formation (named caprock) covers reservoir formation, hydrocarbons stop their progression and accumulate in this structure called trap.

Across geological time, formation features and structures are susceptible to evolve. Already complex with regard to the petrophysical structures (channels) initiated during sedimentary deposition, horizontal layers can be modied by tectonic events (compression, extension) creating folds and 2 submit of anticline structure: type of fold that is an arch-like as illustrated on right Figure 2.2.

Chapter 2. Reservoir modeling context

faults. As saw with hydrocharbon migration, these objects considerably inuence the liquid phase ow. The understanding of the eld allows us to pay attention to these details. Although such small values appear insignicant relative to the global mesh scale, it is crucial to represent for GM modeling and processing.

Considering software environments developed to manage GMs (such as Openow, Gocad, Petrel)

we notice their interest with regard to physical scales, and illustrate it by two following examples.

First, from a numerical point of view, the permeability property is commonly visualized with an exponential color scale instead of linear to focus on small values. Secondly, from structural point of view, it is possible to exaggerate the vertical scale to inate the thin geological layers to make them accessible.

Then in next section, is detailed a mesh structure and its components, commonly used by geologist to represent complex geomodel and tend to preserve geological features.

VM for reservoir modeling

VMs discretize the interior structure of 3D objects, sedimentary basin and reservoir in this study.

They partition their inner space with a set of three-dimensional elements named cells. While pyramid and triangular prism partitions exist, most of the existing VMs are composed of tetrahedral (4 faces)

or hexahedral (6 faces) elements. They are called tets or hexes (sometimes bricks), respectively. A VM composed of dierent kinds of cells, tetrahedra and hexahedra for instance, is termed hybrid.

A VM is described by the location of vertices in 3D space (geometry ) and the incidence information between cells, edges, and vertices (connectivity ). In function of the application domain, VMs also contain physical properties (petrophysical properties in geosciences) associated to vertices, edges, or cells. Note that this mesh is hybrid and unstructured, with both hexahedral and tetrahedral elements.

A non-degenerate hexahedra has 6 faces named quads, 12 edges, and 8 vertices. Depending on incidence information between cells, edges and vertices, hexahedral meshes are either unstructured or structured. The degree of an edge is the number of adjacent faces. An hexahedral mesh is unstructured if cells are placed irregularly in the volume domain, i.e., if degrees are not the same for all edges of the same nature. Unstructured meshes have an important memory footprint, as all the connectivity information must be described explicitly. However, they are well-suited to model complex volumes, Computer-aided design (CAD) models for instance, as shown in Figure 2.4.

An hexahedral mesh is structured if cells are regularly organized in the volume domain, i.e., if This pillar grid also allows to model geological collapses (or erosion surfaces), by using degenerate cells, i.e., cells with (at least) two vertices on one pillar located at the same position (see Figure 2.8 for dierent degenerate congurations).

Globally, the description of the structure represents a large part of overall mesh sizes, between one or two thirds for the GMs we studied. The remaining data consists in mesh properties attached to the cells, describing their activity or various petrophysical properties: continuous or categorical.

Mesh properties Firstly, the cell activity is a Boolean parameter notably used to deactivate cells and model irregular borders, or used to focus on an interest zone.

Secondly, there are categorical (N) and continuous (R) properties. Such properties are spatially distributed in the mesh: each cell is linked with k discrete values and n oating-point values corresponding respectively to the k categorical properties and the n continuous properties. Their distribution is determined during mesh modeling by geostatistical tools (variogram, kriging) and based on eld measurements (typically well and the seismic data). The categorical properties describe for instance the geological formation repartition (rock nature, deposit age). The discrete value is more a label which could corresponds to various parameters, not rationally quantiable or manipulable.

The continuous properties describe petrophysical properties attached to the cell, such as the porosity and the permeability, dened in Subsection 2. and require to preserve a high numerical precision. Considering for comparison an other quantity as the values of pixels of an image. For black and white version, a gray scale is commonly used and featured by a certain binary level, this means a certain number of bits, noted nbits, will be used for each pixel. For a gray scale at 256 binary levels, a pixel intensity value necessitates thus one byte to be coded (log 2 256 ⇔ 8 bit); contrary to a physical value, whose numerical accuracy uses 7 digits and requires on average up to 24 bit.

Using dynamic of reservoir property

As we plan to work on compression, we have to take into account the binary representation. We know that on practice the data is quantied due to the limited precision of computers. Given this limitation, we have chosen to exploit it by giving the advantage of precision to certain ranges and saving bits budget on others. As an example, we consider a large area with an almost constant permeability of 1000 mD, if a value were to increase by 2 mD, it would have no impact on the ow.

After a certain level of permeability, the material allows the liquid to pass through. Conversely small values will have a strong impact on the ow rate, if we add 2 mD to an initial value equal to 10 -5 mD, this could create an artifact and have an impact on the ow. This knowledge can be used to our advantage, instead of an evaluation based on absolute error, we will focus on relative scale. The logarithmic transform will allow to linearize a power law and to preserve a relative deviation.

Referring to the literature, it is known that this approach has already been used for the processing of audio data. It allows coding fewer bits because the transform makes the data more homogeneous.

This treatment, called compandor [START_REF] Li | Compressing and companding high dynamic range images with subband architectures[END_REF], noted Λ , is used to preserve whispering in a loud sound environment. Compandor echoes to expander, noted

Λ

, the inverse method to return to the original distribution.

We propose in this study an original version of the compandor noted Λ α,nbits , parametrized by α and nbits, expressed by Equations 2.4, with the original property denoted P .

Chapter 2. Reservoir modeling context

α can vary between 0 and 1. Equal to one, Λis a simple linear function (stay in absolute scale), but tends to logarithm functions with α approaching zero (moving to relative scale).

The factor 2 nbits -1 permits to convert the real values to discrete values distributed between 0 and (2 nbits -1). Thus, the choice of nbits allows us to clearly set the data precision.

Λ α,nbits (P ) = round λ α (P ) -λ α (min P ) λ α (max P ) -λ α (min P ) × (2 nbits -1) , (2.4)
with λ α (P ) = log (P + 1) if α = 0

(P +1) α -1 α if α > 0
Whereas the expander Λ α,nbits is expressed as follows:

λ α (P ) = exp (P ) -1 if α = 0 (αP + 1) 1 α -1 if α > 0 P = Λ α,nbits (Λ α,nbits (P )) = λ α ( Λ α,nbits (P ) (2 nbits -1) × (λ α (max P ) -λ α (min P )) + λ α (min P )). (2.5)
To demonstrate the compandor eect, we use a vector composed of eight values between 0 and 20 000 displayed on the rst line of the table 2.1. The data is a representative sample of the permeability property later introduced on the Figure 3.20, partly composed of small values lower than hundred and sharing the same upper limit. By applying the Λ α,nbits , we illustrate the interest of using an adapted preprocess to preserve precision of the lowest values at a comparable nbits. We start with 2 bits and incrementally increase nbits up to 4. Binary writing of a value is vertically scaled from the most signicant bits (MSB) to the least signicant bit (LSB). The nbits increment adds an LSB layer and increases the accuracy of the results by using the Λ α,nbits to return to the original distribution.

We test on the rst column a linear function by using α equal to one, and compare results applying logarithm function by using α equal to zero on the second column. Focusing on the original value 10, we observe that for a linear preprocess, the result is always equal to 0, whatever the nbits used.

In comparison, by applying Λ 0 and by increasing the nbits, result for 10 is rened from 26.10 to 15.90 to 13.00. Even considering 10 000, half of the maximal value, result from Λ 0,4 use is closer to the original than the result of using Λ 1,4 , respectively equal to 10 334.10 and 10 666.70.

Then, we consider the Λ 0 helps to better preserve the precision of the values of the sample, as regards the values of the lower half, by passing from a linear scale to a relative scale. What is an opportunity in our scientic context as in other elds dealing with similar data 3 .

By revisiting the necessary precision to obtain realistic simulations, we decide to propose an adaptive representation combining the numerical with spatial (multiscale) precision. 3 Although it does not seem obvious what data would need this type of treatment. Objectively, we could start with the data visualized using logarithmic scales.

α = 1 α = 0 Λ α,2 [ 0 0 0 0 0 0 2 3 ] [ 0 0 0 1 1 2 3 ] MSB 1 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1 LSB 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 Λ α,2
[0 0 0 0 0 0 13 333. Λ α,3

[ 0 0 0 0 0 0 4 7 ] [ 0 0 0 2 3 5 7 ]
MSB 2 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1

1 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 LSB 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 1 Λ α,3
[0 0 0 0 0 0 1142.9 20 000] [0 0 0 15.9 68.7 1179.8 20 000 20 000] Λ α,4

[ 0 0 0 0 0 1 8 15 ] [ 0 0 1 4 7 10 14 15 ] MSB 3 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1

2 0 0 0 0 0 0 0 1 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 1 0 0 0 0 1 1 1 LSB 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 Λ α,4
[0 0 0 0 0 1333.3 10 666.7 20 000] [0 0 0.9 13.0 100.7 735.8 10 334.1 20 000] Table 2.1: Λ α,nbits eect on 1D data stream, changing α value (column), as the bits number (row).

Chapter 2. Reservoir modeling context

From storage to simulation

Problematic Beyond the diculties of processing (storage, transmission, visualization) large data volumes, some GMs could serve as inputs for simulation (RMs). The process gathers diverse heterogeneous components describing structure to petrophysical properties. Each component is featured by specic dimension and dynamic range. Their combination models a geological object used to simulate the ow of liquid phases (oil, gas, water) occurring inside. However, the simulation is not directly executable on the detailed mesh, but imposes an additional step in workow by creating a supplementary data at lower resolution to run a faster simulation. Several ad-hoc methodologies already exist and reduce RM dimensions by upgridding/upscaling the grid and the attached properties, as detailed later in Subsection 3.1.1. However, this approach does not appear to be the most ecient to data processing as it multiplies the data volumes and pipes in the workow. As regards the multimedia, compression tools based on multiresolution approaches are promising options, but are rarely used for simulations.

Proposition From these observations, we propose an approach inspired by tools from various elds and our geoscience knowledge that could integrate the workow. We are thus working on an multiscale representation for GVMs dealing with components heterogeneity. For this purpose we propose a method adapted to the data dimensions and dynamic range. In addition, we plan to devote a larger share of the binary budget to a particular range, the accuracy of which is required for simulation, this by changing notably the absolute scale into relative.

The change in scale and respect of the physical laws governing our properties inuence our evaluation metrics developed to ensure the preservation of data at renable precision (spatial & numerical). Using adapted metrics we evaluate the method for dierent applications, from storage to visualization and simulation.

Figure 2.9 teasers our proposition for GVM representation. Essence of data is progressively condensed across decompositions into smaller red subpart on gure. It facilitates access to the lower resolutions and order the binary data to get the appropriated numerical precision for simulation.

Thus, we want to know if huge GVMs can be wisely represented to deal with dierent aspects taking into account the context.

Plan Our work is divided into two chapters relating to our main contributions. First, in Chapter 3, we introduce a suitable representation for our GVMs, i.e. capable of arranging the data while reducing its binary quantity for compression purpose. Secondly, in Chapter 4, the representation is integrated into a simulation workow to evaluate the impact of the RM generated at renable precision (still considering spatial & numerical precision).

To begin, we review in Science the existing representations used to manipulate volume data, starting with geosciences and then broadening the spectrum (cf. Section 3.1). By combining the strengths of all methods, we dene a dream scalable representation for our complex GVMs (cf. Subsection 3.1.4). To lay the foundations of our representation we use HexaShrink, a multiscale decomposition tool for GVMs, preliminary developed for visual purpose, introduced in Section 3.2.

We rst apply HexaShrink on a benchmark composed by eight GVMs (cf. Section 3.3). Their lower resolutions are visually compared to the results obtained with a popular geomodeller. Then to evaluate compression performances of the scalable representation we associate HexaShrink with generics encoders in a conservative approach (cf. Section 3.4) and nally with an evolved progressive coder, named zerotree, to evolve toward renable approach (cf. Section 3.5). The quest for the suitable rened precision is our objective to improve compression but not only. In Chapter 4, we also assess the precision within simulation workow.

This second contribution starts with a state of the art on renable compression in simulation (cf.

Section 4.1). In this section, we detail some of the tools that will be later used for a comparative study. In Subsection 4.2.1 we model a complete simulation workow to evaluate our representation rst in an upscaling purpose (cf. Subsection 4.2.2). Secondly the renable precision is applied on mesh properties preprocessed or not by compandor to increase the precision on specic ranges (cf. Subsection 4.2.3). To validate this processing, innovative metrics are suggested based on the eld knowledge and expectations. By combination with standard objective metrics we verify their coherency. We thus seek to identify the parameters to evaluate the precision of data suitable for simulation. We exploit on a last part some data features to improve the process (cf. Section 4.3).

CHAPTER 3

HexaShrink, a multiscale representation for geological volume meshes. Application to multiresolution rendering and storage "Compression... What else?"

Laurent Duval, 2017

Scientic community is interested in various methods for handling large volume data. Regarding the specic eld of geosciences, we rst mention the upscaling/upgridding techniques, daily used in simulation for "simplifying" geological volume meshes. Then, we make an overview of generic tools developed for compressing volume meshes. Finally, we approach the most popular techniques existing for compressing scientic data in general. By gathering the bright sides of each method, we present HexaShrink, a scalable representation dedicated for the geological volume containing attributes and discontinuities. We show that HexaShrink provides nice renderings at dierent levels of detail, but can be also integrated in dierent lossless-to lossy compression workows, to facilitate the storage of geological volume meshes.

29

Chapter 3. HexaShrink, from mesh representation to coding

A variety of representations for volume scientic data

In the rst place, this section presents three ways that we explored to put forward a suitable representation for the GMs, and particularly the RMs. Subsection 3.1.1 gives an overview of the upscaling/upgridding methods, well established in geosciences for years in the context of simulation.

In a nutshell, it consists in generating a coarse mesh from a given highly resoluted mesh, and to use it to drive simulation. By decreasing the number of cells, simulation time is ineluctably reduced.

Second, Subsection 3.1.2 overviews the multitude of compression methods dealing with VMs. The main purpose of these methods is a compact storage, in scientic elds generally not concerned by simulation. Third, Subsection 3.1.3 is devoted to newer compression methods developed for scientic data with oating-point accuracy specically. Finally, Subsection 3.1.4 summarizes the advantages and the limitations of these three approaches in our context. We also discuss the required features of an "ideal" representation for GMs. ), but this could not be performed on the initial mesh.

To reduce the computation times, a coarse grid is generated from the initial data, by upscaling the properties and upgridding the RM. These scale changes involve a merge of connected cells and their attached continuous properties, reducing signicantly the number of cells. The porosity is often upscaled by a simple bulk volume-weighted arithmetic average, while permeability upscaling is dicult, because the property is not additive. A generic approach consists in using arithmetic and harmonic averages for the horizontal and vertical permeability, respectively. This method is called static, performing arithmetic operations on property data. It could be sucient for simpler RMs. [START_REF] Li | Optimal uplayering for scaleup of multimillion-cell geologic models[END_REF] study the spatial distribution of permeability, properties, and facies to determine the optimal merging, and preserve details on the way of the main ow. However other dynamic methods, based on the ow study are recommended for more complex meshes. To name one [START_REF] Durlofsky | Scaleup in the near-well region[END_REF] uses a ner resolution in high uid velocity area, identifying it by the resolution of local well-driven ow problems. A large number of upscaling methods exists, provided in detail by reviews [START_REF] Qi | Major challenges for reservoir upscaling[END_REF]. The upgridding or grid coarsening is related to the upscaling result, adapting the mesh structure by merging the cells whose the properties have been upscaled. Quality of upscaling and upgridding can be assessed according to simulation results, as latter explained in Subsection 4.2.1, or by other metrics taken on upscaled data [START_REF] Preux | About the use of quality indicators to reduce information loss when performing upscaling[END_REF]. The main challenge is to choose the more ecient for a given model, the one that will preserve result accuracy, while allowing a signicant reduction of the cell number and consequently of computing time.

A wavelet-based upgridding/upscaling technique has been proposed by [START_REF] Mehrabi | Coarsening of heterogeneous media: Application of wavelets[END_REF] The irregular grids in geosciences are well known, and represent an elegant way to model reservoir complexities, notably the faults. Nevertheless, the cell indexing is not possible, so the resolution of linearized equation systems for ow simulation are much more complicated.

Discussion These wavelet-based methods do not save the details removed during analysis to get RMs of lower resolution, which could permit a perfect reconstruction by reversing the process (i.e., by synthesis). This pushes the users to store the initial data in parallel.

During the review made on geosciences and the simulation elds, we also note the lack of details on the wavelets used. This contrasts with the current practices in multimedia, where the precise nature and the characteristics of the wavelets are given. We thus consider that, depending on the nature and the regularity of the dierent elements of our grids, a more rigorous choice of multiscale representations can be benecial.

In a more global way, the upscaling techniques essentially serve the interest of the simulation providing an additional lighter RM, promising faster executions and coherent results without taking into account the consequence for storage.

Compression of VMs

VMs are used in many research domains. These are data that can be particularly massive, and often generated in large quantities. Therefore many projects have focused on how to compress them.

The objective common to all the compression techniques presented hereinfater is to reach the best compression ratios.

The basic principle and the most straightforward technique to encode a VM is to use an indexed data structure: the list of all the vertex coordinates (three oating-point values, which amounts to 96 bits par vertex), followed by their connectivity. The connectivity is dened cell after cell, each cell being dened by the set of indexes of the adjacent vertices (8 integers per hex). They thus only provide estimates of an actual compression performance. when meshes are unstructured, the most frequent technique performs a traversal of mesh elements, and describes the incidence congurations with a reduced list of symbols. These symbols are then entropy coded. When meshes are structured, the connectivity is implicit, reducing its cost to zero.

For such meshes, the only additional information to encode are geometrical discontinuities describing faults and gaps.

The basic tools previously presented can be implemented on the ontological structure of meshes, and improved in many dierent ways. Their combination, with the assistance of advanced compression techniques, permits more ecient tetrahedral or hexahedral mesh coding. Previously proposed algorithms are presented below, classied into two categories: single-rate and progressive/multiresolution.

Single-rate mesh compression They lead to a compact mesh representation, most of the time driven by ecient connectivity encoding. The rst method for tetrahedral meshes, Grow & Fold, was presented by [START_REF] Szymczak | Grow & fold: Compressing the connectivity of tetrahedral meshes[END_REF] at the end of the nineties. It is an extension of EdgeBreaker [START_REF] Rossignac | Edgebreaker: Connectivity compression for triangle meshes[END_REF] developed for triangle meshes. The method consists in building a tetrahedral spanning tree from a root tetrahedron. The traversal is arbitrary among the three neighboring tets (cf. Subsection 2.1.2) of the cell currently processed, and 3 bits are needed to encode each cell. The resulting spanning tree does not retain the same topology as the original mesh, because some vertices are replicated during the traversal. Fold and glue techniques are thus needed during encoding to restore the original mesh from the tetrahedron tree. The additional cost is 4 bits, leading to a total cost of 7 bits per tetrahedron.

The cut-border initiated by Gumhold and Straÿer (1998) was adapted to tetrahedral meshes [START_REF] Gumhold | Tetrahedral mesh compression with the cut-border machine[END_REF]. It denotes the frontier between tetrahedra already encoded and those to encode. At each iteration, either a triangle or an adjacent tetrahedron is added to the cut-border.

In this case, if the added vertex is not already in the cut-border, this latter is included by a connect operation, and is given a local index. As the indexing is done locally, the integers to encode are very small, leading to a compact connectivity representation. In addition, two methods are proposed to encode geometry and associated properties, based on prediction and entropy coding.

This method yields good bit rates (2.40 bits per tetrahedron for connectivity) for usual meshes, handles non-manifold borders, but worst-cases severely impact bitrates and runtimes (which tend to be quadratic). ). In the latter, for instance, the compressor does not require the knowledge of the full list of vertices and cells before encoding. The compressor starts encoding the mesh as soon Chapter 3. HexaShrink, from mesh representation to coding as the rst hexahedron and its eight adjacent vertices have been read. For a given hexahedron: i) its face-adjacency is rst encoded in function of its conguration with hexahedra already processed;

ii) positions of vertices that are referenced for the rst time are predicted (spectral prediction from adjacent cells); iii) prediction errors are encoded; iv) data structures relative to vertices, becoming useless (because their incidence has been entirely described) are nally removed from memory. Compared to other single rate techniques, streaming tends to achieve similar compression performances for geometry, but poorer performances for connectivity.

Progressive/multiresolution mesh compression Progressive algorithms (also called scalable or multiresolution) enable the original meshes to be represented and compressed at successive LODs (levels of details). The main advantage is that it is not necessary to decompress a mesh entirely before vizualising it. A coarse approximation of the mesh (also known as its lowest resolution) is rst decompressed and displayed. Then this coarse mesh is updated with the successive LOD (termed higher resolutions) that are decompressed progressively. While they cannot achieve yet compression performance of single-rate algorithms, progressive algorithms are popular because they enable LOD, and also adaptive transmission and displaying, in function of user constraints (network, bandwidth, screen resolution. . . ). [START_REF] Pajarola | Implant sprays: Compression of progressive tetrahedral mesh connectivity[END_REF] are the rst to propose in 1999 a progressive algorithm dedicated to VM compression. This work is inspired by a simplication technique for tetrahedral meshes [START_REF] Staadt | Progressive tetrahedralizations[END_REF]. It simplies a given tetrahedral mesh progressively, by using successive edge collapses [START_REF] Hoppe | Mesh optimization[END_REF]. Each time an edge is collapsed, its adjacent cells are removed, and all the information required to reverse this operation is stored: index of the vertex to split, and the set of incident faces to cut". Thus, during decompression, the LODs can be also recovered iteratively, by using the stored data describing vertex splits. During coding, an edge is selected such as its collapse leads to the minimal error, with respect to specic cost functions. This algorithm gives a bitrate inferior to 6 bits per tetrahedron (for connectivity only).

In 2003, [START_REF] Danovaro | Multiresolution tetrahedral meshes: An analysis and a comparison[END_REF] propose two progressive representations based on a decomposition of a eld domain into tetrahedral cells. The rst is based on vertex splits, as the previous method, the second is based on tetrahedron bisections. This operation consists in subdividing a tetrahedron into two tetrahedra by adding a vertex in the middle of its longest edge. Unlike with vertex splits, the representation based on tetrahedron bisections is obtained by following a coarse-to-ne approach, i.e., by applying successive bisections to an initial coarse mesh. Also, this representation only needs to encode the dierence vectors between the vertices added by bisections and theirs real positions.

This representation is thus more compact, as the mesh topology does not need to be encoded, but it only deals with structured meshes.

Compression of scientic data

Ongoing advances in computer science allows performing simulations on ever larger and detailed models, while their management becomes ever more dicult. This subsection present solutions for compressing 3D experimental and simulation data. By using various approaches: data reduction As always in compression, two categories can be identied: lossless and lossy. 3.1.4 How taking advantage of these three techniques for a dream representation of GVM?

We now discuss the main features borrowed to the methods of upscaling/upgrading, VM compression and scientic data compression. Then, we set the bases of a new representation adapted to our GVMs. Concerning our dream representation, we can gather the following key observations:

Our representation has to be based on an embedded multiscale decomposition scheme, which oers obvious advantages for visualization, storage, but also simulation. One of the avantages is that the decomposed data could not take up more space than the initial data, unlike the upscaled/upgraded data.

Despite their high performance in term of compression, the progressive methods developed for VMs are not appropriate for our data. Indeed, the GMs contain fault networks, involving local geometrical discontinuities that are not managed by these generic compression methods.

Therefore our multiscale mesh representation must take into account the specic morphological features of the initial data, and must preserve them as much as possible across the resolutions.

The continuous properties are 3D oating values, sometimes of high numerical precision, which is not always necessary for visualization or simulation. Thus our representation must propose an encoding of these properties at renable precision, according to the application and the user needs. To assess the eciency of our encoder at renable precision, popular compression tools for scientic data such as ZFP or SZ will serve as reference for comparison study.

Particular orientations are observable for properties, that can be explained by their geological origin. This anisotropy needs to be preserved at various scales, preserving by this way the coherency of the mesh in a multiresolution structure.

Chapter 3. HexaShrink, from mesh representation to coding

During a previous collaboration between IFP Energies nouvelles and I3S, some of these observations have been considered to implement the core of HexaShrink [START_REF] Peyrot | H)exashrink: Multiresolution compression of large structured hexahedral meshes with discontinuities in geosciences[END_REF][START_REF] Peyrot | HexaShrink, an exact scalable framework for hexahedral meshes with attributes and discontinuities: multiresolution rendering and storage of geoscience models[END_REF].

HexaShrink is a multiscale decomposition scheme for hexahedral meshes from geosciences. It combines several wavelet transforms adapted to the specic features of these meshes (detailed in next section). Firstly developed for multiscale renderings, the HexaShrink structure is the starting point of our current research, focused on compression and its impact on simulation results.

Description of HexaShrink

This section introduces HexaShrink, the mesh hierarchical structure on which our work has been based. Some paragraphs of this section are extracted from [START_REF] Peyrot | HexaShrink, an exact scalable framework for hexahedral meshes with attributes and discontinuities: multiresolution rendering and storage of geoscience models[END_REF].

MRA or multiscale approximation can be interpreted as a decomposition of data at dierent resolutions, LOD or scales, through a recursive analysis process. It is called exact, reversible or invertible when a synthesis scheme can retrieve the original data. Inter-scale relationships [START_REF] Chaux | Noise covariance properties in dual-tree wavelet decompositions[END_REF][START_REF] Chaux | A nonlinear Stein based estimator for multichannel image denoising[END_REF] often yield sparsication or increased compressibility on suciently regular datasets. In discrete domains, each analysis stage transforms a set of values (continuous or categorical, in one or several dimensions), denoted by S 0 . The resulting representation consists in one subset S -1 that approximates the original signal at a lower resolution, plus one subset of details D -1 , or a combination thereof. The latter represents information missing in the approximation S -1 . Depending on the MRA scheme, the lower resolution S -1 may represent a coarsening or low frequencies of the original samples, or an upscaling in geosciences (cf. Subsection 3.1.1). The subset D -1 represents renement, fast variation or high-frequency details removed from S 0 . We consider here exact systems, allowing the perfect recovery of S 0 from a combination of subsets S -1 and D -1 .

Hence, a similar analysis stage can be applied iteratively, and perfectly again, to the lower resolution S -1 , in a so-called pyramid scheme. Thus, with the non-positive extremum decomposition level L, and indices 0 ≥ l ≥ L, after an |L|-level multiresolution decomposition, the input set S 0 is now decomposed and represented by the subset S L a (very) coarse approximation of S and |L| subsets of details D L , . . . , D l , . . . , D -1 , representing information missing between each two consecutive approximations.

To get an idea of what a MRA will generate from VMs, Figure 3.2 gives an overview of a decomposition obtained with HexaShrink on one of our GMs. To obtain this decomposition, we consider in the following four dierent MRA avors, all called wavelets for simplicity. They stem from iterated, (rounded) linear or non-linear combinations of coecients, as well as separable (applied separately in 1D on each direction) or non-separable ones. Without going into technicalities here (cf. Subsection 3.2.1), computations are performed using the lifting scheme. It suces to mention that lifting uses complementary interleaved grids of values, often indexed with odd and even indices.

Values on one grid are usually predicted (approximations) and updated (details) from the others.

The main interests reside in reduced computational load, in-place computations and the possibility to maintain exact integer precision, using for instance only dyadic-rational coecients (written as Hence, if a VM at resolution l is composed of

m/2 n , (m, n) ∈ Z × N)
C l i × C l j × C l k cells, C l
i being the number of cells in direction i, the VM of lower resolution will be of dimension

C l i 2 × C l j 2 × C l k 2
, to take into account non-power-of-two sized grids. As several digital attributes are associated to each cell (geometry, continuous or categorical properties), dierent types of MRA are performed separately on the dierent variables dening these properties, as explained in the following sections.

Multiresolution scheme for geometry

Standard linear MRA schemes rely on smoothing or averaging and dierence lters for approximations and details, respectively. To preserve coherency of representation of geometrical discontinuities whatever the resolution a special care is taken to avoid excessive smoothing, while at the same time allowing the reverse synthesis. As the pillar grid format is used (cf. Horizontal 2D morphological wavelet transform The fault segmentation guides the MRA to preserve faults, as much as possible, all over the decomposition process. The fault conguration of 4 associated nodes at resolution l is used to predict the extension of the downsampled fault structure at resolution l -1.

This horizontal prediction is based on the logical function OR (∨), computed on each side of each group of 4 nodes. For instance, a resulting fault node conguration contains a west axis if the fault congurations of the 2 left nodes contain at least 1 west axis, as illustrated in Figure 3.9.

By repeating the procedure for each axis of each resulting node, fault node congurations at lower resolutions are fully predicted. This non-linear and peculiar choice is meant to maintain a directional Chapter 3. HexaShrink, from mesh representation to coding avor of orientated faults for ows; other choices could be devised, depending on physical rules and geological intuitions.

Finally, from this prediction, the node whose conguration minimizes its distance with the predicted one, corresponds to the aforementioned approximation coecient, which will be part of the novel Z matrix at lower resolution l -1. The same procedure can be applied recursively until the wanted resolution. Rounded linear 1D wavelet transform This 1D wavelet transform is applied on the output of the above horizontal 2D morphological wavelet, to analyze the z coordinates of the 4 sets of vertices BDR, FDR, BDL and FDL separately, along each selected pillar. The HexaShrink multiresolution here acts on z-locations, decomposing them at each scale location z l into a subsampled pillar coordinate z l-1 and its associated detail d l-1 . By geomodel construction, coordinate behavior along the pillars is expected to be relatively smooth. This entails the use of a modied, longer spline wavelet. The latter can be termed LeGall (Le Gall and Tabatabai, 1988), or CDF 5/3 (after Cohen, Daubechies and Feauveau [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF])), or biorthogonal 2.2 from its vanishing moments.

The lifting analysis operations Prediction and Update are depicted by Figure 3.6. To retrieve respectively the sets of details d l and the z l coordinates at resolution l -1 from scale l, the following equations are used (∀n ∈ N):

d l-1 [n] = z l [2n + 1] - z l [2n] + z l [2n + 2] 2 , (3.1) z l-1 [n] = z l [2n + 0] + d l-1 [n -1] + d l-1 [n] 4 , (3.2)
where both dyadic integers and rounding are evident. With rounding, lifting schemes can thus manage integer-to-integer transformations [START_REF] Calderbank | Wavelet transforms that map integers to integers[END_REF]. For synthesis, to reconstruct

Description of H S

47 resolution l from resolution l -1, we only have to reverse the order and the sign of the equations:

z l [2n] = z l-1 [n] - d l-1 [n -1] + d l-1 [n] 4 , (3.3)
z l [2n + 1] = d l-1 [n] + z l [2n] + z l [2n + 2] 2 .
(3.4)

Managing externalities: borders and boundaries A pertinent multiresolution on complex meshes requires to cope with externalities that may hamper their handling: oor and ceil borders and outer boundaries (Figure 3.5). First, to keep borders unchanged from the original mesh, throughout all resolutions, the following constraints must be met:

z l-1 [0] = z l [0] , (3.5)
z l-1 [n l-1 k -1] = z l [n l k -1] . (3.6)
Both constraints can be fullled if one satises the following conditions:

Floor border condition to meet (3.5):

d l-1 [-1] = -d l-1 [0] , (3.7)
Ceil border condition to meet (3.6):

d l-1 [n l-1 k -1] = -d l-1 [n l-1 k -2] , (n l k odd) (3.8) d l-1 [n l-1 k -1] = -d l-1 [n l-1 k -2] (n l k even) +4z l [n l k -1] -4z l [n l k -2] .
(3.9)

To complete the MRA of the geometry, the same rounded 1D wavelet is also applied to the sets of oor and ceil nodes of the initial VM, to get the x and y coordinates of the extremities of the remaining pillars at the lower resolution.

Second, the Actnum eld should also be considered to lessen mesh boundary artifacts. Indeed, severe disturbances may appear at lower resolutions if not wisely processed during analysis, as shown in Figure 3.10. A cell is deemed active if and only if its 8 adjacent vertices are active at the resolution l. During our study, we found that one vertex at resolution l -1 could be considered active if and only if its parent vertices selected by the morphological wavelet at resolution l (cf. Subsection 3.2.1) are active. So, a cell at resolution l -1 is considered active if and only if its 8 × 2 corresponding parent vertices are active at resolution l.

Multiresolution scheme for properties

Continuous properties Once the geometry is coded, one can focus on associated continuous properties. For scalar ones, a value p i ∈ R is associated to each cell i in the mesh. Consistently with Chapter 3. HexaShrink, from mesh representation to coding

The geological meshes are produced by petroleum companies during expensive drilling projects and seismic programs. These meshes are realistic, and carefully treasured by companies. Collecting geological meshes thus is a complicated quest, which explains the limited number of meshes in our benchmark. However, this dataset aims to be representative of the variety of meshes found in geosciences.

Our eight meshes are illustrated by Figure 3.11, and their main features are indicated in Table 3.1. They have various geometries (from smooth to fractured), and dierent dimensions (from small to large). To have an order of magnitude, the largest mesh (mesh#8) contains almost 380 times more cells than the smallest one (mesh#3). Also, several meshes contain continuous and/or categorical properties. Mesh#1 is decomposed to the lowest possible resolution (Figure 3.12, bottom), which is not very relevant from a geologist perspective. However, while all properties are almost constant, the lower arch corresponding to an anticlinal on the mesh at original resolution remains perceptible on the nal Lego brick resolution. Looking at the porosity property (middle column), one observes how the values are progressively homogenized on coarser hexes. Concerning the rock type (last column), one observes that the modelet scheme tends to locally maintain predominant categories resolution after resolution, which is very satisfactory.

Figure 3.13 shows our largest mesh, mesh#8, that contains an isolated fault on the left side (the Figure 3.14 confronts meshes mesh#5 and mesh#7 downsampled at power-of-two resolutions with HexaShrink and with the geomodeller SKUA-GOCAD. SKUA-GOCAD, but also PE-TREL, are two geomodellers frequently used in geosciences to handle geological objects and to generate meshes for ow simulation. They include in particular upscaling/upgridding functions. We recall that upscaling/upgridding functions convert the properties and the geometry of a given mesh in a non-reversible manner to obtain a simplied version for ow simulation. Such functions are usually exible yet often ad-hoc. We can see that HexaShrink tends to better preserve faults (colored in red). Figures emphasize an improved preservation of mesh borders, with an ecient management of Actnum throughout resolutions. Some artifacts may appear with SKUA-GOCAD's upgridding, which are automatically averted by HexaShrink, leading to nicer meshes at low resolution. As a summary, HexaShrink, while being fully reversible at dyadic scales only, eciently and automatically manages structural discontinuities in the VMs. It may provide an interesting complement to existing irreversible upscaling/upgridding proposed by several geomodellers.

Conservative compression workow for GMs

We just showed that HexaShrink is a powerful tool for structuring the GVMs into in pyramidal models, emphasizing the continuous zones, and rendering them sparser by transformation. It permits to display nice simplied versions, faithful to the original data. Additionally, HexaShrink should be also suitable for compression. Therefore, we will now study its potential compression performance, by combining it with an appropriate encoder. Our rst experimentation is to combine HexaShrink with generic lossless encoders, to propose a conservative compression workow.

We retain three generic lossless encoders: gzip (1992), bzip2 (1996), and LZMA (1998). They are used by a general audience to compress diversied data. gzip has been proposed by Adler and Gailly.

The rst version was released in 1992 while optimized versions are still updated on zlib library. Its Chapter 3. HexaShrink, from mesh representation to coding Geometry.

Porosity. Rock type.

Original mesh#1.

Resolution -1.

Resolution -2.

Resolution -3.

Resolution -4.

Figure 3.12: Original mesh#1, its attributes, and four levels of resolution generated with Hexa-Shrink.

Conservative compression workflow for GMs

Porosity. Rock type.

Original mesh#7.

Resolution -1.

Resolution -2.

Resolution -3.

Resolution -4. The method yields better compression ratios with slower speed performance. Again four years later in 1998 LZMA (Lempel-Ziv-Markov chain algorithm) makes good compression performance, even better than both concurrents (gzip, bzip2) using a dictionary compression process. Its memory cost as its execution time are however higher. Such preliminary tools are not dedicated to scientic data unlike the following, but they illustrates the variety of approaches, and the complexity to develop an ideal tool for all applications and needs. with HexaShrink, and then encoded. The main objective is to conrm the contribution of Hexa-Shrink in the size reduction whatever the encoder used, but also to determine the most ecient for a further implementation on HexaShrink.

Conservative compression workflow for GMs

Coding performance

We focus our analysis on the mesh#5, a quite complicated and faulty model. Without decomposition, we observe variable CRs according to the coder: from 2.46 for gzip to 3.33 for LZMA. The similar observation is made for the other meshes, except for the mesh#7, in which gzip is slightly more ecient than bzip2. But whatever the mesh, LZMA provides the best compression results.

With only one level of decomposition, we can see a signicant improvement. Always considering the mesh#5, the CR obtained with LZMA increases up to 3.71. A growing gain up to 3.81 is obtained using a second decomposition, yet additional decompositions only provide a marginal improvement.

The HexaShrink decomposition improves the eciency of the coders. This observation could be generalized on the global Benchmark, yet LZMA performance decreased on the mesh#8 using

HexaShrink. This might be explained by certain mesh modeling practice with a geomodeller: horizons of the formation layers are preliminary modeled, before incorporating it into the mesh.

The user rst determines for each layer the number of cells in height, and then the thickness of the layers are computed by interpolation. Hence, locally, an observed value v may arise from a scale s and an oset o relationship v = s × m + o, on some integer index m, easier to compress than v. LZMA's superior capability owes to its capacity to capture complex models of byte patterns.

By contrast, with a wavelet decomposition, the ane relationship is poorly captured throughout approximations, due to the rounding in wavelet lifting (cf. Subsection 3.2.1). Hence, multiscale decompositions may slightly reduce the raw compression performance for meshes presenting initial numerical format artifacts, or illusory oating-point precision. This however does not hamper the usability of HexaShrink for storage and visualization, as the direct access to a hierarchy of resolutions respecting discontinuities is granted, while already providing impressive compression rates of about 89, superior to most results for the others meshes.

Speed performance

The execution time is also a signicant issue, a delay could indeed become critical depending on the application. In the present situation, we may resort to asymmetrical compression-decompression schemes, sometimes termed compress once, decompress many. The encoding can take more time because it is performed only once, to optimize storage and transfer. Rather, decoding should be faster to be performed several times, on demand. The speed performance of the conservative compression workow are highly dependent on the mesh, its complexity, its dimensions and the level of decomposition. For a baseline evaluation, a Java implementation was run on a laptop with Intel Core i7-6820HQ CPU @ 2.70 GHz processor and 16 GB RAM. Each mesh was compressed to the maximum level, and decompressed, twelve times. As the outcomes were relatively stable, they were averaged. of the encoders used. Timings of the analysis and synthesis are also presented. The synthesis stage (i.e., the inverse transform of HexaShrink) is faster than the analysis, up to 4 times for instance for small meshes such as mesh#4, while this dierence reduces with larger meshes. Of the three generic coders, gzip is denitely the fastest in general for coding and decoding, followed by bzip2 and LZMA. LZMA necessitates more times to code the data, but is faster than the bzip2 to decode the largest meshes as mesh#8. This observation supports our previous expectation, promoting a fast decompression to be performed several times against a slower compression, performed only once. Moreover, the LZMA speed performance conrms its eciency already demonstrated with the compression rates.

In-depth analysis of the coding performance

As previously explained, a GM is a composite object made of a geometrical structure and continuous, categorical properties. The components are separately treated by an adapted wavelet, generating for each component a hierarchical structure.

To complement the global coding performance of the conservative compression workow presented in Subsection 3.4. 32 bits or 8 bits, corresponding to standard formats. Applied on the raw data, LZMA eciently reduces the number of bits whatever the considered component: an altitude value (Zcorn), coded on 64 bits only needs 25.8 bits through LZMA. One rst HexaShrink decomposition allows reducing it by 2 bits in more, and a second by 0.4 bits. As previously shown with the global mesh results (the gains become modest after two decomposition levels). The same goes for the coordinates values (Pillar), both constitute indeed the geometry, and wavelet methods have already proven to be very ecient on this kind of data.

Meanwhile, the decomposition applied on the continuous property does not have the expected impact. LZMA signicantly reduces the native 64 bits needed to encode a cell label (up to 6.28), but one HexaShrink decomposition slightly increases (up to 8.46). This limited eect could be explained by a high magnitude and an over numerical precision required by the simulation. Besides the lossless generic encoders do not eciently use the inter-scale redundancy in the multiscale structure generated by HexaShrink. Such observations encourage us to test more evolved coders, such as the zerotree coders. This concept will be developed in the next subsection. Also, we observe that the transformation employed currently in HexaShrink for the cell activity is not the most appropriate. The cell activity consisting in binary values, we believe that a binary wavelet [START_REF] Pigeon | Binary pseudowavelets and applications to bilevel image processing[END_REF]) should be more performant. It will be tested in the future.

Finally, with this benchmark study, we veried the relevance of using HexaShrink in a compression workow. The global performance, as the execution time, obtained on the global meshes are meaningful and promising for our future work. However, the in-depth analysis for each mesh component reveals a disparity according to the type of property, for the continuous ones in particular. More evolved encoders could provide an adapted solution for both, exploiting the correlation

Compression workflow at refinable precision for GMs

59 inter subbands and oering performant compression at renable precision.

Compression workow at renable precision for GMs

The HexaShrink methodology yields a multiscale representation of each mesh component. As illustrated by Figures 3.12 and 3.13, the renderings at dierent scales exhibit an evident visual correlation. This inter subband correlation can be employed to improve the coding performance, with the help of a coder at renable binary precision.

How coding at renable precision?

The notion of coding/decoding (or compression/decompression) at renable precision is illustrated by Figure 3.17. In a nutshell, let us consider a 1D vector containing eight integer values (top-left), from zero to 15. Typically, the binary cost of each value is 4. Each original value can thus be progressively encoded, from the Most Signicant Bit (MSB O ) to the Least Signicant one (LSB O ) plus a sign bit, as shown in the left part of the Figure 3.17 (the index O refers to the original data).

During the decoding stage, we can reconstruct the original vector at renable precision by using only a limited number of bit planes, starting from the MSB O (plus the sign bits). On the left-part of Figure 3.17, we can also see the vector reconstructed with only three bit planes (in that case, LSB O = 1), and with only two bit planes (in that case, LSB O = 1). The lower the number of bit planes, the less precise the reconstructed vector. This principle can be also applied on transformed or decomposed data. In the right part of Figure 3.17, we show for example the original vector on which a simple |2|-level S+P transform (for sequential/prediction, cf. Said and Pearlman (1993),) has been applied to. This transform, that can be considered as a non-linear avatar of the Haar wavelet [START_REF] Haar | Zur Theorie der orthogalen Funktionen Systeme[END_REF][START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF], yields several subsets of values: one subset of approximation coecients, and several subsets of details, two in our case. If s l [k] denotes a sequence of approximation coecients at a given level l, the S+P transform yields two half-length subsequences:

s l+1 [k] = s l [2k] + s l [2k + 1] 2 , (3.10) d l+1 [k] = s l [2k] -s l [2k + 1] , (3.11) with d l [k]
denoting detail coecients at level l.

The resulting multiscale data (shown at top-right, the red bars separating the dierent subsets) can be also fully encoded with 4 four bit planes, and also decoded partially with a limited number of bit planes, before reconstruction. In that case, we can also reconstruct at renable precision.

That being said, some remarks can be done, to emphasize the relevance of using transformed data. The absolute values of the combined multiscale data are globally smaller in amplitude, as the transform captures data regularity at dierent scales. As a consequence, only one symbol "1" remains in the MSB M , and for an approximation coecient. At the opposite, the symbols "1" representing the details in the multiscale data are mainly concentrated in the LSB M , which allows reducing Chapter 3. HexaShrink, from mesh representation to coding

Original Multiscale data 0 0 1 4 7 10 14 15

data 1 11|-2 -6| -3 -3 -1 Sign 1 1 1 1 1 1 1 1 Sign 1 1 | 0 0 | 0 0 0 Most Signicant Bit 3 0 0 0 0 0 1 1 1 Most Signicant Bit 3 0 1 | 0 0 | 0 0 0 (MSB O ) 2 0 0 0 1 1 0 1 1 (MSB M ) 2 0 0 | 0 1 | 0 0 0 1 0 0 0 0 1 1 1 1 1 0 1 | 1 1 | 1 1 0 Least Signicant Bit 0 0 0 1 0 1 0 0 1 Least Signicant Bit 0 1 1 | 0 0 | 1 1 1
Result 0 0 1 4 7 10 14 15 Result 0 0 1 4 10 14 15

Sign 1 1 1 1 1 1 1 1 Sign 1 1 | 0 0 | 0 0 0 MSB O 3 0 0 0 0 0 1 1 1 MSB M 3 0 1 | 0 0 | 0 0 0 2 0 0 0 1 1 0 1 1 2 0 0 | 0 1 | 0 0 0 LSB O 1 0 0 0 0 1 1 1 1 LSB M 1 0 1 | 1 1 | 1 1 0
Result 0 0 0 4 6 10 14 14 Result 0 0 1 3 8 12 12

Sign 1 1 1 1 1 1 1 1 Sign 1 1 | 0 0 | 0 0 0 MSB O 3 0 0 0 0 0 1 1 1 MSB M 3 0 1 | 0 0 | 0 0 0 LSB O 2 0 0 0 1 1 0 1 1 LSB M 2

Compression workflow at refinable precision for GMs

61 the global binary cost. To illustrate this very interesting feature in a context of compression, let us simply calculate the total number of bits required to perfectly retrieve the original vector from the original bit planes and from the bit planes representing the multiscale data. By neglecting the sign bit, 21 and 16 bits are respectively necessary, which attests that a transformed/hierarchical data can facilitate the encoding. Moreover, the same phenomenon occurs when the reconstruction is done at renable precision, as shown in the middle and bottom parts of Figure 3.17.

This methodology can be generalized to data of higher dimensions, using dependencies between the subsets of dierent levels. Since approximations shrink by a factor of two in each dimension (cf.

Equation 3.10), 2 d blocks of values at level l relate to one value at level l + 1. This results in binary trees, whose the inheritance between the subbands are depicted in Figure 3.18 for 2D and 3D data.

Principle of zerotree coding

When multiscale data are parsed from most to least signicant bits according to the binary trees, a signicant quantity of 0 is very often observed. In that case, the binary trees are called zerotrees (ZT). They have been popularized for image compression, early with EZW (Embedded Zerotree Wavelet, [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coecients[END_REF])), later with SPIHT (Set Partionning In Hierarchical Trees or SPIHT, Said and Pearlman (1996)).

The principle of ZT coding is to replace a tree of zeros with a single symbol to reduce the information quantity. To succeed, the root node of the tree should be identied. It corresponds to the lowest scale, while related coecients in higher frequency subbands draw the tree branches.

After an appropriate wavelet decomposition, structured data become parsimonious, generating sub-

Original 5MSB M 4MSB M 3MSB M 2MSB M Res. -1
Res. -2 tribution is 0 to 0.5 for porosity, and 0.0007 to 20 000 mD for permeability, except with nearshore 0 , whose the minimal permeability value is higher (4.7 mD).

The Tarbert formation inspired us three models, named hereinafter nearshore 0 , nearshore 1 and nearshore a (the index a indicates an anisotropic distribution along one direction). The Tarbert formation conditioned one model, named hereinafter fluvial. The four resulting datasets are shown in Figure 3.21. They mainly dier in their spatial distribution: the three nearshore environments present smooth variations, while the uvial environment exhibits sharper contrasts, with distinctive heterogeneous geological objects. This discrepancy between environments primarily allows a wide range of petrophysical behaviours.

Chapter 3. HexaShrink, from mesh representation to coding

Our rst experimentation is to encode the eight datasets presented just above with ZT in lossless mode, and to compare its CRs with those obtained with LZMA in the same context. Table 3.3 shows the CRs obtained with four combinations of transforms and encoders: HS P +LZMA: the transform for continuous properties proposed in HexaShrink (presented in Subsection 3.2.2, page 47) followed by LZMA; HS P +ZT: the same transform followed by ZT in lossless mode;

Haar+ZT: the transform HS P is substituted by the classical 3D Haar transform; CDF 9/7+ZT: the transform HS P is substituted by the popular CDF 9/7 transform.

Regarding the properties of the nearshore 0 environment, the best CR is obtained when LZMA is applied directly to the data, without any transform (line none). One level of HexaShrink decomposition decreases the LZMA eciency from 2.79 to 2.12, and to 2.11 for several levels of decomposition. By substituting LZMA with ZT, we expected an improvement, but nally obtain a poorer result, at 1.93. This fact reveals that the transform proposed in HexaShrink for the continuous properties is not optimal in term of compression performance. The reason is that this transform considerably increases the dynamic range of decomposed data, in order to ensure a fully reversible decomposition. Haar and CDF 9/7 better perform, but remains lower than LZMA after all. Our explanation bases on the environment physiognomy: smoother and more homogeneous.

Therefore the performance of an evolved encoder using the residual dependencies within multi-scale decomposition are less remarkable. However, regarding the three other environments, the classical wavelets Haar and CDF 9/7 combined with ZT outperform HS+ZT again, itself outperforming LZMA alone. Similar results are obtained for the permeability.

We can thus conclude that, globally, ZT seems to be a better solution than generic lossless encoders. Also, the transform proposed in HexaShrink seems to be less appropriate than Haar or CDF 9/7 to get the best compression performance on the continuous properties. On the other hand, it ensures a fully reversible decomposition, contrary to Haar or CDF 9/7. Indeed ZT needs integer values as input, a quantization before using it is thus required. Nevertheless, we will see hereinafter that this slight "degradation" does not necessarily have an impact on the visual quality and on the simulation results, which allows us to consider a ZT-based workow as "quasi-conservative".

Our second experimentation is to use the ZT coder in progressive mode, in order to evaluate its eciency in the context of compression at renable precision. By progressively decoding the bit planes, from the MSB to the LSB, a SNR curve can be drawn in function of the quantity of bits decoded, and thus for dierent precisions. Figure 3.22 shows the resulting SNR curves when the Haar transform is used. Each graph analyzes one of our eight continuous properties presented before. Each marker corresponds to a specic number of bit planes decoded. The rightmost marker of each curve corresponds to the CR obtained in lossless mode (Haar+ZT in Table 3.3), i.e., when all the bit planes are used. For purposes of comparison, the bit budgets corresponding to the CR obtained for the methods LZMA, HS P +LZMA and HS P +ZT in lossless mode (previously presented in Table 3.3) are also represented in each graph (by three vertical straight lines). As expected, a (de)compression workow at renable precision allows reaching bit budgets much more smaller than the lossless methods. In the same time, the SNR inevitably decreases according to the bit budget, indicating data less and less precise, and thus less and less faithful to the initial data.

Compression workflow at refinable precision for GMs

A trained eye would notice the shape of our SNR curves diers from the SNR curves classically obtained for images and videos (whose pixel value is commonly encode on less bits). Usually, the curves increase sharply to gradually reach a plateau at a certain bit budget. At this point progressive decompression can stop, because a higher binary budget would not signicantly improve the objective quality. Considering our results obtained on scientic data, the interpretation is more dicult because curves are incredibly linear. The same trend is noticeable in other scientic studies working on similar topics [START_REF] Cappello | Use cases of lossy compression for oating-point data in scientic data sets[END_REF]. This problem could be addressed by testing other metrics.

To have an idea of the impact of a (de)compression at renable precision in a context of geological data visualization, Figures 3. 23 and3.24 show the renderings of two given permeability properties decompressed at dierent bit budgets. For comparison purposes, the original data are also shown on the top corner left of each gure. On the rst gure we can oberve that the visual degradations are quasi imperceptible until a bit budget of 19.71 bits/value, which corresponds to a much lower bit budget than those obtained with the lossless methods. A similar observation can be done on the second gure until a bit budget of 18.85 bits/value. Moreover, as a teaser of the compandor eect on continuous properties (see next chapter), these gures also show the high quality of the renderings obtained for these two properties with only 12.65 and 12.17 bits/value respectively (bottom corner right of each gure), when a compandor is applied before the compression workow. The visual results are even more satisfactory.

All these results are very promising, and validate the fact that processing the continuous properties of geological meshes with a compression workow at renable precision does not necessarily aect the data severely for certain applications. This is the case for visualization, but we will see in the next chapter that the same observations can be made in the context of simulation.

Conclusion of this chapter

In this chapter we have shown that a suitable multiscale representation such as HexaShrink can improve the compression workow for the storage and the visualization at dierent levels of detail of geological meshes. HexaShrink is the result of a literature review of existing approaches dealing with VMs, in particular geological, and other more general scientic data. This solution gathers dierent wavelet-based transforms adapted to the heterogeneous components of the GMs. Firstly developed for visualization [START_REF] Peyrot | HexaShrink, an exact scalable framework for hexahedral meshes with attributes and discontinuities: multiresolution rendering and storage of geoscience models[END_REF], it decomposes an initial mesh into embedded hierarchical structure. The process is fully reversible and generates lower resolutions while preserving the mesh features across resolutions (as faults and properties coherency).

My rst contribution consisted in validating the skills of HexaShrink in a comparative study with upgridding techniques available in popular geomodellers, on a representative benchmark composed of eight meshes. Then, we used HexaShrink in the context of compression, by combining it with generic encoders. We observed that globally HexaShrink improves the coding performance, Chapter 3. HexaShrink, from mesh representation to coding but the eciency depends on the encoder and on the meshes features. As the mesh is heterogeneous, we also detailed the compression performance according to the mesh components. We noted a general improvement compression by using HexaShrink whatever the component, except for the continuous properties. For this data, HexaShrink combined to an encoder slightly increases the size of the encoded data, compared to using an encoder alone. The reason is that the continuous properties are oating-point data, with a high precision and distributed on a wide range, which makes their compression more challenging.

The second contribution was to experiment an evolved encoder as the ZT coder to better exploit the multiscale representation provided by HexaShrink. One conclusion is that the HexaShrink decomposition developed for continuous properties is not the most suitable for compression, although fully reversible and exact. Therefore we replaced it by classical 3D wavelets, and associated them with success to the ZT coder. We showed that compression at renable precision is possible, by gradually decompressing the data, which permits to reach bit budgets signicantly lower than with generic lossless encoders. Finally, our experimental results showed that, even at limited precision, the visual degradation of the mesh renderings can be negligible.

In next chapter, we continue to demonstrate the benets of our data representation, showing in particular that a renable precision has a limited impact on a simulation workow.

CHAPTER 4

Impact on simulation performances "It is not about the loss of information, it is about doing more science!"

John Dennis, 2013

Lossy compression has progressively integrated simulation. For few years now, it is used for visualization of large simulation outputs, and only recently rst tests combine it with calculation to relieve the workow. However, whatever the step compressed in the workow, it is primordial to evaluate its impact on simulation outcomes. While there are many metrics for visual assessment none metric objectively assess the impact on simulation since expectations dier according to the eld. To optimize the compression in simulation, we thus need to understand the requirements of professionals in geoscience. The discipline uses the simulation results to validate any change in simulation parameters. For this study, we rene the precision or resolution of our mesh properties thanks to HexaShrink. The quality of the rened data is evaluated upstream with objective metrics, and veried downstream after simulation with its results and a subjective validation made by engineers.

74

Chapter 4. Impact on simulation performances

Compression in simulation

In this chapter, we study the impact of renable precision on simulation, both for spatial and numerical precision. We especially focus on the impact of numerical precision induced by lossy compression in our literature review than on data reduction with resolution change. Although the latter is a standard but crude method in science to accelerate the computation process (by subsampling, such as upscaling/upgridding), it has recently been shown that data reduction has a impact on analytical results worse than compression, at equivalent binary quantity. This observation motivates our interest for lossy compression for simulation. The following review, mainly related to studies in the elds of uid mechanics and climatology, conrms our intention to quantify spatial and numerical precision of simulation data.

Until few years ago, scientists were suspicious about the use of lossy compression during a computational process, fearing the loss of signicant data. Usually accepted in the context of visualization, it was contested for simulation and analysis. But this is now changing. Slowly, in conjunction with the exponential increase in digital consumption, researchers are becoming aware of the growing diculties with bottlenecks and the increasing costs of computing resources reported in Section 1.2. Baker et al. ( 2014) explain that daily-generated simulation results are signicantly under-exploited considering their cost of production, due in part to limited storage capacity.

Lossless compression rst appears as a suitable solution, by reducing the data quantity while preserving the exact accuracy of processed data. Unfortunately, this conservative method provides poor performance on scientic data. The compression ratio indeed barely exceeds two [START_REF] Son | Data compression for the exascale computing era survey[END_REF] because of the high precision and the dynamic range of numerical data. Besides, it is known in multimedia that lossy compression allows improving signicantly compression ratios, with mostly undisturbing perceptual changes. Accordingly, our objectives in this chapter are to address the eects of compression on simulation workows and to understand to what extent its impact may be acceptable in simulation. Finally, we question the need to maintain the high precision of simulation data.

To illustrate the heavy process of simulation, we describe the Computational Fluid Dynamics (CFD) approach, used in our case study (cf. Subsection 4.2.1). Similar to geosciences, other elds using simulation are strongly impacted by processing huge data quantities, especially in climatology. Here, compression is increasingly looked at, as evidenced by the sharing of large datasets in open source to study the impact of lossy compression. However, the majority of studies still use compression more to improve visualization and storage than to integrate it into their workow. Our review focuses on integration work and classies the studies according to their elds, databases and codes. Besides, part of the review refers to quality measures commonly used to study the impact of compression. Subsection 4.1.2 is the core of our review: compression and evaluation methods are described in an applicative context. We note a growing use of compression and the establishment of routines in the choice of evaluation tools and methods. These routines are mostly borrowed from the research teams that develop the most popular compression tools for scientic data, namely SZ and ZFP. Referencing their studies can help tracking additional work and identifying innovative thinking.

Compression in simulation

Application for uid dynamic

Simulation is a numerical tool used by researchers and engineers to study and generate multiple cases, in order to better understand a complex system (gas injection, combustion, turbulent uid behavior etc.). It may provide a detailed vision of an inaccessible phenomenon, at dierent levels of resolution.

Its use extends to various elds, ranging from engineering design (vehicle, aerospace, etc.) to the study of the formation of the universe, using molecular dynamics (MD). At a smaller scale, MD can also simulate the behavior of materials and chemical reactions. But as mentioned in the introduction, our scope of interest is more focused on geosciences and in particular on applications using CFD for simulation. Biology domain also use simulation resources to learn about gene expression by simulating probabilistic and statistical events.

Fluid dynamics example

The study of uid motion by CFD in a limited volume involves solving systems of partial dierential equations (PDEs) using numerical methods as nite dierence, nite volumes, nite elements, or spectral methods. The scientic domains using CFD dier by the nature of uid and the context, such as climatology considering an air volume (atmosphere) for weather forecasting or climate change studies. Projects in this domain are carried out by leading international laboratories like NCAR, because the atmosphere is a complex area in interaction with adjacent layers, dened in the introduction of Chapter 2.

In a straightforward case of uid dynamics, as laminar ow into cylinder pipe, the issue could be analytically solved. It means the exact solution can be computed for any given point in the control volume, with pressure and velocity values as results. But considering complex problems, an analytical solution does not necessarily exist. In this case, we try to approximate the ows by a physical law formalized by parameterizable mathematical equations (Navier-Stokes equation for the study of turbulent ows, Darcy equation for the study of ows in porous media). The most complex situations may require additional components. For example, the modeling of the transport of chemical elements in the atmosphere (such as pollution) additionally uses an advection-diusion equation, which exponentially increases the numerical and temporal cost of resolution. For reasons of economy, the modeling of a system of equations is the result of a compromise between representation accuracy and computational resources.

After modeling the system of equations, the solution is calculated on the control volume. To obtain results at dierent points, the volume is discretized into small blocks (cells introduced in Subsection 2.1.2). The physical parameters of uid (velocity, pressure, etc.) are calculated for each cell. Their values are conditioned by the values of nearby cells interacting with each other.

The solution would certainly be more accurate by using an unstructured mesh but the structured case remains an attractive option because the mesh generation, as well as simulation code, are simpler. Finally, CFD obtains a discretized approximation of the real ow relevant to scientic exploration/analysis.

Depending on the size of the mesh and the precision required, the size of the output of large simulations could reach several terabytes or even petabytes, knowing that the time dimension multiplies the volume of data by the number of time steps. Furthermore, the number of outputs (velocity, Chapter 4. Impact on simulation performances pressure, etc.) tends to increase with the complexity of the system of equations discussed earlier in this section.

Selected domains

As mentioned above, dierent scientic elds use CFD simulation. The uid mechanic is the parent discipline. It includes several applicative domains (such as climatology, branches of geology, etc.) studying uid ows in various volumes of interest (from the atmosphere to reservoir formations). To solve the data processing problem and optimize the use of computing resources, several elds now embrace compression.

Tables 4.1 and 4.2 summarize information pertaining to datasets founding in the literature: domains, format, size and code function. The names of these data are written in bold. In the eld of uid mechanics, a rst comment is to note that the majority of them have sought to integrate compression into computational codes (such as LULESH, Nek5000 etc.) instead of applying compression to databases (Johns Hopkins Turbulence). Fluid mechanics being the mother discipline, we can foresee that the integration of compression into simulation is an ultimate goal, from which the ndings will later be applied to other areas. As a second comment, looking at the example of CMIP in Table 4.1, we observe an increasing evolution of volumes generated within the same project over time. From 2000 to today, the few terabytes of scientic data have become several petabytes. Nevertheless, the cost of simulation and the energy required to transfer/store large volumes of data are underestimated and contrary to claims. Data are indeed likely to be exchanged and copied many times since they are freely accessible to maintain transparency on their analysis. CESM-LE produced in 2005 2.5 petabytes of raw data, post-proceeded to obtain 170 terabytes, which compose CMIP Phase5 [START_REF] Meehl | Overview of the coupled model intercomparison project[END_REF] in netCDF format. While CMIP Phase6 [START_REF] Eyring | Overview of the coupled model intercomparison project phase 6 (CMIP6) experimental design and organization[END_REF] is estimated to be between 20 and 40 petabytes. To give an idea and visualize what its storage implies, consider a standard external hard drive whose capacities are between 2 and Other meteorological objects, such as hurricanes, are tested to assess compression impact (ISAB-EL). Hence datasets in climatology and meteorology generally consist in 2D, 3D or 4D grids modelling the evolution of the global or local atmosphere. Temperature, pressure, humidity parameters and velocity vectors compose the information stored in grid cells.

Climatology -oceanography

Compression in simulation

With much fewer public attention than climatology, other domains generate equivalent or even greater data volumes. This review aims at broadening the scope by considering additional domains and projects in cosmology (HACC/NYX), quantum computers for chemistry.

Geology The domain is rich and comprises several sub-disciplines using simulation. In reservoir engineering, our case study, ow simulation studies the motion of phases (gas, oil and water) in porous media as detailed in Subsection 4.2.1. Furthermore, uids can also be the geological materials, as it is the case for convective process in Earth mantle/core or magneto-convection in the outer core [START_REF] Schmalzl | Using standard image compression algorithms to store data from computational uid dynamics[END_REF]. If the task is to study the rupture of geological material, other simulation approaches than CFD could be used and focus on small objects, such as drill cores [START_REF] Bouard | Generation of DEM packings from RingMesh geological models[END_REF] or on larger objects with earthquakes. While, for the moment, little work is still devoted to the use of lossy compression for simulation, such a practice is more frequent for the transmission and storage of seismic data.

An important part of the publications listed here results from collaborations with simulation projects. Some datasets may be condential and unreferenced. For the sake of repeatability, we listed only works whose data origins are known and referenced. The ones marked by an asterisk (*), are referenced by Scientic Data Reduction Benchmarks By identifying the specic needs in rate, ratio, and accuracy for each application, they promote the notion of a specialized compressor to meet particular requirements.

To accelerate the execution of quantum simulation in chemistry, Gok et al. to simulate a turbulent pipe, or the eect of wind on a wing. They apply a compression based on the Chebyshev transform to the outputs. This preliminary study visually evaluates the impact and invariance of results through an statistical analysis. Otero et al. discusses the fact that 70 % of the data could be truncated, while preserving accuracy for scientic needs. In the article conclusion, compression for checkpointing/restart gures as future works. Indeed, this step represents an accessible bottleneck in the workow, easily optimized: by reducing the quantity of processed data, the simulation becomes faster. The ultimate purpose would involve integrating the method into the simulation code to optimize, for instance, the memory allocation. But for the moment, the majority of studies are examining the feasibility of lossy compression by visualizing and quantifying its impact through objective metrics.

Objective evaluation

We now focus on the methods used to evaluate the impact of lossy compression in simulation. As with image and video compression, the impact on simulation can be assessed using conventional objective measures. However, these measures may not be relevant since they are not adapted to the scientic application.

To optimize scientic data visualization, the TTHRESH project developed by Ballester-Ripoll et al. ( 2019) evaluates the impact of compression on volume data from the Johns Hopkins turbulence database. The choice of classic objective measures (such as relative error, RMSE, and PSNR) seems relevant if it is limited to a visual assessment. The same argument is questionable when considering simulation data. Therefore it would be wise to assess the eects of data compression in simulations using physically motivated metrics", as stated by [START_REF] Laney | Assessing the eects of data compression in simulations using physically motivated metrics[END_REF]. Their study uses, among others, fpzip compression, and integrates it into three simulation codes. Among them, there are a Lagrangian shock-hydrodynamics code from LULESH and an Eulerian higher-order hydrodynamics turbulence modeling from Miranda. Metrics motivated by physics are symmetry shock radius for the rst case, and Rayleigh-Taylor instability and spectrum of perturbations in the second case. This work is a landmark study, being the rst applying lossy compression to the Chapter 4. Impact on simulation performances physics state of simulations as a strategy for mitigating the data movement bottleneck expected on future systems". He considers that compression ratios under 3:1 do not aect the simulation results.

The usage of greater compression ratios may be dependent on applications.

Nevertheless, it is undeniable that objective metrics are easily implemented and fast to compute (because roughly computed from point to point dierences as explained in Subsection 1. Such a crossing between dierent evaluation methods is classical in multimedia. By combining classical objective metrics and subjective evaluations (like MOS dened in Subsection 1.3.3), we could obtain reasonable objective thresholds for which a subjective alteration is acceptable for human perception. For scientic data, a subjective evaluation depends on the competent opinion of specialists who interpret the quality of (de)compressed data.

Subjective evaluation

To subjectively assess the performance of several compression tools, Baker et al. (2019) collaborate with climatologists on the project CESM-LE. Four distinct variables are compressed (2D surface temperature, 2D clear sky net solar ux, 3D grid box averaged cloud liquid number, 2D convective precipitation rate) at renable precisions. Thereafter, they are subjectively evaluated by climatologists, who are asked whether the data are the same, or dierent from the original. To identify the relevant objective metrics, i.e., correlated with expert opinion, ten objective metrics are tested, including MSE, PSNR, SSIM (as dened in Subsection 1.3.3). For this study, SSIM obtains the closest results to the subjective opinion of specialists. The SSIM threshold required to preserve the data precision in this study is higher than the SSIM required for medical data to ensure diagnosis [START_REF] Baker | Toward a multi-method approach: Lossy data compression for climate simulation data[END_REF].

In conclusion to this state of the art, we demonstrated that compression in simulation is a relatively recent concern. Although the publications on this topic are still in limited number, they are carried by renowned laboratories, and the number of related papers rises. Some works use classical methods (tools, quality assessment) borrowed from multimedia data compression, while Similarly with the rst one, we focus on the application of simulation and consider metrics specic to the eld of reservoir engineering. We also specify subjective criteria, as in the

Evaluation of multiresolution progressive compression method

83 second publication, that we crossed with objective methods to validate our compression workow on an extensive simulation benchmark.

Evaluation of multiresolution progressive compression method

To evaluate the impact of a lossy compression stage in a simulation workow, we rst design a case study (introduced in Subsection 4.2.1). It includes an input mesh, called lundi, which is none other than the mesh#6 in the benchmark introduced in Section 3.3. This mesh has been already processed with HexaShrink and generic encoders for lossless compression. This enables us to make the link with a complete simulation workow.

Subsections 4.2.2 and 4.2.3 present simulation results using our scalable progressive representation. We perform tests according to renable resolutions and renable precisions. Results are objectively and subjectively assessed to validate the coherency of the dierent metrics. Finally, the simulation performance obtained from our compression workow at renable precision is compared with those obtained from the two most prominent compression tools: SZ and ZFP.

Proposed uid simulation workow

Our simulation workow represents a typical case of reservoir engineering, as detailed in the three next subsections. First, we describe the input mesh. Second, we shortly present the simulation of liquid phase ow through the reservoir. Third, we detail the simulation results, which are analyzed to optimize the exploitation of the reservoir. Since simulation runs are generally repeated a large number of times for sensitivity studies, the aim is to accelerate their calculation by decreasing the number of cells. The method is validated if the outputs of the simulation are not modied or if the change remains acceptable. To evaluate our representation of the mesh at renable precision and resolution, this argument is decisive, as explained in Subsection 4.2.1.

Versatile RM

Our RM is called lundi 3 . If the global structure of the model has been briey described for the compression benchmark (cf. Section 3.3), it is further detailed in this section. We designed a global subsurface morphology integrating three continuous vertical faults. They are not aligned along grid axes, and possess dierent osets, to emulate mildly complex environments. This morphology may also challenge compression algorithms. This model can be discretized with dierent grid dimensions.

For our evaluation, we choose a grid with 128 × 128 × 32 cells to allow reasonable simulation times, with respect to our extensive comparative workow. The average cell size in this lundi realization is about 1.7 × 1.7 × 0.95 m 3 . This shape ratio for the tiles is common in (sedimentary) geology for modeling mesh cells. It is consistent with progressive horizontal ne deposits of material over the years.

lundi topography arises from a realistic geological context. It models a side of an anticline structure (classical in hydrocarbon trap reservoir study), as spotlighted by Figure 4.1. Therefore the most elevated corner of lundi corresponds to the top of the anticline (at 3360 m depth). The 3 lundi refers to the Atlantic pun in the Icelandic language (macareux moine in French), a pelagic and relatively pacic bird in the rather competitive world of Petrel and Skua seabirds, or associated software. The red dot indicates the water breakthrough.

opposite corner on the diagonal is situated 50 m below. Our benchmark is thus based on a quarter of the anticline modeled by lundi. The global shape of lundi as well as its structural discontinuities, have been thought as a challenge for multiscale visualization and faithful upgridding for simulation.

Two-phase ow simulation

We use for our simulation the black oil model [START_REF] Thomas | Reservoir simulation of variable bubble-point problems[END_REF], and limit the study to two liquid phases: oil and water (we ignore gas phase for the sake of simplicity). The system can be described by physical laws (conservation of mass and Darcy law). It can be modeled by a mathematical system composed of non-linear partial dierential equations. Finally, unknown system parameters, such as pressure and saturation, are computed for each cell, for each time step.

Two-phase ow simulation is thereafter performed on lundi. It is meant to predict oil production in a two-phase reservoir, driven by water injection. At the initial time, the two phases in the reservoir are horizontally stratied, with oil above water. We simulate a quarter ve-spot conguration inspired by (Lie, 2016, Chapter 5.4.1). It involves two wells, as shown by Figure 4.1: one producer P 1 and one injector I 1 (among the four ones {I i } localized on the anticline base). Water is injected by I 1 in the lower part of reservoir. The water pressure pushes the oil through the reservoir up to the producer P 1 (distant from 300 m). Injector pressure and producer rate remain constant, respectively set at 300 bars and 300 m 3 per day.

Reservoir production overtime

Oil production is evaluated according to various reservoir properties and injection settings. One of the reference indicators is the estimated water cut (WC), or the ratio of water produced in a well compared to the total liquid volume. It is recorded over a period of time, for a certain time step, and can be illustrated by the typical WC shape in Figure 4.1-right. Furthermore, output analysis can also focus on other parameters known to reservoir engineers, such as well pressure or oil production, to be sensitive to the slightest change in the input data.

WC is a valuable source of information for decision anticipation for eld exploitation. The inection point marked by a red point on the curve Figure 4.1 is the water breakthrough. It denotes the rst water arrival to the production well. After this date, extracted liquid will progressively contain more and more water. Hence, to avoid eliminating water by expensive post-treatment or surface equipment, this is interesting to delay this date by changing the exploitation method.

Simulation is a powerful tool to predict the water breakthrough and to determine the best option, by testing dierent methods and conditions. Hence simulation can be run a large number of times with varying parameters (reservoir sensitivity analysis). It therefore seems reasonable to limit its calculation cost by using a lower resolution grid.

Well production metric

To check whether standard upscaling is suitable for simulation, reservoir engineers carefully inspect the simulation outcomes and compare them to outcomes obtained with the reference RM. This approach is common to validate the accuracy of a RM at a lower resolution. We extent it to the assessment of compression at renable precision through a complete comprehensive compressionsimulation workow, illustrated in Figure 4 provided by simulation of the original RM. We denote by WC the water cut curve provided by simulation of compressed RMs. We separately address the impact of a global renable resolution of mesh (for upscaling/upgridding), and that of renable precision of properties (here permeability).

The dierences between WC and WC curves should be evaluated with respect to the expected Objective evaluation of renable representations With HexaShrink, we generate dierent resolutions of the RM, and then several versions of the RM at full resolution lled by continuous properties at renable precision. During this work, we mainly focus on permeability for renable precision.

When considering an upscaling on a RM, several quality indicators, related to petrophysical properties, have been proposed to evaluate the resulting information losses [START_REF] Preux | About the use of quality indicators to reduce information loss when performing upscaling[END_REF]. When considering a precision renement, the conventional objective metrics are classically based on dier- The two-phase ow simulation is applied on the dierent lundi resolutions, resulting in distinct WC curves. They are compared with the reference curve obtained from the original RM. The results in the nearshore 0 environment are displayed on Figure 4.5. The results for the three other environments are available in appendices (page 124). First, we observe that the WC curve simulated on the reconstructed resolution res-0 is identical to the reference WC whatever the environment.

Evaluation of multiresolution progressive compression method

This conrms that the HexaShrink decomposition is fully reversible for simulation. Then, we observe that the shape of the reference WC curve is modied depending on the lundi environment (nearshore 0 , nearshore 1 , nearshore a , fluvial). The water breakthrough in fluvial environment happens 30 days after the water breakthrough observed with nearshore 0 environment. Thus, the petrophysical properties signicantly aect the simulation outcomes.

Considering the resolution res-1 in the nearshore 0 environment, with a number of cells reduced by eight, its WC curve stays close to the WC curve reference. Then, according to our subjective metric (explained in Subsection 4.2.1) its WC curve is acceptable. Compared to the equivalent WC curves obtained in three other environments, the outcome obtained with nearshore 0 is the best. This dierence may be explained because of the high permeability of nearshore 0 environment (mostly green according to the permeability color scale). In that case, the modication of permeability would have less impact on the simulation. This allows us to focus mainly on the impact of the change in structure due to the use of lower resolutions. We infer that the structure of the mesh res-1 remains suitable to preserve the quality of simulations. Thus, HexaShrink is a promising upgridding technique, while further improvements are necessary for property upscaling.

At equal quality, the HexaShrink-based simple upscaling is computationnaly ecient. While the original mesh required more than three hours for simulation, only ten minutes are needed with the resolution res-1 on the IFP Energies nouvelles supercomputer Ener440. For comparison purpose, the same simulations run on a scientic laptop (with Intel Core i7-6820HQ CPU @ 2.70 GHz processor and 16 GB RAM) would take much longer: seventy-two hours on the full resolution, and four hours on the resolution res-1.

The results presented in this section highlight the power of upscaling/upgridding in saving simulation time but also the sensitivity of the simulation to the features of the input mesh. A mesh scaling is a complex process that requires care in nding the optimal approach for each component of the RM. From the literature review, we know that the best methods are dynamic approaches presented in Subsection 3.1.1. Unlike uniform merging of cells (eight by eight), the dynamic approach only uses cells outside the ow path to limit the impact on ow dynamic. As ecient upscaling methods use Local Grid Renement (LGR), we used it to correct the position of the wells.

As specied when describing the simulation in Subsection 4.2.1, the vertical wells cross through the cell columns located at opposite corners of the mesh. Considering the original data, {I i } is located at surface coordinates (1,1). Using lower resolutions, {I i } is still located at the (1,1) but is spatially translated because of cell dimensions increase. LGR method has been applied around the wells to avoid this issue. It retains the full original resolution on cell columns crossed by wells. indicates the data with the best metric evaluation.

Next, we look at objective metrics whose assessment is correlated to visual expectations and especially to simulation validation, supporting an enlightened use of compression for simulation.

Successively considering the objective metrics MRE, SNR, nRMSE, Λ-SNR, we note several anomalies. As an example, the SNR is equal to 78.78 dB with Λ 1 , and it is superior to the SNR obtained with Λ 0 (58.23 dB). This would mean that the Λ 1 is more appropriate than Λ 0 . Considering our previous conclusions, it nally appears that the standard SNR assessment is contrary to visual appreciation. We note the same anomaly using nRMSE. However, Λ-SNR and MRE indicate the Λ 0 is of better quality, as visually approved. Considering now the simulation results, Λ 0 yields identical mention, while the simulation outcome is only considered as correct with Λ 1 . Subjective evaluation of simulation results is in this case correlated to visual intuitions. Because simulation is a decisive aspect, we observe that Λ 0 nally generates higher quality data, in terms of visualization but also of simulation. Data quality is foreseeable regarding Λ-SNR and MRE, while SNR and nRMSE provide misleading assessments.

Simulation validates visual quality and benet of the use of the Λ 0 compandor.

Λ-SNR and MRE corroborate visual appraisal and simulation results.

This example questions the blind use of some metrics (SNR, nRMSE) on scientic data (with their particular distributions), whose assessment could be counter-intuitive.

Evaluation of multiresolution progressive compression method

93 Whatever the objective metric or the compression method, we observe that as expected the global quality decreases with the bit budget (MRE and nRMSE increase, SNR and Λ-SNR decrease). Depending on the relative position of the curves, we can rank the eciency of the methods.

Considering MRE and nRMSE (graphs on the right), the lower curve identies the best method. For SNR evaluations (graphs on the left), this is the highest curve. As previously noted, the evaluations provided by SNR and Λ-SNR are conicting despite their mathematical relation. In the same logic, we consider that the Λ-SNR assessment is the most coherent because it is validated by subjective results, while SNR is not consistent with what simulation would suggest. Indeed, regarding subjective evaluation provided by marks color of the SNR curves (bottom-left graph), the simulation results based on Λ 1 are acceptable for a bit budget around ten bits. At equivalent bit budget, Λ 0 yield identical simulation results. We thus expected that the SNR curve for Λ 0 would be higher, synonymous of better quality. However, it is the reverse. This conrms that the SNR leads to faulty reasoning on these data for reservoir simulation.

Finally, by focusing on consistent objective metrics (top graphs), the higher the Λ-SNR or the lower the MRE, the greater the chances of obtaining suitable simulation results. This seems to provide an objective limit (underlined by the blue dotted lines), below which the data are no longer suitable for simulating with enough accuracy. For a better understanding of the graphics, triangles are attached to the y-axis. They represent a balance between subjective and objective evaluations.

Their red tip tends toward poor data quality, while their blue base indicates higher quality.

In the remainder of our report, we only focus on Λ-SNR graphs considered as more readable. Therefore Figure 4.8 keeps the upper left graphic from Figure 4.7 and adds on the right the Λ-SNR graph for fluvial permeability (the results for all the properties are available in appendices). In the same logic as before, a graduated evaluation scale (gured by the color triangle) can be drawn, linking objective evaluation to the subjective one. This allows us to concur that the identical threshold varies according to the property and its environment. Considering the fluvial environment, identical simulations are obtained for Λ-SNR quality evaluated above 55 dB, while 45 dB is sucient for the nearshore 1 environment. The dierence between the two thresholds can be explained by their level of complexity. The fluvial environment would require a higher objective quality level to preserve channel objects.

In all the studied cases, the use of the Λ 0 compandor enhances the renable representations of the permeability. Considering the Figure 4.9, we can better perceive this compandor eect. First of all, whatever the value of the alpha parameter used for Λ , the data processed by Λ are distributed between 0 and 2 nbits-1 . If the shape of the distribution remains unchanged using α = 1 (linear), histograms are better spread using α = 0 (logarithm). Indeed the height of the rst bar decreases for fluvial (bimodal: channels and oodplain) and tends to disappear for nearshore 1 . The small values, originally concentrated between 0 and 250 mD, are dispersed to several distinct bars. The use of Λ 0 gives them more weight, and preserves them during compression. Concretely, from a physical point of view, it is necessary to preserve the lower permeability values, to save impermeable barriers and to maintain ow. Results for all the other environments are available in appendices (page 130).

The identical threshold lies between 40 dB and 55 dB. It means that whatever the environment we certify that, above a Λ-SNR equal to 55 dB, simulation results are identical, and those for a considerably reduced bit budget. Initial distribution is between [0, 20 000], then Λ α,25bits using α = [0, 1] distributes the properties between [0, (2 25bits -1)] (bottom).

Chapter 4. Impact on simulation performances

Λ-SNR

Comparative study with state-of-the-art

To challenge our method and check the positive eect of Λ 0 , we add to our graphs the curves of the probably two most popular compression tools for scientic data : SZ and ZFP (introduced in Section on their compression performance. In addition, the point wise relative error option of SZ was also been tested. It integrates a logarithm transform mapping before prediction, and is akin to our Λ 0 approach.

Taking a step back, the results of the other tools conrm the previous observations, notably the identical threshold (dotted blue line) at 45 dB for nearshore 1 environment, and at 55 dB for fluvial environment, below which the data precision is no longer sucient for simulation. We observe that Λ 0 , as the logarithm mapping, also improves the compression performance of both compression tools. Focusing on the fluvial permeability (right graph), the data processed by Λ 0 (solid lines) obtain identical simulation results for Λ-SNR above 60 dB. At similar objective/subjective quality, the regular methods (dotted lines) require about 5 bits more per value than the same method using a Λ 0 , demonstrating the interest of a compandor stage. We observe that Λ 0 , as the logarithm mapping, also improves the compression performance of both compression tools. Focusing on the fluvial permeability (right graph), the data processed by Λ 0 (solid lines) obtain identical simulation results for Λ-SNR above 60 dB. At similar objective/subjective quality, the regular methods (dotted lines) require about 5 bits more per value than the same method using a Λ 0 , demonstrating the interest of a logarithm compandor stage.

Λ-SNR

More specically, we observe that our method gives results comparable to the other tools (SZ and ZFP). With Λ 1 our alternative is just below ZFP, the best method (without Λ 0 use). Lindstrom's tool is indeed known to be the most ecient on volume data. Like our alternative, it combines a transform with a ZT coder. Thus, it seems consistent to obtain comparable trends.

ZFP gets its superiority from the sub-blocking process explained in Subsection 3.1.3. The equalization of the exponents (oating-point writing) within each subblock absorbs high variations of data scales. This is also one of the benets of the logarithm function. Therefore combining Λ 0 with ZFP increases the compression performance but not as much as it does for other methods that take the lead. Considering SZ, with Λ 0 we obtain results equivalent to those obtained by SZ with the point wise relative error option.

Although not evaluated in term of execution speed for HPC implementation, our approach (regarding its compression performance) provides suitable results comparable to state of the art compression tools for simulation. Therefore, incorporating precision layer to HexaShrink consti- To conclude, continuous properties can be generated at renable numerical precision without visual noticeable degradation, while the simulation outcomes remain identical to the reference. To guide the user through the compression process, we identify objective metrics correlated to professional expectations. Our tests show performant compression provides identical simulation results using less than the half of the binary quantity originally used for permeability property (nearshore 1 environment): among 12 bits per value instead of 34 bits. Moreover, the use of compandor functions seems to provide consistent predictions in addressing the properties distribution (here permeability), by better taking into account the physical laws that inuence the continuous properties. This example shows that the knowledge of the data (in particular their specic distribution) can considerably improve their processing. Such observations were drawn whatever the geological environment chosen for the simulation of our model lundi.

Lastly, a comparative study with SZ and ZFP tends to corroborate the above observations.

There seems to be a minimum quantity of binary information to recover the essence of the RM.

This value depends on the data, since our subjective/objective thresholds vary according to the property environment. But even using dierent tools, there seems to be an agreement to an objective threshold (blue dotted lines in graphics of Figure 4.10). We aim at nding heuristics on the essential binary information to preserve for lossy compression, in the spirit of the entropy bound for lossless source coding. Delving into the extensions of the above contributions, taking advantage of the data anisotropy or investigating the bit plane encoding may oer novel insights.

Exploiting spatial properties distribution with variograms

A variogram, dened from the Equations 4.5 and 4.6, is a geostatistic tool used to describe the spatial variability/continuity of properties in specic directions. It is dened as the variance of the dierence between values at two points separated by a distance h, equivalent to a number of cells.

γ(x i , x j ) = 1 2 Var(Z(x j ) -Z(x i )) .

γ(h) = 1 2N (h) N (h) i=1 (Z(x i + h) -Z(x i )) 2 . (4.6)
The (uniform) distance separating points x i , x j is denoted by h, while N (h) denotes the number of pairs of points (separated by distance h).

Along the three main directions ı, , and k, three 1D variogram curves can be computed. As illustrated in Figure 4.12, the curve is described against its sill (asymptotic limit), and its range (distance in which the dierence of the variogram from the sill becomes negligible). Estimating the variogram could gure as a preprocessing to compression, to determine some preferential directions in the data. An anisotropic wavelet decomposition could then be applied to benet from this feature. This kind of approach has been already proposed in Christophe et al.

h γ(h)

(2008) for hyperspectral data. This work uses wavelet packets on approximation and detail subbands localized along directions pointed by the study of variograms.

Our geoscientic data also presents anisotropy. Here we focus on the petrophysical properties introduced in Subsection 3.5.3, whose the anisotropy orientation should be along the vertical axis k but is not necessarily intuitive, considering the various cell dimensions. The variogram of the porosity in the nearshore 1 environment is displayed in Figure 4.13 (the other environments can be found in appendices, page 128). The range distance initially expressed in meters (4.20 × 4.20 × 1

Chapter 4. Impact on simulation performances meters), is converted into number of cells, because of property discretization and average dimensions cells. We observe that the range along ı and  are ten times higher than the range obtained along k. This means that along  for instance, two cells separated by less than 55 other cells are likely to be correlated. Above that distance, a weaker link is expected along the k direction. Such a spatial distribution is common and these variograms are comparable to those obtained on the Tarbert formation (from SPE10 model), equal to (14,32,3).

Because the range computed in the k direction (depth) is short, we assume that the information changes fast. Classical dyadic decomposition can leave behind information on detail subbands that are still correlated. Preliminary results with anisotropic decomposition Figure 4.14 illustrates the two approaches compared in this section. The volume data on the left (preprocessed by Λ 0,25bits ) is decomposed with a 3D dyadic wavelet, generating a structure displayed in the middle of the gure. For this experimentation, we used 25 bits for Λ (instead of 34 bits) and two decomposition levels (instead of ve) for sake of simplicity on the gure. This method is compared to the wavelet packet decomposition illustrated on the right, more suitable to handle data anisotropy. The previous study of variograms recorded a very short range along k, showing sudden spatial changes and concentration of information in this specic direction. Despite the rst dyadic decomposition, the LLH 1 detail subband is still correlated. We can recognize in this subband an approximation of the original data at a lower resolution, rather than detail data. It should therefore look more like noisy data, and be composed of small values. At this place, the wavelet packet decomposition uses two additional decompositions to better decorrelate information on this subband.

To assess the interest of using a wavelet packet transform to take into account the data anisotropy, we compare the weighted entropy [START_REF] Shannon | A mathematical theory of communication[END_REF] of the two approaches. The entropy designates In case of decomposed data, the weighted entropy is generally used. It consists in weighting the entropy of each subband by the ratio between its number of cells and the total number of cells:

wH = N k=1 nb cells in subband k nb cells in data H(k).
Computed on the entire data preprocessed by Λ 0,25bits , the weighted entropy is close to 19 bits per value, whatever the environment (all the results are in appendices, page 128). Considering the property in the nearshore 1 environment, the classical 3D dyadic decomposition reduces the weighted entropy from 18.92 to 15.62 bits per value. With the wavelet packet decomposition, the weighted entropy reaches 14.99 bits per value, which is thus lower.

As demonstrated previously, our evaluation method introduced in Section 4.2 is complete and relevant to assess compression in a simulation workow. We therefore use it in the same way in this section to assess the anisotropic wavelet decomposition. For this experimentation, we made again ve decomposition levels and 34 bits, to be compared with former results. Results displayed on So far, the performance of our alternative method was close to the concurrence. The anisotropic approach, based on Λ 0 , a wavelet packet decomposition and an adapted ZT structure [START_REF] Christophe | Compression des Images Hyperspectrales et son Impact sur la Qualité des Données[END_REF]) surpasses the concurrence results. This experimentation conrms the positive eect of Λ 0 compared to Λ 1 : solid lines are always above dashed lines. In addition, the objective thresholds determined earlier, correlated to subjective evaluation, are still valid (dotted blue lines). In the fluvial environment, the best result was held by SZ (point wise relative option). The lower data precision able to satisfy the simulation workow is generated with a budget of 8.92 bits, and assessed by a Λ-SNR equal to 53.77dB. For a lower bit budget of 8.83 bits, the new anisotropic method (black solid line) generates a data at higher objective quality, equal to 56.77 while maintaining identical simulation results.

In short, wavelet packet decompositions show promising performance, and would require in-depth experimentations.

Thresholding at necessary precision

The objective of this last experimentation is to determine the precision required for suitable simulation data. Up to now, our results designate thresholds by using adapted objective metrics (blue dotted lines). Nevertheless, the threshold value varies according to the data type. Therefore our research pursues a quest to address the minimum information, regardless of the handled data. able precision decreasing the bit budget (a mark per precision) using our method (wavelet packet or dyadic decomposition), SZ and ZFP (using or not compandor), objectively evaluated by Λ-SNR.

Λ-SNR

Subjective appreciation of simulation results is represented by a color code.

Chapter 4. Impact on simulation performances

We focus henceforth on encoding parameters, in particular on the ZT activity measured for each bit plane and below dened. The notion of bits per bit plane has been particularly studied by [START_REF] Wang | Compression ratio modeling and estimation across error bounds for lossy compression[END_REF] to describe the ZFP activity, while others have detailed the encoding process of SZ.

The standard ZT coder [START_REF] Shapiro | Embedded image coding using zerotrees of wavelet coecients[END_REF] uses four symbols to encode wavelet transformed data (POS, NEG, ZTR, IZ). They determine the status of a coecient for a particular bit plane considering its previous state in the higher bit plane. As explained in Subsection 3.5.2, the coecients are successively visited and compared to a decreasing threshold equal to 2 bit plane . If the value of the coecient is higher than the threshold, the coecient is considered signicant, and is included in the list of signicant coecients. From this bit plane, the coecient is encoded using POS and NEG symbols. On the contrary, if the coecient is lower than the threshold, it means that all its bits considered so far were zero. Depending on its location in the decomposed structure, a coecient gures as the ZT root (encoded with a ZTR symbol) or belongs to a ZT (by consequence it is implicitly represented by a root element) while IZ is for isolated zero.

Consequently, the state of a coecient changes continuously using lower thresholds. From unexpressed, it is potentially identied as a ZTR or IZ and nally becomes signicant. We dene the activity as the number of coecient state changes for a bit plane, with nSIG the number of signicant symbols (POS, NEG), and nZ the number of zero symbols (ZTR, IZ) at the i th bit plane. The activity is relative to the total number of coecients noted C , and dened by: ZT activity

(i) = nSIG(i) -nSIG(i + 1) + nZ(i) -nZ(i + 1) C (4.7)
Figure 4.16 shows the curve of ZT activity according to the bit plane for the permeability in the nearshore 1 environment, preprocessed by Λ 0,34bits and decomposed by the CDF 9/7 wavelet. This curve describes the progressive encoding process of binary data from the MSB to the LSB.

Consequently, this reverses the quality progression (Λ-SNR, bit budget) drawn so far, that showed the precision decreasing from left to right using less bit budget along the x-axis. Now, reading along the same axis involves using more and more binary planes, thus increasing the precision of the reconstructed data.

We interpret the ZT activity as a clue on the encoding evolution and allows identifying dierent phases. At beginning, using a partial number of MSB M planes, ZT are few but constitute large structures (comprising a large number of zero coecients). Consequently the starting activity is low because only few nZ and nSIG are encoded. Progressively, the number of ZTs tends to increase by splitting precedent large ZTs into smaller ones, while ever more coecients become signicant. For the example of Figure 4.16, the ZT activity takes o from the 36 th bit plane and peaks at the 31 th bit plane. Then, the ZT activity decreases to zero as well as the number of ZTs, since all coecients become signicant using an ever lower threshold. From the 22 th bit plane, we observe that the ZTs vanish. However, the ZT can continue to encode bit planes, even though they no longer have a structure. Our postulate is that, beyond this bit plane, information consists of incoherent data that does not contribute to the increasing of the precision for the reconstructed RM.

Besides, if considering the renderings, as soon as the activity returns to zero, i.e., once the 22 th bit plane is processed, the (de)compressed data is already very close to the original data visually (indexed by 8 in the Figure ). The decoding of the subsequent bit planes will only add very light visual details. If considering now the simulation accuracy, as soon as the activity returns to zero, our simulation results are already identical. This can be seen thanks to the color of the marks on the activity curves that allude to the subjective evaluation, as previously done in Subsection 4.2.3.

Finally, it would seem that we learn more about the data precision required to have a good visual quality and also an unbiased simulation outcomes by interpreting ZT activity than by analyzing objective metrics such as the curve Λ-SNR, displayed at the bottom of the gure.

Figure 4.17 shows the curve of ZT activity according to the bit plane for the permeability in the nearshore 1 environment, but this time preprocessed by Λ 1,34bits . The curve is quite similar than with Λ 0,34bits . A peak is observable between the 36 and 18 th bit plane, although less pronounced.

But the conclusion is the same: as soon as the activity returns to zero, the visual quality is good, and the simulation accuracy is already considered as identical. It conrms the fact the ZT activity seems to be a relevant information to assess the quantity of data required for simulation, whatever the compandor features.

To summarize, in this last section, we showed that the compression workow could be probably improved, rst by taking into account the anisotropy of some continuous properties. It could be simply done for instance by changing the dyadic decomposition by a wavelet packet decomposition, which enhances the decorrelation of some subbands neglected by an isotropic approach. Second we showed that the decrease of ZT activity during (de)compression seems to indicate the minimal data quantity required for obtaining accurate simulation outcomes even with a limited precision (as for displaying very nice reconstructions, with negligible visual losses). Finally, we also showed that, even if the bit planes subsequent to the peak of activity do not further improve the quality of the decompressed data, the curves of the objective measurements (Λ-SNR) continue to increase linearly, demonstrating an illusory" precision. Actually, the linearity of this curve does not allow us to distinguish the signicant data from the rest considered as noise by [START_REF] Natarajan | Filtering random noise via data compression[END_REF]. 4.17: Correlation between ZT activity (top-curve) and objective evaluation by Λ-SNR (bottom-curve) of permeability from nearshore 1 environment generated at renable precision by increasing bit planes number for reconstruction while our compression alternative using Λ 1 .

CHAPTER 5

Conclusions & perspectives

Conclusions

This work evaluated the relevance of compression for reservoir meshes, from visualization to simulation through lossless and lossy storage, along four chapters. The rst chapter provided context on the growing concern of data handling in data-intensive science. In a second chapter, we focused on geosciences and volume meshes used to model objects and simulate geological phenomena and reservoir behaviors. The third chapter reviewed methods employed to deal with huge datasets in various scientic elds. We then introduced a comprehensive multiscale approach for the dierent components of geological volume meshes. Its evaluation for visualization and compression shows the benets of combining embedded scales with both generic and renable data encoding techniques.

The impact of rening spatial resolution and numerical precision on a complete simulation workow is nally assessed in Chapter 4 on the lundi reservoir model. We thereafter provide more details on our contributions.

Compression in the domain of multimedia is very commonly used to reduce data quantity and improve the distribution of sound, images and videos. Compression standards (jpeg, mp3) have gradually invaded our daily digital lives, without end users being really aware of their benets and implications. Compressing digital data often involves the concatenation of regular components following a logical scheme. The singularity of a compression algorithm often lies in the nature of its components, chosen to meet specic needs: execution speed, compression ratio, quality. It is hardly optimal on all fronts. Hence, it is therefore essential to identify user needs, to nd a balance between theoretical advances and practical concessions. This allows to propose adaptable and perfectible models for various data sources. Contrary to multimedia, data formats in sciences is not well standardized, as they are not in-Chapter 5. Conclusions & perspectives tended, originally, to be exchanged amid a wide audience. However, their quantity exploded as well in the last decade, in a quest for better reality representation, leading to increased model precision and resolution. Technological advances on computational resources and facilities (high-performance computing, gigantic data centers) allows the processing of ever more detailed data. Yet, this system is gradually touching its limits, and bottlenecks are appearing in data storage, transmission, or computing. To address these growing issues, researchers have being re-examining compression for scientic data, in a variety of application elds with very dierent data nature and dimensions. In geosciences, volume meshes are complex data made of heterogeneous elements (from structure to properties). Their reduction to better handle their growing dimensions and therefore reduce the computation time of reservoir ow simulations is already a research subject. Still, the variety of models and the assessment of quality at dierent steps of a simulation workow still requires attention.

We demonstrated in the third chapter that compression can signicantly reduce the quantity of scientic data by studying diverse meshes and dierent compression tools, with approaches ranging from lossless to lossy. In the spirit of upscaling methods deployed in geosciences, we base our generic hexahedral mesh representation on HexaShrink, a multiscale decomposition oering dyadic intermediate resolutions. It is principally non-expansive: lower resolutions are embedded in a single structure that does not increase the number of binary objects. It is based on dierent kinds of discrete wavelet decompositions, adapted to each mesh component. We take a special care in the visual preservation of discontinuities, like faults. The relevance of this scheme for visualization is veried on eight meshes, in a comparison with lower-resolution grids obtained from known geomodellers.

HexaShrink allows a sparser representation of regular data variations. Combined with several generic encoders, this hierarchical decomposition is thereafter proved ecient for storing meshes in a lossless manner. Their size can be decreased by a factor of two to ten, while faithfully preserving their exact information. A detailed analysis for each mesh component shows huge dierences in compression ratios. Contrary to structure data, continuous properties were hardly compressible. This is corroborated by the literature, and related to the dynamic range of petrophysical data with classically-used scientic oating-point formats.

This observation encouraged us to pursue our work on compressing continuous properties with an evolved progressive coder, named zerotree. It aims at better handling multiscale decompositions by exploiting remaining redundancies in the transformed data. Besides being adapted to spatial renement across resolutions, it also parses data in a binary-depth order. Contrary to the lossless methods used before, data can be generated at renable precision, allowing a progressively lossy compression, from the most to the least signicant bits.

The quality of reconstructed data was evaluated on a wide range of bit-per-cell values, using objective metrics and visual clues. A focus was laid on continuous properties of lundi, a reservoir model designed for the purpose of a second comparative benchmark. We demonstrated that a wavelet transform associated to zerotree coding was again competitive in performance with respect to generic encoders. It saved at best about ten bits per value at single oating-point precision (32 bits). It however requires a strict control of degradation, to be used with condence in simulation. Indeed, in the fourth chapter, we demonstrated that compression may used very little negative impact on simulation, if properly guided. It is analyzed on the lundi reservoir mesh across a cus-

Perspectives

111 tomizable base of simulation parameters. Motivated by the diculties to compress permeability data, we complemented the wavelet/zerotree compression with compandor, a scalar data transformation inspired from physical principles. We also proposed a novel family of compandor-modied objective metrics. Although simple, they proved well-suited to objectively measure the subjective quality of well production evaluated on compressed meshes.

To be able to conduct the above proof-of-concept for a compression scheme with renable resolution and precision, the HexaShrink decomposition was chosen isotropic. Expecting increased performance with the knowledge of reservoir data orientation, and guidance on the actual precision needed for accurate simulation, we nally explored two complementary directions. The rst one uses anistropy information obtained from variograms, and suggested that sparser decompositions (hence better compression) could be achieved with wavelet packets instead of wavelets. The second one consists in investigating the binary activity of bit depths in zerotree encoding. Looking at the binary precision of the data, we may detect when encoding has already attained a limit of structured information that can be used to compress data. Below that limit, the binary bit planes add little to no information. They thus carry a noise-like content, and could be discarded as unuseful.

These preliminary results however remain solid intuitions, showing promises to be conrmed in future works.

Perspectives

We proposed a comprehensive assessment for the validation of compression methods to be used in a simulation workow. It included lundi, a versatile reservoir model, a combination of subjective evaluations and novel objective metrics. The multiscale HexaShrink representation is usable throughout a complete simulation workow (visualization, lossless and lossy compression, model simulation). It provided overall positive performance with respect to state-of-the-art tools. However, while solid, this work could be complemented in the following directions.

The dataset for comparative evaluation could be extended. The lundi mesh could be lled with dierent property distributions, additional topological discontinuities and especially increased in volume 1 . This motivates our goal to share those datasets with the simulation community, to allow other independent benchmarks to be run.

By studied other volume meshes and dierent property natures, we could further investigate the use of companding as a useless preprocessing step, and the associated compandor metrics for objective/subjective assessment.

As geoscience meshes are composite, a nal compressed bitstream for all the objects should be obtained. This requires to implement progressive coding on every component: activity, corners, pillars, discrete properties, etc. Complexity may arise from the interdependence of the latter: depending on a rock type, and its thickness, its continuous altitude might be quantized at lower precision, without aecting notably the main simulation features. Being able to exploit jointly their impact has potential in improving compression ratios.

We have obtained limited results regarding the acceleration of simulation: only the rst lower resolutions from HexaShrink can be used for upscaling. We may expect improvements by embracing anisotropy, and notably dissociate the decomposition along depth. This would increase the 1 We use a limited size of 128×128×32 for eciency reasons, to be able to performed a large quantity of simulations.

x

  lors de l'encodage zerotree (bit rate per bit plan). Ce qui nous permet de valider expérimentalement notre hypothèse initiale, à savoir : qu'il existe une quantité minimale d'information capable de rendre compte d'une donnée scientique. Cette intuition fondamentale testée sur un nombre réduit de cas pratiques nécessiterait davantage de tests. Pour conclure nous avons démontré dans cette étude l'ecacité d'HexaShrink : un outil de décomposition multi-échelle développé pour améliorer la visualisation d'une donnée de grandes dimensions et faciliter son traitement : transmission, stockage. La méthode de décompositon basée ondelettes est complétée d'outils de compression progressive évolués, capables de générer une donnée scientique à précision réduite. Produite à partir d'une quantité d'information réduite, nous avons montré qu'elle était susante pour les besoins de la simulation. Il est envisagé d'en faire usage dans de futurs travaux an d'optimiser un worow de simulation HPC.

Figure 1 :

 1 Figure 1: Decorating my living room while writing.

Figure 1 . 1 :

 11 Figure 1.1: π number represented with IEEE-754 oating-point standard, equal to 3.1415927 based on single precision.

1. 2 . 5 Figure 1 . 2 :

 2512 Figure 1.2: On left, a picture of the IFP Energies nouvelles supercomputer, called ener440, at Lyon -Solaize, with details of its back (right picture) and its cooling doors equipped with cold water circuit 13 degrees.

Chapter 1 .

 1 General introductionsky to focus on stars and objects of the cosmos.In genomics, technological advances have reduced the cost of DNA sequencing. From an estimated cost of 100 M$ in 2000 for the sequencing the entire genome of a human, this cost decreased to 1000$ starting in 2016. This lower cost has made it possible to study more and more cases. As example, project of United Kingdom's Biobank plans to sequence the genomes of 500 000 volunteers and follow them for decade to give a more complete picture of society-wide health according to[START_REF] Pavlichin | The quest to save genomics: Unless researchers solve the looming data compression problem, biomedical science could stagnate[END_REF]. However, the data generated by the sequencing of an human genome represents ten to hundreds of gigabytes of data. Such gures demonstrate why this eld promises to be one of the most expensive in terms of storage.Changing of scale: Digital Sky Survey (DSS) illustrates the scientic data consumption considering the astronomy. It is a huge atlas of sky pictures started in 2000, [which] collected more data through its telescope in its rst week than had been amassed in history of astronomy (according to the Economist). Obviously, quantity of data for scientic usage did not increase progressively, but exponentially. Considering the current world's largest digital sky survey (PAN-STARS, 2019), it collected 1.60 PB of data during the four rst years, which represents 30 000 times the total text content of Wikipedia"3 .Far from being reassuring, the largest simulations currently generate up to several petabytes of data in a single run, as the example of the code developed by NCAR [National Center for Atmospheric Research] for analysis of climatology presented in Subsection 4.1.1. The size of its outputs (CMIP [Coupled Model Intercomparison projects] current version of Phase6) has been multiplied by one million in the space of 20 years from 0.50 TB to hundreds of petabytes in 2020.

Figure 1 . 3 :

 13 Figure 1.3: Classical compression/decompression pipeline.

Figure 1 . 4 :

 14 Figure 1.4: Lena (512x512 grayscale image) generated by linux batch conversion, from original (on left) to jpeg versions, by keeping 10% (middle) and 5% (right) of the quality, and reducing size of image from 163 Kbytes to 4 Kbytes

Figure 1 . 5 :

 15 Figure 1.5: Classical objective quality metrics organized around the norms in their relative form Lr p or absolute form L p .

Figure 2 . 1 :

 21 Figure 2.1: Illustration of composite geoscientic workow. First step (left) is the geologists work.

(2. 2 )Figure 2 . 2 :

 222 Figure 2.2: Illustration of geological structures, with uvial deposit environment (left) and schematic geological cross section of basin, and resources localization (right).

(2. 3 )

 3 In addition, certain ranges of values have a strong impact on ow simulation. As small ones, whose space accumulation can create impermeable structures. Their placement is inherited from the time of sediment deposition.Deposit environments example Depending on the deposit environment, very dierent structures could be initiated, as illustrated with uvial example on left Figure2.2. Dierence between two pictures of left-gure is the river water level, leading to distinct deposits at the same place. In the upper left-gure, low water ow digs a channel and distributes the coarse and light sediments according to the velocity felt in dierent parts of the channel. In case of ood illustrated in bottom left Figure, water exits the riverbed and oods the neighboring plain, called ood plain. A ne deposit of sediments will form there the future silt and alluvium formations. These rocks have low porosity & permeability and constitute natural impermeable border along channels. Made of coarser elements, channels system is composed by sand formations with high porosity & permeability. It therefore forms a preferential ow way in the subsoil in case of uvial system.

Figure 2 . 3 :

 23 Figure 2.3: Example of a GM used (left); same mesh with the associated porosity property (right).

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: CAD model dened by an unstructured VM.

Figure 2 . 6 -Chapter 2 . Reservoir modeling context 1 node with 1 vertex 1 node with 2 equal vertices 1 node with 4 equal vertices 1 node with 8

 2628 (a) provides an illustration of a fault-free volume. On Figure 2.6-(b), we see that this structure allows to describe for instance a vertical fault (by positioning vertices dierently about the node), while preserving the Cartesian grid. The most popular data structure for structured hexahedral meshes with geometrical discontinuities is the Corner Point Grid tessellation of an Euclidean 3D volume: a structure developed by Ponting (1989) and still currently used (Røe and Hauge, 2016; Lie, 2016). This structure is often termed pillar grid. It is based on a set of vertical or inclined pillars running from the top to the bottom of the geological model. A cell is dened by its 8 adjacent vertices (2 on each adjacent pillar, equal vertices 1 noeud à 4 sommets de position P et 4 sommets de position P' (a) Free-fault area. 4 vertices at position P and 4 vertices at position P' (b) Area with a vertical fault.

Figure 2 . 6 :

 26 Figure 2.6: A fault-free and a fault area.

Figure 2 . 7 :

 27 Figure 2.7: An hexahedron, according to the pillar grid structure.

1 . 1 .

 11 Such parameters have a huge amplitude range (a) Single. (b) Two opposed. (c) Two adjacent.

Figure 2 . 8 :

 28 Figure 2.8: Degenerate hexahedral cells due to a single collapsed pillar (a) or two dierent collapsed pillar locations (b), (c).

Figure 2 . 9 :

 29 Figure 2.9: Teaser for data volume processing. Transformations are successively applied, to condense information essence (approximation) in red subpart, whose lower resolutions are built from Res.-0 to Res.-3.

  VMs multiresolution decomposition based on wavelets[START_REF] Jacques | A panorama on multiscale geometric representations, intertwining spatial, directional and frequency selectivity[END_REF] was proposed by[START_REF] Boscardín | Wavelets bases dened over tetrahedra[END_REF] for tetrahedral meshes. It is based on the tetrahedron subdivision scheme[START_REF] Bey | Tetrahedral grid renement[END_REF] that transforms a tetrahedron into 8 sub-tetrahedra, by introducing 6 new vertices on each edge. After analysis, the input mesh is replaced by a base tetrahedral mesh, and several sets of detail wavelet coecients. Although coecients corresponding to dierences between two resolu-3.1. A variety of representations for volume scientific data35 tions could be encoded for mesh synthesis, this work does not provide an actual compression scheme. Chizat (2014) proposed a prototype for a multiresolution decomposition of geoscientic hexahedral meshes with the pillar grid structure (cf. Subsection 2.1.2). His main contribution resides in a multiresolution analysis (MRA) that partially manages geometrical discontinuities representing the faults. It can be achieved by using a morphological wavelet transform (cf. Subsection 3.2.1). This non-separable transform enables the preservation of some fault shapes at dierent resolutions, as depicted in Figure 3.1.

Figure 3 . 1 :

 31 Figure 3.1: Dyadic non-separable multiresolution rendering on a simple geological mesh Chizat (2014).

(

  [START_REF] Cunningham | Unifying linear dimensionality reduction[END_REF], data clustering[START_REF] Geist | A survey of high-performance computing scaling challenges[END_REF]) and data indexation (Kraska et al., 2018). Recently Li et al. (2018) and Duwe et al. (2020) published surveys on compression of scientic data in the context of HPC, which indicates a growing interest for this topic. Note that a major driving application for compression concerns climatology-related simulations (meteorology, oceanography), where the data generated tend to be the heaviest ones (cf. Subsection 4.1.1).

Figure 3 . 2 :

 32 Figure 3.2: A |L|-level decomposition of a geological mesh obtained by MRA (mesh#6 from Figure 3.5). From the initial mesh (top-left), an approximated version of the initial mesh (S L ) plus |L| sets of details (D i ) are obtained. The number of sets of details increase with the level of decomposition |L|, but the global quantity of data remain stable.

2 3 -Figure 3 . 3 :

 333 Figure 3.3: Analysis and synthesis stages for VMs.

  Subsection 2.1.2), vertices are inevitably positioned along pillars. So, our multiresolution scheme for geometry information only focuses on: the z coordinates of the 8 vertices associated to each node. According to the naming convention presented in Figure 3.4, those 8 vertices can be dierentiated according to their relative positions [Back (B)/Front (F), Bottom (B)/Top (T), Left (L)/Right (R)]; the x and y coordinates of the nodes describing the low (bottom) and high (top) extremities of all the pillars (the x and y coordinates of intermediary nodes being implicit). The nodes are called hereinafter the oor and ceil nodes, respectively.The GVMs can exhibit very irregular boundaries. Hence, a Boolean eld called Actnum may be associated to each cell to inactivate its display (and its inuence during simulations as well). It (a) A node and its 8 surrounding cells. view of the node into its 8 vertices.

Figure 3 . 4 :

 34 Figure 3.4: Vertex naming with the Back (B)/Front (F), Bottom (B)/Top (T), Left (L)/Right (R) convention.

Figure 3 . 5 :

 35 Figure 3.5: Mesh#5 (in yellow) has inactive cells (in red) to describe its boundaries using the Actnum eld.

Figure 3 . 6 :Figure 3 . 7 :Figure 3 . 8 :

 363738 Figure 3.6: HexaShrink multiresolution scheme for geometry: (left) input grid and its top view; (middle) output from the non-separable, non-linear 2D morphological wavelet based on a fault segmentation (based on the top view); (right) non-linear 1D wavelet transform along pillars (orange lines).

Figure 3 . 9 :

 39 Figure 3.9: Prediction of a fault node at resolution l -1 from the four parents' conguration at resolution l, orange ovals denoting ∨ operands.

Figure 3 .

 3 Figure 3.11: Our benchmark composed by eight meshes.

Figures 3 .

 3 Figures 3.12 and 3.13 show respectively the decompositions of mesh#1 and mesh#8 provided by HexaShrink. The dierent resolutions are arranged in rows, by decreasing scale. The rst column represents the mesh without any attribute. The second and the third columns represent the same mesh onto which a continuous and a categorical property is mapped, respectively.

Figure 3 . 54 Chapter 3 .

 3543 Figure 3.13: Original mesh#7, its attributes, and four levels of resolution generated with Hexa-Shrink.

55 structure is based

  on a combination of LZ77 and Human algorithms. bzip2, developed four years later, is in direct competition with gzip, and uses the Burrows-Wheeler transform[START_REF] Burrows | A block-sorting lossless data compression algorithm[END_REF], taking advantage of recurring patterns, and nishing also by an Human algorithm.

Figure 3 .Figure 3 .

 33 Figure 3.15 presents the execution times of our conservative compression workow in function

Figure 3 .

 3 Figure 3.16: Binary cost of each component on mesh#5 in function of the number of successive HexaShrink decompositions.

Figure 3

 3 Figure 3.17: A textbook case with a 1D vector to show the interest of encoding transformed/decomposed data at renable precision. (Left) the original vector coded/decoded at renable precision. The index O refers to the original data. (Right) a multiscale version of the original vector (obtained by transformation) coded/decoded at renable precision. The index M refers to the multiscale data.

LL LH 2 LH 1 HL 2 HH 2 HL 1 HH 1 Figure 3 .

 2122113 Figure 3.18: Inheritance between the subbands generated by wavelet-based decomposition on 2D and 3D grids.

Figure 3 .n e a r s h o r e 0 n e a r s h o r e 1 n e a r s h o r e a fl u v iFigure 3 . 21 :

 31321 Figure 3.19: 2D data at three dierent resolutions (Haar wavelet) and rened numerical precision: from the original with 8 × 8 values in 031 with 5MSB M (top-left) to a two-fold coarser resolution coded with 2MSB M (bottom-right).

Chapter 3 .Figure 3 .

 33 Figure 3.22: SNR curves of four compression workows. Three of them are conservative (LZMA, HS P +LZMA, HS P +ZT in lossless mode), while Haar+ZT is used in progressive mode. The marks are obtained thanks to a progressive decoding of the bit planes. The four red circles in the bottom-right graphic refer to visual representations of Figure 3.23.

Chapter 3 .Figure 3 . 23 :Figure 3 .

 33233 Figure 3.23: Original uvial permeability property generated at renable precisions by using lower bit budgets of data (preprocessed by Λ ) processed by Haar+ZT combination. Quality is evaluated here by SNR metric. Precisions noted from 1 to 4 correspond to the red circles on Figure 3.22.

  Currently, the simulation domain the most aected by issues in data management is climatology. Unsurprisingly it is the most interested by compression, given the growing number of publications. Simulations made in the domain are meant to understand past, present and to forecast future climates. Studies mentioned in the review are aliated to well-established computing centers: National Center for Atmospheric Research (NCAR) in Denver, German Computing Center (Deutsches Klimarechenzentrum, DKRZ) in Hamburg and the weather machine at the Los Alamos National Laboratory, New Mexico. Yet, the majority of the examples of simulation data provided in open-access are from NCAR. Among their projects, we can cite ASR [Artic System Reanalysis], and POP [Parallel Ocean Climate] (initiated by Los Alamos Laboratory and revised by NCAR) meaning for compression. The project the most often encountered in this literature review is unquestionably the CESM-LE [Community Earth System Model-Large Ensemble] simulation code. The generated data are open-source, and the results of its analysis are discussed during IPCC Intergovernmental Panel on Climate Change (GIEC [Groupe d'experts intergouvernemental sur l'évolution du climat] in french), which sets every four or ve years. This gathering of scientists aims at alerting public opinion on climate change.

77 4 terabytes

 4 for a thickness of about 2 cm. To store the largest version of CMIP Phase6 it would be necessary to pile ten thousand disks of 4 terabytes, with a stack height of 200 meters. There are several versions of CMIP (or side projects: CESM-ATM, ASR, CAM-SE). The set represents the most tested data to assess compression. Allison Baker, notably, is one of the rst researchers who investigated the subject of lossy compression in climatology. She tested varied tools and encouraged the development of methods inspired by multimedia. Lossy compression tools (such as fpzip, ISABELA, GRIB2, wavelet-zerotree) were used in Baker et al. (2014) to question the relevance of over-precision when transmitting, storing and analyzing data. In this article, she notably underlines the sensitivity of complex scientic codes. With dierent supercomputers, CESM code will not give the same bit-for-bit results". Therefore, it seems relevant to ask what is interesting in preserving the full precision of the data. Because it appears that limited changes were tolerable, results are subjectively evaluated by experts. It motivates Baker et al. (2017) to pursue this line of research looking to dene acceptability by using the subjective opinion of climatologists. SWOT [Surface Water and Ocean Topography mission] deals with experimental data registered by satellites (NASA, CNES) to provide information on large rivers, lakes, and reservoirs. Although the data are not generated by simulation, we mention his mission because it is a current project of large dimension. Moreover it served to assess compression performance of generic lters available on netCDF format (Delaunay et al., 2019).

2 :

 2 1 recently constituted to assess compression methods. The initiative in the context of Exascale Computing Project came mainly from Argonne and Livermore laboratories teams, leading actors on the subject. 1 https://sdrbench.github.io/ Chapter 4. Impact on simulation performances Simulation codes, mini-app 4.1.2 Impact of compression on application The usage of lossy compression through simulation workows became popular in the last ve years. This practice knows an increasing success with the introduction of new lossy compression tools for scientic data, rstly with ZFP (Lindstrom, 2014), followed by SZ (Di and Cappello, 2016; Liang et al., 2018b; Tao et al., 2017), already introduced in Subsection 3.1.3. Behind these projects are research teams on compression from Lawrence Livermore National Laboratory (LLNL) and Argonne National Laboratory (ANL), led respectively by Peter Lindstrom and Franck Cappello. They produced a large number of publications to test and compare performances of compression tools on scientic data from diverse domains. As an evidence of the growing interest for the topic, since 2016, more than a dozen of graduate students joined ANL team to work on SZ projects, to apply and adapt the method to various steps of simulation for dierent domains. They contributed to publishing more than twenty papers at various IEEE Conferences on Big Data -for High Performance Computing, Networking, Storage and Analysis -Parallel Distributed processing Symposium etc. mainly dedicated to HPC. The reference article of Cappello et al. (2019) compiles seven lossy case studies and points the main simulation bottlenecks (listed in Subsection 1.2.2, and denoted in italics in the rest of the section).

2 .

 2 (2018) adapted SZ into SZ-PaSTRI. They improved SZ prediction by pattern identication. In the same vein, Liang et al. (2018b) considerably ameliorated compression performance on cosmology simulation (HACC) by reducing storage footprint. They combined SZ with a logarithm mapping transform and changed the tool name into SZ_vlct. The method was nally incorporated into SZ 2.0 (selecting point wise relative option). In heavy simulation, checkpointing consists in creating periodically temporary les to backup intermediate steps in case of simulation failure. Lossy compression is there used to reduce le size, in order to relieve the overall process without degrading simulation results. This simulation step has been the subject of many works from leading laboratories, testing the impact of lossy compression on temporary les. In the example of Calhoun et al. (2018), at ANL, SZ is used in PlasComCM code to reduce the overall checkpointing time. Most recently, Underwood et al. (2020) developed FraZ to optimize existing tools (SZ, ZFP, MGRAD). Since compression codes are rich in options, performance varies depending on option selection To reduce these dierences, FraZ optimizes the xed rate option, which is less ecient, compared to the xed accuracy. Aliated to ANL, we note that simulation data used for the work are shared at the Scientic Data Reduction Benchmark (such as Hurricane, HAAC, CESM, NYX).From LLNL, Dienderfer et al. (2019) theoretically established the bound of the error introduced by ZFP integrating operators during its compression. Continued by Fox et al. (2020), the method theoretically addresses forward and backward error for simulation and validates the process by using PDEs that model diusion and advection phenomena. We note that Lindstrom et al. (2016)collaborated on an isolated work related to geology. In this eld, compression methods are familiar because they have been historically used on seismic data (but barely for simulation, as already said at the end of the Subsection 4.1.1). This time, instead of applying the method to the experimental data, the compression is applied to an operator used for full 3D seismic waveform tomography (f3dt) simulation.Beyond the studies carried out by LNL and ANL, SZ and ZFP are commonly used as benchmarks in other laboratories to compare their performance with their proper compression tools. Zhang et al.(2019), for instance, compare SZ & ZFP to a method based on block decomposition supplemented by DCT and followed by adaptive quantization. They demonstrate the viability of their solution to restart a simulation from a lossy state, by checking error propagation. The method is applied with success to CMIP5 and Nek5000. As seen with the previous example, the rest of the section deals with tools based on wavelet transforms or using similar decompositions to integrate it into their workow. SSEM, for instance, is a wavelet-based method developed by[START_REF] Sasaki | Exploration of lossy compression for application-level checkpoint/restart[END_REF] during a collaborative work between Tokyo Institute of Technology and LNL to deal with checkpoints for a climate application: NICAM. In another study, Nek5000 is again used by[START_REF] Otero | The eect of lossy data compression in computational uid dynamics applications: Resilience and data postprocessing[END_REF]

  3.3). It appears resourceful in this case to combine such evaluations with more costly physical analysis methods (introduced by Laney et al. (2014) for instance) and check their correlation. For example, to determine the optimal wavelet kernel for turbulent ow visualization, Li et al. (2015) in a comparative study evaluate classical error metrics such as nRMSE in parallel with the enstrophy. This parameter, in uid mechanics, computes the kinetic energy by considering the volumes of the highest values. In the same vein, Pulido et al. (2016) compare the impact of dierent multiresolution transforms (B-splines, Daubechies, Coiet, curvelet, surfacelet etc.) and evaluate alterations on turbulence datasets by considering vorticity, isosurface, curvature etc. in parallel to classical PSNR.

  others study the domain requirements to propose tailored solutions. More recently, exascale projects in simulation emerged, giving rise to the need to process mass-produced data. The quantity of data to be produced in the future encouraged the development of compression tools for scientic data and the creation of specic benchmarks, highlighting the advent of a new discipline marked by the necessity of rationalizing the simulation data. Considering this review, our work will be between the works of Laney et al. (2014) and of Baker et al. (2019).

Chapter 4 .Figure 4 . 1 :

 441 Figure 4.1: (Left) illustration of a quarter ve-spot model. Blue dotted lines depict motion of water injected from corners by four wells noted I i up to producer P 1 on the top of anticline. Our model lundi (quarter of anticline) is highlighted in brown and limited by vertical red planes. It only considers two wells (I 1 and P 1 ) among ve. (Top-right) water cut curve measured over a year at P 1 .

Figure 4 . 2 :

 42 Figure 4.2: Objective and subjective evaluation steps throughout our compression-simulation workow.

Chapter 4 .Figure 4 . 3 :

 443 Figure 4.3: Subjective evaluation of simulation results, using region boundaries around the water cut curve (around WC curve from nearshore 0 environment) to dene a qualitative ranking: 'Identical' (in blue), 'Correct' (in green), 'Acceptable' (in yellow), 'Uncorrect' (in orange) and 'Aberrant' (in red).

  87ences between the cell values of the original property denoted P (c) (c being the indexation of the cells in the mesh) and the cell values lled by a rened property denoted P (c). These dierences are generally averaged, and can be considered relative to the original values or squared etc., generating a large number of potential metrics, previously introduced in Figure1.5 in Subsection 1.3.3. However, they might be ill-suited to petrophysical properties. In the following, we only concentrate on a limited number of objective metrics and propose a novel family of quality indicators combining classical metrics and compandor. The classical metrics we select are: MRE, nRMSE and SNR (respectively on Equations 4.1, 4.2 and 4.3).

(4. 3 ) 2 C

 32 The novel family of metrics is obtained by appending the compandor operator lambda to the classical metrics. For instance, the Λ-SNR equation becomes: c=1 Λ(P (c)) -Λ( P (c))

(4. 4 ) 4 . 2 . 2

 4422 HexaShrink evaluation in upscaling/upgridding The HexaShrink representation can perform up to 5 levels of decomposition on the lundi reservoir model. The rst three resolutions lled by nearshore 0 properties (permeability and porosity) are shown on Figure 4.4. Structurally, the faults remain visually coherent across resolutions. The uniformly reduced number of cells across resolutions is detailed in Table 4.3. Properties results of HexaShrink decompositions for our four environments are represented in Figures 5.11 and 5.12 in appendices (page 126).

Chapter 4 .Figure 4 . 4 :Figure 4 . 5 :

 44445 Figure 4.4: (Top-row, left to right) reconstruction of lundi and its rst three resolutions obtained with HexaShrink. (Center and bottom) the same resolutions lled with the permeability and porosity data, themselves decomposed.

Λ 0 Λ 1 Figure 4 . 6 :

 146 Figure 4.6: From original nearshore 1 permeability (top) are decompressed two versions of the data at renable precision, both using 18 MSB M processed by compandor: (left) with Λ 0 or (right) with Λ 1 . Data are objectively, subjectively assessed and classical λ compandor objective metrics. (*)

Figure 4 .

 4 Figure 4.7 summarizes the MRE, SNR, nRMSE, and Λ-SNR curves we obtained for nearshore 1 permeabilities. Dashed curve refers to Λ 1 (linear compandor) while the solid curve denotes the Λ 0 use (logarithm compandor). The bit budget depends on the precision targeted during the (de)decompression, corresponding to a number of bit planes decoded. These graphs also indicates

Figure 4 . 8 :Figure 4 . 9 :

 4849 Figure 4.8: Permeability from nearshore 1 (left) and fluvial (right) environments generated at renable precision decreasing the bit budget (a mark per precision) using our method (compandor), objectively evaluated by Λ-SNR. Subjective appreciation of simulation results is represented by a color code.

4. 2 )

 2 . We use their recent versions, and compile them with their most ecient options according to (Underwood et al., 2020): SZ, 2.1.3 version, with PSNR compilation option; ZFP, 0.5.5 version, with accuracy compilation option. Results are displayed on Figure 4.10 and complement those presented in Figure 4.8. Dashed lines represent compression results for the regular use of the tools, while solid lines assess the eect of Λ 0

Figure 4 . 10 :

 410 Figure 4.10: Permeability from nearshore 1 (left) and fluvial (right) environments generated at renable precision decreasing the bit budget (a mark per precision) using our method, SZ and ZFP (with or without compandor stage), objectively evaluated by Λ-SNR. The subjective appreciation of the simulation results are still represented by the color code.

  tutes a promising research axis for scalable representations in simulation workows, including all the components of a RM. The recent paper of Hoang et al. (2021) actually oers an unied encoding of resolution and precision for scalar elds.

Figure 4 .

 4 Figure 4.12: Variogram curve in one direction, for an exponential model (among spherical, gaussian, cubic etc. models).

Figure 4 .

 4 Figure 4.13: 1D variogram curves of the permeability property in the nearshore 1 environnement.

Figure 4 .

 4 Figure 4.14: Dyadic decomposition (middle) vs wavelet packet decomposition along k (right). The initial data (left) is the permeability property in nearshore 1 environment preprocessed by Λ 0 . The weighted entropies for each data is indicated in the bottom part of the gure.

Chapter 4 .

 4 Impact on simulation performances the asymptotic number of bits per symbol required to transmit the minimal quantity of information in a lossless compression perspective. For a data constituted by C discrete values, each value c has a probability of occurrence equal to P c . The entropy, in bits per value, is dened by H = -C c=1 P c log 2 (P c ).

Figure 4 .

 4 Figure 4.15 in black thus complement the graphs of Figure 4.10.

Figure 4 .

 4 Figure 4.15: Permeability from nearshore 1 (left) end fluvial (right) environments generated at ren-

Figure 4 .

 4 Figure 4.16: Correlation between ZT activity (top-curve) and objective evaluation by Λ-SNR (bottom-curve) of permeability from nearshore 1 environment generated at renable precision by increasing bit planes number for reconstruction while our compression alternative using Λ 0 .

Figure

  Figure 4.17: Correlation between ZT activity (top-curve) and objective evaluation by Λ-SNR (bottom-curve) of permeability from nearshore 1 environment generated at renable precision by increasing bit planes number for reconstruction while our compression alternative using Λ 1 .

Figure 5 . 4 : 2 ( 5 2 ( 5 E 1 Figure 5 . 5 :

 542525155 Figure 5.4: Water cut curve in blue obtained on the original data: near shore environement, serves as reference. The water breakthrough is spotted by rst red point.

Figure 5 . 6 :

 56 Figure 5.6: Subjective evaluation of simulation results, using region boundaries around the water cut curve (around WC curve from nearshore a environment) to dene a qualitative ranking: 'Identical' (in blue), 'Correct' (in green), 'Acceptable' (in yellow), 'Uncorrect' (in orange) and 'Aberrant' (in red).

Figure 5 . 7 :

 57 Figure 5.7: Subjective evaluation of simulation results, using region boundaries around the water cut curve (around WC curve from fluvial environment) to dene a qualitative ranking: 'Identical' (in blue), 'Correct' (in green), 'Acceptable' (in yellow), 'Uncorrect' (in orange) and 'Aberrant' (in red).

Subsection 4 . 2 . 2 :Figure 5 . 8 :Figure 5 Figure 5 . 10 : 1 Figure 5 .

 42258551015 Figure 5.8: Overlay of water cut curve simulated on ne grid lundi mesh, (nearshore 1 environment)and WC simulated on its lower resolutions (up to res-3) generated by HexaShrink.

1 Figure 5 .Figure 5 .

 155 Figure 5.12: Visualization of the lundi lower resolution from HexaShrink lled by continuous properties, with porosity on upper part and permeability on lower part, for the four distinct geological environments: "nearshore a " (left), "fluvial " (right).

Figure 5 .

 5 Figure 5.14: Three variograms of the mesh directions computed on continuous property on nearshore a environment.

Figure 5 .

 5 Figure 5.15: Three variograms of the mesh directions computed on continuous property on fluvial environment.

Subsection 4 . 2 . 3 :Figure 5 .

 4235 Figure 5.16: Permeability from nearshore 0 (top on the left) end nearshore 1 (top on the right) nearshore a (bottom on the left) and fluvial (bottom on the right) environments generated at renable precision for decreasing bit budget (a mark per precision) by our alternative (wavelet packet or dyadic decomposition), SZ and ZFP (using or not Λ 0 ), objectively evaluated by Λ-SNR; highlighted by subjective appreciation of simulation results (color code).

Figure 5 .

 5 Figure 5.17: Permeability from nearshore 0 (top on the left) end nearshore 1 (top on the right) nearshore a (bottom on the left) and fluvial (bottom on the right) environments generated at renable precision for decreasing bit budget (a mark per precision) by our alternative (wavelet packet or dyadic decomposition), SZ and ZFP (using or not Λ 0 ), objectively evaluated by MRE; highlighted by subjective appreciation of simulation results (color code).

Subsection 4 . 3 . 2 :Figure 5 . 133 Figure 5 .

 43251335 Figure 5.18: Correlation between ZT activity (top-curve) and objective evaluation by Λ-SNR (bottom-curve) of permeability from nearshore a environment generated at renable precision by increasing bit planes number for reconstruction with our compression alternative using Λ 1 .

Figure 5 . 135 Figure 5 . 21 :

 5135521 Figure 5.20: Correlation between ZT activity (top-curve) and objective evaluation by Λ-SNR (bottom-curve) of permeability from nearshore 0 environment generated at renable precision by increasing bit planes number for reconstruction with our compression alternative using Λ 1 .

Figure 5 . 137 Figure 5 . 23 :

 5137523 Figure 5.22: Correlation between ZT activity (top-curve) and objective evaluation by Λ-SNR (bottom-curve) of permeability from fluvial environment generated at renable precision by increasing bit planes number for reconstruction with our compression alternative using Λ 1 .

  Using dynamic of reservoir property . . . . . . . . . . . . . . . . . . . . . . . Description of HexaShrink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	xiv xviii			Contents Contents
	2.1.3 Compression vocabulary 2.2 From storage to simulation HexaShrink, from mesh representation to coding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . compression ratio CR levels of details LOD MRA multiresolution analysis	23 26 29
	3.1	A variety of representations for volume scientic data	. . . . . . . . . . . . . . . . .	30
	3.1.1 3.1.2 3.1.3 Preprocessing Upscaling & upgridding techniques Compression of VMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Compression of scientic data . . . . . . . . . . . . . . . . . . . . . . . . . . Λ compandor expandor Λ CHAPTER 1 30 31 35
	. . . . . . . . . . . . . . . . . . . . . Scientic projects with oversized data . . . . . . . . . . . . . . . . . . . . . . Infobesity: source variety & growing production 1.2.1 1.2.2 Bottleneck & issues in data workow . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Energetic, economic and ecological cost . . . . . . . . . . . . . . . . . . . . . Numerical diet for (scientic) data? . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Legacy of multimedia standards . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Data compression components . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 How taking advantage of these three techniques for a dream representation of GVM? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusions & perspectives 1.2 1.3 3.2 5.1 Meshes & components VM volume mesh GM geologic mesh GVM geologic volume mesh RM reservoir mesh CAD computer-aided design P generic continuous property Objective evaluation metrics Lr relative euclidean distance MRE mean relative error maxRE maximum relative error L euclidean distance MAE mean absolute error maxAE maximum absolute error MSE mean square error RMSE root mean square error nRMSE normalized root mean square error SNR signal noise ratio compandor-signal noise ratio Λ-SNR Tools Simulation HPC high performance computing CFD computational uid dynamic WC curve water cut curve Conclusions . . Acronyms 4 4 5 6 7 8 8 9 39 109 discrete wavelet transform DWT CohenDaubechiesFeauveau CDF 9/7 zerotree coder ZT coder General introduction
	1.3.3 2 Reservoir modeling context Classical performances evaluation peak signal noise ratio PSNR SSIM structural similarity index measure . . . . . . . . . . . . . . . . . . . . . . . .	10 15
	Binary writing MSB LSB	most signicant bit least signicant bit		20
			xiii xvii

1.1 Data feeds scientic workow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.1 Data types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2 Digital infrastructures for data handling . . . . . . . . . . . . . . . . . . . . 2.1 From reservoir formation to numerical model . . . . . . . . . . . . . . . . . . . . . . 17 2.1.1 A primer in geosciences for reservoir . . . . . . . . . . . . . . . . . . . . . . . 17 2.1.2 VM for reservoir modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2.1 Multiresolution scheme for geometry . . . . . . . . . . . . . . . . . . . . . . . 42 3.2.2 Multiresolution scheme for properties . . . . . . . . . . . . . . . . . . . . . . 47 3.3 Visual results: decompositions obtained with HexaShrink . . . . . . . . . . . . . . 49 3.4 Conservative compression workow for GMs . . . . . . . . . . . . . . . . . . . . . . . 51 3.4.1 Coding performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.4.2 Speed performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4.3 In-depth analysis of the coding performance . . . . . . . . . . . . . . . . . . 57 3.5 Compression workow at renable precision for GMs . . . . . . . . . . . . . . . . . . 59 3.5.1 How coding at renable precision? . . . . . . . . . . . . . . . . . . . . . . . . 59 3.5.2 Principle of zerotree coding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.5.3 Coding performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Impact on simulation performances 73 4.1 Compression in simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.1.1 Application for uid dynamic . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.1.2 Impact of compression on application . . . . . . . . . . . . . . . . . . . . . . 80 4.2 Evaluation of multiresolution progressive compression method . . . . . . . . . . . . 83 4.2.1 Proposed uid simulation workow . . . . . . . . . . . . . . . . . . . . . . . 83 4.2.2 HexaShrink evaluation in upscaling/upgridding . . . . . . . . . . . . . . . . 87 4.2.3 Continuous properties, which precision? . . . . . . . . . . . . . . . . . . . . . 90 4.3 Deepening the method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.3.1 Exploiting the anisotropy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 4.3.2 Thresholding at necessary precision . . . . . . . . . . . . . . . . . . . . . . . 102 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 "Data is the new oil" Clive Humby, 2006

  Its results were transmitted and analyzed by searchers in local on their desktops. In 2011, for a bandwidth capacity of 1 MB s -1 , the transfer of a data sets of size between 56 GB and 5.60 GB (daily production) could catch from 12 hours to 50 days, according to the report of[START_REF] Woodring | Revisiting wavelet compression for large-scale climate data using JPEG 2000 and ensuring data precision[END_REF].

	1.2. Infobesity: source variety & growing production	7
	storage/memory size of supercomputers equipment. The latter has indeed increased tenfold over
	the last ten years, while bandwidth has only improved its capacity by a factor of 2.5. Therefore, as
	shown by Cappello et al. (2019), the HPC performance of the largest supercomputers is currently
	limited by data transfer during simulation, which creates bottlenecks (as for checkpointing/restart
	approach explained in Subsection 4.1.2). We could also mention the need for in situ visualisation
	(mainly studied), code execution acceleration etc. improved using recent compression tools later
	discussed in our literature review Subsection 4.1.2.
	These bottlenecks mentioned in HPC studies are proofs of the limitation of a digital ecosystem
	facing the rise of data production. Their deleterious eect is also perceptible at a more global level
	by studying energetic, economic and ecological indices.

As example to illustrate past transmission diculties, we can cite a simulation code, named POP [Parallel Ocean Program] (Smith and Gent, 2002), developed by Los Almos National Laboratory and run on supercomputers.

Although the bandwidth capacity has been steadily increased (reaching 2.50 TB s -1 for IBM Summit, the world's second largest supercomputer), it did not improve enough compared to the 3 https://hubblesite.org/contents/news-releases/2019/news-2019-12.html

  [START_REF] Kunkel | Exascale storage systems an analytical study of expenses[END_REF] estimated at this time at 144 000e per year the price of storing 5.60 PB of data dedicated to earth Science system. Among processing methods classically studied to save space, he concluded compression is the most promising for the laboratory needs.

	Chapter 1. General introduction
	1.3 Numerical diet for (scientic) data?
	Given the growing increase in data production in HPC, scientists address the question by various
	research axes. They notably explore approaches based on compression, deduplication, indexation
	(Kraska et al., 2018) and clustering (Geist and Reed, 2016). Advances in those domains are carried
	out by prominent research centers, for instance in the USA (Lawrence Livermore National Lab., Oak
	Ridge National Lab., Argonne National Lab. National Center for Atmospheric Research (NCAR)),
	in China or Japan (Sakai et al., 2013; Kolomenskiy et al., 2018). Large HPC facilities are being
	developed, notably with exascale computing, able to perform billiars ops. In particular, it is being
	promoted at the International Conference for High Performance Computing, Networking, Storage,
	and Analysis, on November 2020. The processing of their outputs has given rise to numerous studies.
	More recently, the emergence of bottlenecks during simulation, mentioned in Subsection 1.2.2, has
	led to new questions to overcome them. Scientists seek methods that can be integrated into the
	scientic workow (in situ technologies) by satisfying divers applications that comprise it.
	So far we pointed several bottlenecks along the simulation workow considering problematic data
	volumes, but situation is worst than we believe. The eld suers from a lack of standards, especially
	in experimental and simulation sciences. Therefore a universal solution is not practicable. Conversely
	multimedia data (pictures, videos, sound) have classical formats, leading to the development of
	compression standards that have invaded our daily lives. From there we identify a classical scheme
	of compression tools, and common methods to evaluate their eciency.
	1.3.1 Legacy of multimedia standards
	Compression facilitates the dissemination of information by reducing data size. It is worth noting
	that multimedia compression supported the deployment of the web. Interestingly, the most used
	standards for 1D signals (audio, known as MP3) and 2D images (pictures, with JPEG) were already
	developed by the beginning of the 1990's. In audio and pictures, individual digital les are relatively
	small. Standards have evolved more steadily for 3D multimedia data instantiated by video (space
	and time). While storage capacity and network transmission have witnessed an impressive increase
	over the last decade, data quantity (including duplication) grows exponentially, and disk transfer
	rates tend to plateau. For instance, a Full HD (high-denition) two-hour movie (1080 lines, 1920
	pixels per line, 50 images per second) needs a transmission capacity of 2.49 Gbit s -1 for streaming
	without compression. The entire video le requires 2.49 × 3600 × 2, or 17.91 Tbit for storage! There
	is no need to mention that video streaming would be impossible without compression. Two types of
	compression are traditionally distinguished: with perfect information preservation, or allowing some
	precision loss.

Various studies in particular led by DKRK (Deutsches Klimarechenzentrum -German Climate Computing Center) (Hübbe et al., 2013; Kunkel et al., 2014; Kuhn et al., 2016) seek comparing various compression tools, and economical return by storage optimization.

  . Its eciency has been investigated in many works, and notably on RMs[START_REF] Rasaei | Upscaling and simulation of waterooding in heterogeneous reservoirs using wavelet transformations: Application to the SPE-10 model[END_REF]. The basic concept consists in performing a wavelet transformation on the permeability property, producing blocks of details. Independently, their values are compared to a pre-dened threshold to determine which block is considered as signicant during the uid ow. If a block is considered as non signicant, the cells constituting this block are merged. Otherwise, the cells are preserved, to conserve the full resolution in the regions of interest for simulation. Later, Babaei (2013) em-

ployed wavelets to provide a coarsening operator during simulation. Hence, the simulations are driven by adapted grids, which reduces the computation times, while preserving satisfactory simulation results. A more recent article

[START_REF] Rezapour | Upscaling of geological models of oil reservoirs with unstructured grids using lifting-based graph wavelet transforms[END_REF]

generalizes the method to irregular grids by using graphs, as later explained in Subsection 3.1.3 with the work of

[START_REF] Iverson | Fast and eective lossy compression algorithms for scientic datasets[END_REF]

.

  Chapter 3. HexaShrink, from mesh representation to coding values. Hence, it becomes possible to encode each coordinate with an integer index, instead of a 32-bit oating-point value. It is common to quantize the coordinates with 12 or 16 bits, reducing the geometry information by a compression factor of 2.60 or 2 respectively. Quantization inevitably

	introduces an irreversible loss in accuracy. Visualization typically tolerates precision loss (as long
	as visual distortion remains negligible), unlike some numerical simulations requiring a priori more
	precise computations.
	Prediction (as well as related interpolation methods) further improves the geometry compactness.
	Predictive coding resorts to estimating the position of a vertex from already encoded neighbor
	vertices. Prediction errors (dierences between predicted and actual positions) are generally smaller
	in amplitude and sparser, which makes their entropy coding (which codes dierently frequently
	occurring patterns) ecient (Salomon and Motta, 2009, p. 63 sq.). Regarding the connectivity,

To reduce the memory footprint or make the transmission of VMs faster, well-known techniques exist, for instance the quantization of vertex coordinates. It consists of constraining the vertex coordinates to a discrete and nite set of

  [START_REF] Isenburg | Compressing hexahedral volume meshes[END_REF] are the rst to deal with hexahedral VMs. The connectivity is encoded as a sequence of edge degrees in a way similar to Touma and Gotsman (1998) for triangular meshes via a region-growing process of a convex hull called hull surface. It relies on the assumption that hexahedral meshes are often highly regular, which implies that the majority of vertices are shared by 8 cells. It involves an almost constant edge degree all over the mesh, which signicantly decreases the entropy of the connectivity information. A context-based arithmetic coder[START_REF] Witten | Arithmetic coding for data compression[END_REF] is then used to encode the connectivity at very low bit rates, between 0.18 and 1.55 bits per hexahedron. Regarding geometry, a user-dened quantization rst restricts the number of bits for coordinates, and then a predictive scheme based on the parallelogram rule encodes the position of vertices added during the region-growing process. To globally optimize the geometry compression, a Minimum Spanning Tree minimizing the global prediction error for the whole mesh is computed. This method is more ecient than prior ipping approaches whose traversal does not depend on the geometry, but solely on the connectivity.

	Streaming compression is a subcategory of single rate compression, dedicated to huge data that
	cannot t entirely in the core memory. A particular attention to I/O eciency is thus required,
	to enable the encoding of huge meshes with a small memory footprint. Isenburg and coworkers

[START_REF] Krivograd | A hexahedral mesh connectivity compression with vertex degrees[END_REF]

propose a variant to

[START_REF] Isenburg | Compressing hexahedral volume meshes[END_REF]

that encodes the vertex degrees number of non-compressed hexahedra around a given vertex instead of the edge degrees. On the one hand, this variant achieves better compression performances than

[START_REF] Isenburg | Compressing hexahedral volume meshes[END_REF]

for dense meshes. On the other hand, it only deals with manifold meshes, and the algorithm is complex as interior cells are encoded after border cells (it involves many specic cases to process when encoding the connectivity).

Lindstrom and Isenburg (2008) proposed an original algorithm for unstructured meshes called

Hexzip. This algorithm is considered as fully lossless, because the initial indexing of vertices and hexahedra is preserved. For this purpose, connectivity is encoded directly in its indexed structure, by predicting the eight indices of an hexahedron from preceding ones. This technique is suitable because hexahedral meshes generally have regular strides between indices of subsequent hexahedra. A hash-table is then used to transform the index structure into a very redundant and byte-aligned list of symbols, that can be compressed eciently with gzip

[START_REF] Deutsch | DEFLATE compressed data format specication version 1[END_REF]

. Concerning the geometry, spectral prediction

[START_REF] Ibarria | Spectral predictors[END_REF]

is used. This algorithm is faster and less memory intensive than

[START_REF] Isenburg | Compressing hexahedral volume meshes[END_REF]

as the connectivity is not modied. It handles non-manifold meshes and degenerate elements. On the other hand, bitrates are higher because of the lossless constraints. Unlike methods presented above, Chen et al. (2005) focus on geometry compression for tetrahedral meshes. The authors proposed a ipping approach based on an extension of the parallelogram rule (initially proposed for triangle meshes (Touma and Gotsman, 1998)) to tetrahedra. It consists in predicting the position of an outer vertex of two face-adjacent tetrahedra, with respect to the other vertices.

are the rst to propose streaming compression for VMs (extended from his method for triangular meshes ): for tetrahedral meshes

[START_REF] Isenburg | Streaming compression of tetrahedral volume meshes[END_REF]

, and then for hexahedral meshes

[START_REF] Courbet | Streaming compression of hexahedral meshes[END_REF]

  Chapter 3. HexaShrink, from mesh representation to coding Lossless techniques Generic lossless methods are used by broad audience to compress a wide As already noted in the context of lossless compression, it is fairly common to test popular lossy compression methods on scientic data although not adapted. JPEG and JPEG2000 developed for multimedia data have been for example tested on climatology by[START_REF] Hübbe | Evaluating lossy compression on climate data[END_REF][START_REF] Woodring | Revisiting wavelet compression for large-scale climate data using JPEG 2000 and ensuring data precision[END_REF][START_REF] Baker | A methodology for evaluating the impact of data compression on climate simulation data[END_REF], and uid mechanics[START_REF] Schmalzl | Using standard image compression algorithms to store data from computational uid dynamics[END_REF].Computational uid dynamics (CFD) is powerful tool for simulating turbulent ows. However, computational time can be high, due to the large scale of simulated data. To address this problem,[START_REF] Nakahashi | High-density mesh ow computations with pre-/post-data compressions[END_REF] proposes a compact support called BCM (Building Cube Method) to decrease the number of cells. It assemblies sub-cubes (Cartesian grids) at variable resolutions to create an unstructured grid, with local renement adapted to geometry and ow features. To diminish the size of binary data and relieve the computation, run-length coding is applied to the BCM grids.The work of[START_REF] Sasaki | Exploration of lossy compression for application-level checkpoint/restart[END_REF] improves this method by applying the lter CDF 9/7 to each cube.Wavelet coecients are then quantied at various thresholds, and then zerotree coded (the method is presented in Subsection 3.5.2).[START_REF] Bradley | Wavelet transform-vector quantization compression of supercomputer ocean models[END_REF] worked on wavelet based solution with quantization step for NCAR ocean models. Its dimensions were limited, but already problematic to repatriate the simulated data from supercomputers. Hereafter open source libraries as QccPack Fowler (2000) permit simple implementations of methods based on wavelets. Chapter 3. HexaShrink, from mesh representation to coding progressive lossy compression by truncation of the less signicant bits. The impact of such a reduction is notably studied on dierent simulation codes (LULESH: Lagrangian shock hydrodynamics, Miranda: High-order Eulerian hydrodynamics, details in Subsection 4.1.1) by Laney et al. (2014) as detailed in Subsection 4.1.2. ISABELA (In-situ Sort-And-B-spline Error-bounded Lossy Abatement) by Lakshminarasimhan et al. (2011) provide an adaptation to noisy data. To prepare it, a smoothing spatial preprocessing is applied. Then, the data is approximated using B-splines or wavelets. Additionally, the temporal dimension is also exploited by temporal pattern identication. However, the decompression can not be performed on the entire data, but only on subsets, which complicates its usage. Iverson et al. (2012) develop a compression model for unstructured meshes in geosciences, by conversion of the grid into a graph, taking advantage of the locality. In such a structure, the graph is dened by N nodes of the mesh, the set of edges E connecting two nodes. A zero weight indicates if the connection does not exist, otherwise a weight proportional to the connection importance is used. It can correspond to a distance if considering the mesh geometry, but can deal with multiple other properties. In two years, Ainsworth et al. (2018, 2019, 2020a,b) produce four sequential works on multilevel techniques for compression and reduction of scientic data (MGRAD project). From univariate to multivariate cases, Ainsworth et al. (2020b) nally manages 2D and 3D meshes. This method employs an orthogonal decomposition based on the Riesz basis properties. In this study, airfoils are studied to evaluate preservation of aerodynamic force and coecients of pressure. Other applications are mentioned in Subsection 4.1.2 to evaluate the impact of compression in simulation. .1. A variety of representations for volume scientific data 39 et al. (2019) are interested in intrinsic metrics of both compression tools, to better understand processes. In addition to netCDF and HDF, new libraries specic to lossy compression emerge for large scientic data. For example, SCIL (Scientic Compression Library) (Kunkel et al., 2017) notably integrates ZFP and SZ. Several compression components can be combined according to the needs of users. CubismZ, created by Hadjidoukas and Wermelinger (2019) is another recent example. This library proposes a parallel implementation of methods based on wavelets, or tools such as ZFP, SZ and fpzip. Finally, lossy compression has been often employed in the context of simulation, but researchers remain cautious about the risk of data degradation. Beyond the size reduction, compression promotes discussions toward a rationalization of ever growing datasets, on terms of required precision and discretization. Even if compression did not seem to pervade the eld of simulation in geosciences (except in Iverson et al. (2012) and Lindstrom et al. (2016)), an awareness on the need to manage data at a reduced precision already exists, if considering the simplication brought by upscaling.

	variety of data, regardless of data features (dimension, type, ...). Mainstream methods although
	generic, such as gzip or LZMA (tested in Subsection 3.4.1), may be ecient in term of compression
	ratios and execution times. Nevertheless, lossless methods dedicated to scientic data are now
	available.
	FPC (Floating-Point Compression) is a fast and lossless algorithm conceived to optimize the
	transmission and the storage of oating-point numbers Ratanaworabhan et al. (2006). Based on
	prediction techniques, FPC compresses better and faster than other algorithms available at the time
	of its release. But its eciency decreases with the randomness rate of the data.
	Later, many scientists, particularly climatologists, adopted two standard formats to store multi-
	dimensional data in an unied library interface: NetCDF (Network Common Data Form) and HDF5
	Lossy compression allows achieving higher compression gains, but implies to control the data
	degradation strictly. However, there is no metric established in the context of simulation, because Coming back to uid mechanics, we can mention other works using wavelets on various structured
	expectation vary according to the application. In this case, quality metrics borrowed from the simulation data, from 2D to 4D (Wilson, 2002; Schmalzl, 2003; Kolomenskiy et al., 2018) or point
	visualization eld are found sucient, and are used to appreciate the distortion involved by lossy clouds (Salloum et al., 2018). Loddoch and Schmalzl (2006) proposed an equivalent approach for
	compression (maximum point wise error, root mean squared error (RMSE), peak signal to noise compressing volume uid dynamic datasets by replacing wavelets by other transform and such as
	ratio (PSNR)). DCT, B-splines, etc. and predictors such as Lorenzo, etc.

(Hierarchical Data Format 5)

[START_REF] Koranne | Handbook of Open Source Tools[END_REF]

). Yet, in this format the original data are not available anymore, and must be decompressed entirely to be used. HDF5 provides compression components to develop tools designed specically for climatology, such as MAFISC

[START_REF] Hübbe | Reducing the HPC-datastorage footprint with MAFISCmultidimensional adaptive ltering improved scientic data compression[END_REF]

.

A sequence of lters is applied to the data and completed by generic encoding tools.

Other compression formats have been developed specically for HPC. For instance, ISOBAR

[START_REF] Schendel | ISOBAR preconditioner for eective and high-throughput lossless data compression[END_REF]

proposes to identify and lter hard-to-compress datasets in order to better target the compression needs.

Lossy techniques Lossless methods ensure a perfect reconstruction of the initial data, but their compression performance is limited. Despite this limitation, many researchers are reluctant to use lossy compression, fearing that the losses will bias or distort simulation results. However, as highlighted ironically by

[START_REF] Kipnis | Analog-to-digital compression: A new paradigm for converting signals to bits[END_REF]

, data discretization in simulation is already a form of quantization, which leads to a lost of information.

Today, popular formats such as HDF5 progressively incorporate powerful and easily implementable lossy compression tools (Linear Packing, Layer Packing (Silver and Zender, 2017), Bit shaving, Bit Grooming (Zender, 2016), Digit Rounding (Delaunay et al., 2019) etc.) that remove last bits of values considered as noise and/or irrelevant data. These irreversible and drastic components are completed by lossless techniques that rearrange bits (Shue) and coders (DEFLATE, Zstandards). The results obtained by Delaunay et al. (2019) on meteorological and oceanic experimental datasets made of oating values are promising. Their degradation is positively assessed by objective metrics (SNR, mean relative error, mean absolute error). Lossy compression is now a common place in climatology (Baker et al., 2017, 2019), where the data volume for the forecast are the most massive. One of our oldest references in simulation data compression ((Bradley and Brislawn, 1993)) already dealt with storage of data generated by supercomputers from NCAR. In Clyne et al. (2007), VAPOR enables compression and visualization of massive computational datasets. It also facilitates analysis thanks to its interface. This work has been motivated by the data complexity of two particular applications: turbulent plume dynamics and current sheet formation. VAPOR is based on a wavelet decomposition (mainly Haar) that generates a hierarchical representation. Last version of VAPOR (Li et al., 2019) guarantees that the visualization package [is] tailored to analyze simulation data in earth system. To further investigate lossy compression based on wavelets, Li et al. (2015) worked on a volume data set characterizing turbulent ows. The eciency of several wavelets (Haar, CDF 9/7, CDF 8/4) are evaluated for visualization, with classical objective and proper eld metrics: the best option seems to be a combination of the lter CDF 9/7 and coecient prioritization (instead of approximation). The method has been also implemented on spatio-temporal data (Li et al., 2017b).

Wavelet based methods inspired many works on scientic data compression in various elds, from seismic in geosciences to tomography in medical. Most of the works deals with experimental data, but very few with simulation data. However, in recent years, simulation elds tended to catch up. As mentioned above with our oldest reference on the topic, we observed that the interest of integrating compression in a simulation workow is not new.

Fpzip

[START_REF] Lindstrom | Fast and ecient compression of oating-point data[END_REF]

) is an online compressor for oating-point values, implementable on the I/O of diverse applications. It manages various data through its multidimensional prediction scheme based on the Lorenzo predictor. Depending on the usage, the scheme proceeds After fpzip, Lindstrom (2014) introduces ZFP, available on HDF5, ADIOS library. The concept is notably inspired by xed-rate texture compression methods used in graphics hardware. It permits a random access to compressed oating-point data, by splitting into blocks at variable precision. Floating-point representation is exploited by dierentiating the exponents and the mantissa in IEEE-754 format. DCT is applied to the mantissa elements within subblocks, and the resulting coecients are zerotree coded. The exponents are distinctly stored. A current work aims at determining a closed form expression for bounds on the error introduced by the three compression options of ZFP (Dienderfer et al., 2019). SZ (Cappello et al., 2019) is a predictive method for multidimensional oating point data. A rst step consists in attening the data into a single dimensional array. If appropriate, a log-mapping transform can be applied (in the version 2.0 (Liang et al., 2018a)). The array is then read progressively, and the values are predicted from the previous ones, looking for the best tting among various deterministic models (mean-integrated Lorenzo predictor, curve tting scheme constant, linear or quadratic). Quantization, Human coding and gzip complete the compression workow. Discussion SZ and ZFP appear today as the most competitive compression tools for scientic data. Complementary tools are trying to exploit their shortcomings. For example, FRaZ (Underwood et al., 2020) works on improvement of user modes (xed-rate, absolute error bound etc.) and Wang

3

  Yes 22.73 MB 46 % Porosity, Permeability 6 524,288 128 × 128 × 32 Yes 64.27 MB 100 % Porosity, Permeability 7 5,577,325 227 × 95 × 305 Yes 274.57 MB 97 %

	3.4. Conservative compression workflow for GMs	51
	Mesh index # Cells 1 93,600 80 × 45 × 26 No 4.62 MB 100 % Characteristics Dimension Faults File size Actnum 2 1,000,000 100 × 100 × 100 No 42.46 MB 100 % 3 36,816 59 × 39 × 16 Yes 1.46 MB 100 % 4 210,000 100 × 100 × 21 Yes 7.88 MB 20 % 5 450,576 149 × 189 × 16 Porosity Properties Continuous Porosity 8 13,947,600 240 × 295 × 197 Yes 580.94 MB 100 % Porosity	Categorical Rock type Rock type Rock type
	Table 3.1: Our collection of geological meshes with their ontological characteristics and geological
	properties.

diagonal crest shape) and a faulty block on the right. Even at the coarsest level, corresponding to a downsampling by 2 4 × 2 4 × 2 4 , these two structural discontinuities are still present, while keeping a good shape delity, globally. Concerning the attributes, the decompositions are also adequate.

Table 3 .

 3 2 presents the compression ratios obtained with gzip, bzip2, and LZMA applied to the outputs of HexaShrink. Each mesh component is encoded independently. A compression ratio (CR) is the ratio between the size of the raw data and the size of the compressed/encoded data. The higher the compression ratio, the more ecient the compression/encoder. The lines none indicate

the CRs obtained when the raw data are directly encoded (HexaShrink is disabled). The lines with a number i indicate the CRs obtained when the mesh components are decomposed i times

Table 3 .

 3 Chapter 3. HexaShrink, from mesh representation to coding 2: Coding performance of our conservative compression workow. HexaShrink is combined with gzip, bzip2 and LZMA, for dierent levels of decomposition.

	Mesh	Level	gzip	bzip2	LZMA
		none	3.73	4.98	6.43
	1	1	5.62	6.07	7.52
		24	5.67	6.126.13	7.427.44
		none	3.23	8.41	10.12
	2	1	6.49	10.82	11.81
		26 7.487.58 12.7513.03	13.35
		none	2.67	2.99	3.63
	3	1	3.88	4.70	5.24
		24 4.034.05	4.924.93	5.475.48
		none	1.83	1.89	2.21
	4	1	2.64	3.06	3.48
		24	2.76	3.223.23	3.643.65
		none	2.46	2.55	3.33
	5	1	3.14	2.83	3.71
		24 3.253.26	2.912.92	3.803.81
		none	1.88	2.25	3.04
	6	1	2.70	3.17	3.71
		26 2.842.86	3.393.42	3.903.93
		none	2.32	2.25	3.04
	7	1	3.31	3.53	4.44
		26 4.144.24	4.484.68	5.545.73
		none	3.20	5.98	12.52
	8	1	5.42	7.07	8.90
		27 5.806.72 7.6310.12 9.0510.23

Table 3 .

 3

	3: Comparative (near) lossless coding performances on porosity (top tabular) and perme-
	ability (bottom tabular) property from our four environments with compression ratios at dierent
	decomposition levels combining HexaShrink with LZMA, zerotree, and Haar wavelet transform
	with ZT.

  .2.

		The reference for well production is the WC curve
	Original RM	Analysis Renable compression * P Synthesis	RM at renable resolution/precision
	P	Objective metric	P
		Simulation workow
	WC	Subjective metric	WC

Table 4 .

 4 3: Dimensions of lundi and its lower resolutions generated by HexaShrink.

	4.2. Evaluation of multiresolution progressive compression method	89
	Resolution	Fine scale model size	Total cell number
	0	128 × 128 × 32	524288
	-1	64 × 64 × 16	65536
	-2	32 × 32 × 8	8192
	-3	16 × 16 × 4	1024
	-4	8 × 8 × 2	128
	-5	4 × 4 × 1	16

  Λ-SNR (top-left), MRE (top-right), SNR (bottom-left), nRMSE (bottomright). Subjective appreciation of simulation results is represented by a color code.

		Compression methods
						Λ 1,34bits CFD 9/7 -ZT				Λ 0,34bits CFD 9/7 -ZT
		Validity water cut curve
					Identical			Correct			Acceptable		Uncorrect		Aberrant
											80												0.5
											60												0.4
											40 45	Λ-SNR											0.2 0.3	MRE
											20												0.1
																						0.03
	20	18	16	14	12	10	8	6	4	2 100 0	0	20	18	16	14	12	10	8	6	4	2	0 4	0 •10 -2
											80												3
											40 60	SNR											2	nRMSE
											20												1
	20	18	16	14	12	10	8	6	4	2	0	0	20	18	16	14	12	10	8	6	4	2	0	0
			bit budget (bits per value)						bit budget (bits per value)
	Figure 4.7: Permeability from nearshore 1 generated at renable precision for decreasing bit budget
	(a mark per precision) by our alternative (changing compandor parameter), objectively evaluated
	by four quality metrics:

Table 5 .

 5 2: Weighted entropies computed on permeability property of four environments decomposed by dyadic wavelets or wavelet packets.

			Weighted entropy
		Volume data	Dyadic wavelet transform	packets wavelet transform
	nearshore 0	18.48	15.58	14.97
	nearshore 1	18.92	15.61	14.99
	nearshore a	18.92	15.61	14.99
	fluvial	18.81	15.62	15.00

https://www.statista.com/statistics/254266/global-big-data-market-forecast/

https://www.businesswire.com/news/home/20190627005451/en/Global-Big-Data-Market-Witness-CAGR-19.7

https://www.grandviewresearch.com/press-release/global-high-performance-computing-hpc-market

https://www.marketsandmarkets.com/Market-Reports/Quantum-High-Performance-Computing-Market-631.html

[START_REF] Salloum | Adaptive wavelet compression of large additive manufacturing experimental and simulation datasets[END_REF] Now all these market growth predictions are subject to change upwards in the light of recent events (majority of references edited in 2019, pre-Covid19).

It is, therefore, necessary to have a good knowledge of the tool to obtain the best compression: what is dicult to ask from users of non-computer science backgrounds.

Remerciements

Pour ce travail de thèse, j'ai pu bénécier de l'aide et du soutien de nombreuses personnes appar-

the handling of cell blocks C of 2×2×2 cells throughout scales, we use an adaptation of the well-known Haar wavelet. The resolution l -1 is a scaled average of cells at resolution l. The approximation coecient p l-1 is thus the average value of the related eight property coecients {p l 1 , p l 2 , . . . , p l 8 }.

The seven details required for synthesis are dierences with respect to the approximation coecient:

p l n ; d l-1 n = p l n -p l-1 , ∀n = 1 .

To deal with real-valued (oating-point) properties, and avoid accuracy imprecision due to the divide operator, we introduce the following modications. First, reals are mapped into integers up to a user-dened precision, here with a 10 6 factor. Second, we disable the division by using a sum.

The analysis system thus becomes:

and the synthesis system turns into:

Approximation and coecients are stored as is. To recover the accurately scaled values, the division operator should however be applied as a simple linear post-processing.

Categorical properties We complete the global mesh multiresolution decomposition with an original categorical-valued scheme called modelet [START_REF] Antonini | Method of exploitation of hydrocarbons of an underground formation by means of optimized scaling[END_REF]. We assume that a mesh cell category belongs to a set of classes Ω 0 = {ω 1 , ω 2 , ..., ω W }, taking discrete values. The cell block

Visual results: decompositions obtained with H S

49

C l = {p l 1 , p l 2 , . . . , p l 8 } thus contains, at resolution l, integers indexing categories from Ω l . They take values in a subset of Ω. The multiresolution scheme is expected to produce, at lower resolutions, discrete values in embedded subsets:

It may happen that the above denition does not yield a unique maximum. If two or more categories dominate a cell block, a generic approach consists in taking into account its rst block cell neighborhood (the surrounding 26 cells, except at mesh borders and boundaries). We aect the dominant value in the rst neighborhood to p l-1 . In case of a draw again, the second-order surrounding can be used, iteratively. In practice for the presented version of HexaShrink, we limit to the rst-order neighborhood, and choose the lowest indexed category when the maximum is not unique. Equipped with this unique lower resolution representative value, we proceed similarly to Subsection 3.2.2 for details, by using dierences between original categories and the mode. As classes are often indexed by positive integers, a slight motivation allows to get only non-negative indices. By avoiding negative values, one expects a decrease in data entropy of around 5 %, which benets to compression.

We thus change the sign of a detail coecient if and only if it generates a value out of the range of {ω 1 , ω 2 , . . . , ω W }, and then control this condition during reconstruction. So, all details {d l-1 n } for a cell block C are determined by:

).

During synthesis, the coecients {p l i } are obtained thanks to the closed-form equation:

HexaShrink performs a decomposition of the mesh into a pyramidal structure, the inverse process perfectly reconstructs the data from the lower scale to the initial mesh, integrating details progressively. Adapted transforms are distinctly implemented on the various components, dealing with their dimensions and features as previously explained. The tool has rstly been designed for a visualization purpose, using thus multimedia concepts. Later, positive compression results demonstrated the HexaShrink capability as a compression tool. The decomposition generates sparser data. Then, encoding steps permit to reduce the data quantity by reordering the data, to a more compact storage size. The prime contribution consists in analyzing the process and experimenting the method on various meshes, in order to validate an appropriate generic approach.

Visual results: decompositions obtained with HexaShrink

Geosciences oer a huge variety of meshes, because of their various dimensions, structures and properties. We propose to validate HexaShrink on a representative benchmark of eight meshes.

By this way we expect to experiment a generic methodology, suitable for various meshes.

Chapter 3. HexaShrink, from mesh representation to coding trees mostly composed of zero or close to zero coecients. By property inheritance, an insignicant coecient in a particular subband is likely to have insignicant descendants in related subbands.

The principal mechanism for data coding consists in localizing signicant root coecients above a given threshold. Then bit plane per bit plane, descendant coecients below the threshold are coded globally as zerotrees. The progressive compression process starts with a threshold close to the maximal coecient. When a wavelet bit plane (corresponding to amplitudes between 2 q-1 and 2 q -1) is encoded, the threshold is reduced to describe the lower bitplane. By rst coding signicant features, and progressively thereafter smaller details, this approach better takes into account the inherent wavelet decomposition structure than external lossless coders. The transmission of a ZT coded data could be interrupted at any time, resulting in a (de)compression at renable precision.

The truncated data would correspond to a more or less accurate approximation of the initial data. This is lossless/conservative when all bit planes are transmitted.

Combining a ZT coder with HexaShrink seems to be a promising way to encode GMs in a progressive way, to nally provide a compression/decompression workow at renable precision, and at renable resolution. Aside pure compression gains, we are interested in the potential of ZT to contribute to simulation eciency. Indeed, works are currently devoted to nding alternatives to Because the mesh#6 is complete and has ideal dimensions for the reservoir simulations presented in the next chapter, unlike other meshes 1 , all our tests with the ZT coder will be performed on this mesh. In addition, we only present the CR obtained on continuous properties, as we showed in the previous section that these data are the most delicate to compress.

The continuous properties porosity and permeability of the mesh#6 are inspired by the tenth SPE 1 Seven other meshes in the benchmark are sedimentary basins and have kilometric spatial dimensions. The average distance of inuence between two wells is 300 m, so the dimensions of the basins are much too large to consider simple cases.

Chapter 3. HexaShrink, from mesh representation to coding

Comparative Solution Project (also known as SPE10) [START_REF] Christie | Tenth SPE comparative solution project: A comparison of upscaling techniques[END_REF]. We experimented two formations (cf. Figure 3.20): a Tarbert formation, prograding near shore environments, and a Uperness formation, prograding uvial environments.

Chapter 4. Impact on simulation performances

However this care does not seem to signicantly improve simulation results.

As observed, changes in mesh resolution are not simple. Such manipulations can signicantly impact the simulation results. The HexaShrink inherent structure constrains possibilities for upscaling improvement, yet with a simple dyadic decomposition, lower resolutions may provide consistent results and accelerate simulation. In the following, we pursue our quest of essential data required for simulation, we explore the impact of numerical precision and improve the compression ratios.

Continuous properties, which precision?

In the previous chapter (cf. Section 3.5), we showed that a ZT coder oers promising compression performance on the continuous properties, in comparison with generic encoders. We now pursue this methodology to investigate the impact of renable precision on simulation outcomes. In Section 3.5.3, we show that the quality of data compressed with our workow at renable precision is evaluated by SNR between 150 dB and 200 dB (see Figure 3.22). It is therefore questionable whether such a quality is adequate and necessary for scientic applications, given that the values dynamic are much higher than those encountered in multimedia. If sucient, we can further go into compression using lower numerical precisions. Then we will try to dene a minimum precision threshold, i.e., the minimum binary quantity (bit budget) for suitable simulation.

Renable representation with zerotree

Our method is applied to permeability properties. As other scientic data, they are hardly compressible because of their distribution. A specicity of permeability is a large range of possible values (typically from 0.0007 to 20 000 mD as SPE10 model) with a large proportion of very small values.

We especially investigate the use of preprocessing permeabilities with compandors (cf. Subsection 2.1.3) noted Λ α . We compare the performance using objective metrics as well as subjective validation through the simulation workow. Our principle observations are summarized in italic hereinafter.

On Figure 4.6 are displayed three versions of nearshore 1 permeability at dierent precisions. The original one and two others only using the eighteen MSB M of data decomposed by CDF 9/7 wavelet [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF]) and then encoded with ZT. They dier from the use of a compandor, applied before transformation as a preprocessing. The version on the left uses Λ 0 (logarithmic compandor) while the right one uses Λ 1 (linear compandor). We observe that the left version is better preserved, for instance around the blue patch with low values. On the right picture, the same patch appears noisy with the use of Λ 1 . It concurs that the lower values are better preserved by compression if using Λ 0 . Besides, for a similar MSB M , the overall bit budget is lower with Λ 0 , but also the quality is better, visually. This demonstrates that a ZT coder is more ecient on this preprocessed data, better decorrelated by the wavelet transform when the compandor Λ 0 is included in the compression workow.

It appears Λ 0 increases the compression performance of a ZT, and better preserves lower values, visually.

Deepening the method

We present in this section partial investigations and intuitions that appear to be promising for optimizing the RM compression. The anisotropy of our meshes rst inspired us. This feature particularly observable on properties can be studied using a classical statistical tool in geosciences, called variogram. It allows determining preferential orientations in the distribution of petrophysical properties. In the next subsection, we test a packet wavelet decomposition along these orientations, to exploit this feature. Such a method should increase compression performances without degrading the simulation results.

In the second subsection, we focus on the ZT-based encoding stage. The study of its intrinsic parameters provides information on the data structure and could advise on the minimum amount of binary data required for simulation. Their analysis could even be more complete than the objective study performed so far.

Exploiting the anisotropy

For this experimentation, our approach is still evaluated on the continuous properties of the four geological environments, generated on purpose with anisotropic behavior. The distribution of these properties is conditioned by the depositional environment: sediments, initial sedimentary basin topology, sediments transport, diagenesis [START_REF] Fowler | Dissolution/precipitation mechanisms for diagenesis in sedimentary basins[END_REF] as explained in Subsection 3.1.3.

Such parameters are directionally dependent, subject to marine/uvial current during the deposit and more generally by the gravity, the pressure from the upper layers, or other constraints during all its genesis. The sedimentary structure is stratied, composed of horizontal layers of varying thickness and composition. By denition, this structure is thus spatially anisotropic along ( ı, ), with a potentially dierent behavior along a depth ( k), as shown in Figure 4.11. Consequently, considering that these properties are isotropic is probably suboptimal for compression.

Isotropic

Isotropic on surface Anisotropic through the depth Anisotropic in three directions Finally, a major line of research would reside at the interface of simulation and compression.

Simulation codes use complicated renement schemes, potentially on multiresolution grids. They can be sensitive to the data precision, hence the use of extended oating-point formats. They are not, to the best of our knowledge, devised to use the prioritization of information that can be aorded by compression tools: prominent information at specic resolutions, precision that can be rened upon numerical accuracy bounds, etc. Such a blending of compression pipelines and simulation workows could yield gains in both data storage and computational time.

List of Figures During the rst phase, extracted uid does not contain any water. Therefore water saturation is zero or close to zero. There, only tiny constant error E 0 is tolerated by reservoir engineer. During the middle phases, water saturation sharply increases to progressively reach a plateau during last phase. Acceptable errors, respectively noted E 1 and E 2 temporally grows, by cumulation of e i and a delay for two last periods, as reported in the Table 5.1.