Student Member, IEEE Emanuel Zgraggen

Alex Galakatos

Andrew Crotty

Senior Member, IEEE Jean-Daniel Fekete

Tim Kraska

How Progressive Visualizations Affect Exploratory Analysis

Keywords: Exploratory analysis, interactive visualization, progressive visualization, scalability, insight-based evaluation !

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

T HE literature [1], [2], [3], [4] often states that a delay of one second is the upper bound for computer responses after which users lose focus on their current train of thought. In order to ensure a highly interactive environment for data analysis, a visual data exploration system should therefore strive to present some actionable and understandable artifact for possible query over any dataset within a one second threshold. It is important to note that this artifact does not need to be the complete or most accurate answer, but it should be an answer that allows users to keep their attention on their current task.

Traditionally, visual data exploration systems employ a strategy whereby user-issued queries are offloaded to a database management system (DBMS), and the results are displayed once the complete answer is computed. We call this the blocking approach: a user's current train of thought is blocked until the query result is fully computed. Even with modern hardware and state-of-the-art DBMSs that can process millions of data points per second, this traditional approach suffers from the basic issue that some datasets will still be too "big" to yield results for interactive user queries within one second.

A wide variety of strategies, including precomputation, prefetching, sampling, and progressive computation, have been proposed to overcome this fundamental limitation. However, each of these approaches comes with its own set of advantages and challenges. In particular, progressive • E. Zgraggen, A. Galakatos, A. Crotty computation, where data is processed incrementally in small chunks, offers an interesting tradeoff between result accuracy and computation speed. Moreover, unlike many of the other approaches, progressive computation also provides a number of natural opportunities to incorporate user feedback and computational steering. The research community has recently regained interest in progressive computation, attempting to analyze and exploit some of the peculiarities of this approach [START_REF] Fisher | Trust me, i'm partially right: incremental visualization lets analysts explore large datasets faster[END_REF], [START_REF] Stolper | Progressive visual analytics: User-driven visual exploration of in-progress analytics[END_REF], [START_REF] Fisher | Incremental, approximate database queries and uncertainty for exploratory in Large Data Analysis and Visualization (LDAV)[END_REF], [START_REF] Fisher | Exploratory visualization involving incremental, approximate database queries and uncertainty[END_REF], [START_REF] Fekete | Progressive analytics: A computation paradigm for exploratory data analysis[END_REF], [START_REF] Rosenbaum | Progressive refinement: more than a means to overcome limited bandwidth[END_REF]. However, the effects of progressive computation-and progressive visualization-on user behavior and knowledge discovery in exploratory settings have not been studied in detail. The aim of this work is to investigate how progressive visualizations affect users in exploratory settings. To this end, we design and conduct an experiment where we compare three different visualization conditions: (1) instantaneous, (2) blocking, and (3) progressive. The instantaneous visualization condition acts as a stand-in for hypothetical systems where all queries, independent of the dataset size, return and display accurate results within a strict latency constraint. On the other hand, blocking visualizations simulate traditional systems where results are displayed only after the full dataset has been processed. Blocking visualizations are limited by the throughput of the underlying computation engine; that is, the wait time increases proportionally to the size of the dataset. Finally, progressive visualizations, where data is processed incrementally in small chunks, present approximate results to the user at different points during the computation. To compare each of these strategies, we capture interaction logs and verbal data from the participants using a think-aloud protocol. We then extract several knowledge discovery and user activity metrics from these recordings, and we analyze these metrics to test how progressive visualizations compare to the alternatives. For the purpose of our study, we picked simple uncertainty visualizations and update strategies and are specifically not testing different variations in those.

The main contributions of this paper are two-fold. First, we find that progressive visualizations significantly outperform blocking visualizations in almost all knowledge discovery and user activity metrics. Second, our results show that, surprisingly, progressive visualizations do not differ substantially from the best case scenario of instantaneous visualizations across many key metrics (e.g., insights per minute, insight originality, visualization coverage percentage). These findings suggest that progressive visualization is a viable solution to achieve scalability in visual data exploration systems.

RELATED WORK

Our research builds upon related work in the areas of Big Data Visual Analytics and Latency in Computer Systems.

Big Data Visual Analytics

Visual data analysis, or the task of gaining insights from a dataset through visualizations, is an interactive and iterative process where users must frequently switch between a wide range of distinct but interrelated tasks. While the specific set of tasks that recur in visual data analysis, as well as the tools that support them, are relatively well understood [START_REF] Amar | Low-level components of analytic activity in information visualization[END_REF], [START_REF] Heer | Interactive dynamics for visual analysis[END_REF], the constantly increasing volume of data has forced the interaction paradigm away from interactive approaches back to large-scale batch processing [START_REF] Fisher | Interactions with big data analytics[END_REF].

Thus, we seek to address this fundamental conflict in visual exploratory data analysis. On one hand, we want to provide an experience where users can actively steer the exploratory process, allowing them to see results for each action without the distraction of a slow and unresponsive interface. On the other hand, constantly growing amounts of data and increasingly complex analysis techniques make it impossible to compute and deliver accurate responses within latency thresholds that are suitable to keep users on their current train of thought.

The research community has proposed several approaches to address this problem, each with its own set of advantages and challenges. We introduce a simple use case to illustrate these different approaches. Imagine a visual data exploration system that processes tabular data, allowing users to create simple aggregated histograms over this data by selecting different attributes. Furthermore, histograms are linked together, where selections in one visualization trigger filtering operations in others, such as in GraphTrail [START_REF] Dunne | Graphtrail: Analyzing large multivariate, heterogeneous networks while supporting exploration history[END_REF], PanoramicData [START_REF] Zgraggen | Panoramicdata: Data analysis through pen & touch[END_REF], and Vizdom [START_REF] Crotty | Vizdom: interactive analytics through pen and touch[END_REF]. Figure 1 shows an example where the user only wants to see the histogram of income for a specific age range.

Conceptually, the simplest way to implement a data exploration system to support these interactions is through a blocking approach. After the user requests a particular histogram, the system scans the full dataset to perform the required aggregation and displays the result only after processing all of the data points. Using this approach, the user cannot see results until after the entire dataset has been processed. Response times of blocking systems are therefore proportional to the size of the dataset. Several commercially available visual analysis tools, including Tableau [START_REF]Tableau[END_REF] (and its research predecessor Polaris [START_REF] Stolte | Polaris: A system for query, analysis, and visualization of multidimensional relational databases[END_REF]), Spotfire [START_REF]Tiboco Spotfire[END_REF], and Microsoft Power Pivot [START_REF]Microsoft Power BI[END_REF], use a blocking approach.

Precomputation provides an opportunity to reduce the latencies of the blocking approach by performing some of the computation up front before a user begins exploration. The system can perform this precomputation during a loading phase and simply return a result from its cache when the user issues a request. This approach requires a substantial time and processing investment during the initial loading phase but can completely eliminate visualization latency during data exploration. However, a major drawback of precomputation is the potentially enormous number of both visualizations and exploration paths the user can take. In particular, the number of possibilities depends both on the characteristics of the data as well as the exploration tasks supported by the system. Permitting arbitrarily filtered histograms (e.g., as shown in Figure 1) drastically increases the number of required precomputations. Even with these issues, many systems have successfully used precomputation to improve the user experience. For instance, some commercial DBMSs support the creation of online analytical processing (OLAP) cubes [START_REF] Chaudhuri | An overview of data warehousing and olap technology[END_REF], where the data is pre-aggregated along selected dimensions. Other research systems introduce improvements to these techniques using more advanced algorithms and data structures [START_REF] Lins | Nanocubes for real-time exploration of spatiotemporal datasets[END_REF] or by exploiting modern hardware (e.g., GPUs) [START_REF] Liu | imMens: Real-time visual querying of big data[END_REF].

Similar to precomputation, prefetching incrementally precomputes results during user exploration rather than computing all possible results a priori. Prefetching approaches typically either limit the degrees of freedom given to the user or employ an intelligent oracle that predicts the user's most likely subsequent actions, such as in a map application where users can perform only pan and zoom actions. Since the user can only transition to a very limited number of possible next states (i.e., panning left/right/up/down or zooming in/out), the system can therefore use the time that a user spends examining a specific region to prefetch neighboring tiles that might be requested next. Both Fore-Cache [START_REF] Battle | Dynamic prefetching of data tiles for interactive visualization[END_REF] and Semantic Windows [START_REF] Kalinin | Interactive data exploration using semantic windows[END_REF] use prefetching by exploiting locality in the user's exploratory behavior to predictably prefetch chunks of data in anticipation of subsequent. Other prefetching based systems use models to estimate the most likely action the user will perform next. These predictive models can be built using a wide variety of different information. For example, Doshi et. al. [START_REF] Doshi | Prefetching for visual data exploration[END_REF] propose and compare different techniques that incorporate a user's interaction history, while Ottley et. al. [START_REF] Ottley | Personality as a predictor of user strategy: How locus of control affects search strategies on tree visualizations[END_REF] show that certain personality traits affect user exploration strategies.

Prefetching systems need to provide a fallback strategy for cases when the prediction fails, most commonly by reverting to a blocking approach for cases when the user issues a query for which the result has not yet been prefetched. Furthermore, even with prefetching, the result might take too long to compute if the user only spends a short amount of time between interactions.

While all of the approaches discussed thus far are guaranteed to return completely accurate results for all queries, sampling takes a fundamentally different approach by trading speed for accuracy. Instead of computing the completely accurate result, the system uses a small sample of the data to compute the histogram and presents it to the user with some quality metric for the approximation (e.g., confidence intervals, error bars). Systems that use sampling can return results for all user queries within the specified interactivity time constraints without needing to precompute any query results. The data management community has heavily explored the development of approximate query engines. For example, BlinkDB [START_REF] Agarwal | Blinkdb: queries with bounded errors and bounded response times on very large data[END_REF] is a SQL engine that uses sampling to answer queries over large datasets and allows users to specify accuracy or response time requirements. DICE [START_REF] Kamat | Distributed and interactive cube exploration[END_REF] similarly supports subsecond latencies for large datasets using a variety of optimizations including sampling.

Although sampling-based approaches can provide users with a quick overview of the dataset, they also introduce a completely new set of challenges. For instance, rare (but potentially important) datapoints might not be captured in the sample. Additionally, even experts sometimes have trouble interpreting statistical accuracy metrics [START_REF] Cumming | Inference by eye: reading the overlap of independent confidence intervals[END_REF], [START_REF]Understanding the new statistics: Effect sizes, confidence intervals, and meta-analysis[END_REF]. Ferreira et. al. [START_REF] Ferreira | Sample-oriented taskdriven visualizations: allowing users to make better, more confident decisions[END_REF] introduced specialized visualizations that mitigate some of these issues, but further research is necessary in order to apply their findings to different visualization types and analysis tasks.

Progressive systems are an extension of sampling-based approaches that incrementally compute results over increasingly larger samples in order to provide more accurate results to the user over time. This concept is well known in the graphics domain and widely used on websites to improve user experience when loading high-resolution images [START_REF] Harrison | Evaluation of progressive image loading schemes[END_REF]. A down-sampled, low-resolution image is displayed quickly and replaced with the high-resolution version when fully downloaded. This concept was previously explored in the data management community [START_REF] Hellerstein | Interactive data analysis: The control project[END_REF], where most of the subsequent research focused on creating progressive versions of well known DBMS operations (e.g., joins [START_REF] Haas | Ripple joins for online aggregation[END_REF], aggregations [START_REF] Condie | Mapreduce online[END_REF]). More recently, progressive approaches have gained interest among the HCI and visualization communities [START_REF] Fisher | Trust me, i'm partially right: incremental visualization lets analysts explore large datasets faster[END_REF], [START_REF] Stolper | Progressive visual analytics: User-driven visual exploration of in-progress analytics[END_REF], [START_REF] Fisher | Incremental, approximate database queries and uncertainty for exploratory in Large Data Analysis and Visualization (LDAV)[END_REF], [START_REF] Fisher | Exploratory visualization involving incremental, approximate database queries and uncertainty[END_REF], [START_REF] Rosenbaum | Progressive refinement: more than a means to overcome limited bandwidth[END_REF], [START_REF] Fekete | ProgressiVis: a Toolkit for Steerable Progressive Analytics and Visualization[END_REF]. Progressive Insights [START_REF] Stolper | Progressive visual analytics: User-driven visual exploration of in-progress analytics[END_REF] is an example of a progressive system that allows the user to analyze common patterns in event sequence medical data. Fisher et. al. [START_REF] Fisher | Trust me, i'm partially right: incremental visualization lets analysts explore large datasets faster[END_REF] presented sampleAction, a tool that simulates progressive queries over large datasets, and discuss and analyze the implications of such queries and different confidence metrics through quantitative case studies. While this approach has similar drawbacks to sampling (e.g., bad for finding outliers, does not work with ordered data, requires statistical sophistication from users), many authors [START_REF] Fisher | Incremental, approximate database queries and uncertainty for exploratory in Large Data Analysis and Visualization (LDAV)[END_REF], [START_REF] Fisher | Exploratory visualization involving incremental, approximate database queries and uncertainty[END_REF], [START_REF] Crotty | Vizdom: interactive analytics through pen and touch[END_REF], [START_REF] Fekete | ProgressiVis: a Toolkit for Steerable Progressive Analytics and Visualization[END_REF] advocate for its usefulness in exploratory data analysis. For example, it allows users to decide what level of accuracy they need and provide opportunities to inject user-steerability into algorithms and computations. However, the effects of this approach in terms of user performance and behavior have not yet been analyzed in detail.

Latency in Computer Systems

As many have argued [START_REF] Heer | Interactive dynamics for visual analysis[END_REF], [START_REF] Hanrahan | Analytic database technologies for a new kind of user: the data enthusiast[END_REF], the goal of visual data exploration systems is to operate at a rate that matches the pace of data analysts. Systems should therefore attempt to keep query response times below thresholds that will make users lose focus on their current task. A study by Liu et. al. [START_REF] Liu | The effects of interactive latency on exploratory visual analysis[END_REF] shows that latencies of 500ms have significant effects on user performance in data exploration scenarios. In other domains, including web search [START_REF] Brutlag | Speed matters for google web search[END_REF] and video games [START_REF] Beigbeder | The effects of loss and latency on user performance in unreal tournament 2003 R[END_REF], even smaller latencies (300ms and 100ms, respectively) can negatively influence user performance. Frameworks proposed by Nielsen and others [START_REF] Seow | Designing and engineering time: The psychology of time perception in software[END_REF], [START_REF] Nielsen | Powers of 10: Time scales in user experience[END_REF], [START_REF] Card | The information visualizer, an information workspace[END_REF], [START_REF] Shneiderman | Response time and display rate in human performance with computers[END_REF] suggest that responses within one second allow users to stay on their current train of thought, while response times over ten seconds exceed the average attention span. We use these models to justify the different delay times used in our experiment.

EXPERIMENTAL DESIGN

The aim of this work is to investigate how progressive visualizations influence users during exploratory data analysis, as well as how these techniques compare to blocking and instantaneous visualizations. We use a full-factorial 3 visualization conditions (blocking, instantaneous, progressive) × 2 dataset-delay conditions (6s, 12s) × 2 dataset-order conditions (123, 312) experiment. Our experiment is based on work by Liu et. al. [START_REF] Liu | The effects of interactive latency on exploratory visual analysis[END_REF] and Guo et. al. [START_REF] Guo | A case study using visualization interaction logs and insight[END_REF], which suggest using a hybrid evaluation approach that uses system logs and insight-based metrics coded from think-aloud protocols in order to analyze both (1) user interactions and (2) analysis performance. We expect that users will generate more insights per minute with instantaneous visualizations than with progressive ones (H1) and that users will generate more insights per minute with progressive visualizations than with blocking ones (H2). Furthermore, we anticipate user activity levels to be higher with instantaneous visualizations than with progressive ones (H3) as well as higher with progressive visualizations than with blocking ones (H4). This section provides a detailed description of the experimental design.

Visualization Conditions

In order to understand how progressive visualizations influence users in terms of knowledge discovery and interaction behavior, we compare them against two baseline visualization conditions: (1) blocking and (2) instantaneous. The instantaneous condition represents a hypothetical ideal scenario where the system always return query results within the time constraint regardless of dataset size. The blocking condition represents the other extreme, which is how many current visual data exploration systems operate. For blocking visualizations, query results are displayed only after the computation over the entire dataset concludes, with visualization latencies scaling with dataset size. The progressive condition bridges the instantaneous and blocking conditions by displaying an approximate result as soon as possible and then incrementally refining the results over time. Eventually, the full and completely accurate result is displayed, as in the blocking approach, but initial approximate results are still returned as soon as possible.

Figure 2 shows a schematic illustration of the differences between these three conditions. Assuming that a given hardware platform requires x seconds to compute and display the full query result, then a system operating in the hypothetical instantaneous mode will return results immediately; a blocking system will use the full x seconds before displaying any results; and the progressive system will show inaccurate results that are incrementally refined throughout the x seconds.

Datasets

We use three datasets from different domains in our experiment. The first dataset (DS1) contains information about cars [START_REF] Lichman | UCI machine learning repository[END_REF] (8 attributes: 7 quantitative, 1 nominal), such as "acceleration" and "horsepower." The second one (DS2) includes data about wines (7 attributes: 5 quantitative, 2 nominal), with attributes including "type", "country of origin" and several different ratings. Finally, the third one (DS3) is a subset of the 1994 US census [START_REF] Lichman | UCI machine learning repository[END_REF] (9 attributes: 3 quantitative, 6 nominal). We use a fourth dataset (Titanic) to introduce participants to the system.

The datasets contain 10,000 data points each. We introduce a dataset-delay factor that has two possible values (6s or 12s) and is used to artificially delay the computation of visualizations for the progressive and blocking conditions. In the blocking case, the user will have to wait either 6s or 12s before any visualization is displayed, while the computation spans the entire time period (6s or 12s) with 10 incremental updates in the progressive case. While 6s and 12s still represent relatively small times (it is common to have datasets large enough that computations can take hours), we chose these values for two reasons: (1) these delays are above and below the 10 second threshold that is often stated as the average attention span [START_REF] Nielsen | Powers of 10: Time scales in user experience[END_REF]; and (2) they are small enough to make an in-lab user study feasible.

System

We created an experimental system specifically for our study. The UI, shown in Figure 3, consists of a list of the current dataset's attributes (a) and four visualization panels (v1-v4). Users can drag and drop attributes onto the axes of the visualization panels, and tapping on an axis cycles through different aggregation functions (e.g., count, average). Our system supports two visualization types: (1) bar charts and (2) 2D-histograms. Visualizations are always binned and never show individual data points, which allows us to fix the time required to display a visualization, even for arbitrarily large datasets. That is, the number of visual elements to render is decoupled from the size of the underlying dataset.

Selecting a bin in a visualization triggers a brushing operation in which the selected data points are highlighted in all other visualizations. In the example (Figure 3), the user selected the rightmost bar in v2. All other visualizations now shade bins in two different colors: blue indicating the overall amount and purple the amount of data points that corresponds to the user's selection. Bar height and rectangle area are scaled to reflect the number of data points that match the selection. A textbox (c) permits finer-grained brushing control through arbitrary Boolean statements (e.g., highlight all wines where price < 35 and vintage > 2010). Similarly, a second textbox (b) supports filtering operations through Boolean statements to select a specific subset of the data. Brushing and filtering operations require all affected visualizations to recompute their results from scratch, regardless of the current visualization conditions. For example, selecting an additional bar in v2, and thereby changing the currently applied brush, will force all other visualizations (i.e., v1, v3, and v4) to recompute. Depending on the visualization condition, v1, v3, and v4 will either show the results of the new brushing action instantaneously, a loading animation until the full result is available, or incrementally updated progressive visualizations. The system does not perform any form of result caching.

Upon startup, the system loads the selected dataset into memory. The system randomly shuffles the data points in order to avoid artifacts caused by the natural oder of the data and to improve convergence of progressive computations [START_REF] Dobra | Turbo-charging estimate convergence in dbo[END_REF]. We approximate the instantaneous visualization condition by computing queries over entire dataset as quickly as possible. Initial microbenchmarks for a variety of queries over the small datasets used in the experiment yielded the following measurements: time to compute a result ≈ 100ms and time to render a visualization ≈ 30ms.

Although not truly instantaneous, a total delay of only ≈ 130ms is well below the guideline of one second, therefore allowing us to simulate the instantaneous condition. We simulate the blocking condition by artificially delaying the rendering of a visualization by the number of seconds specified through the dataset-delay factor. In other words, we synthetically prolong the time necessary to compute a result in order to simulate larger datasets. While a blocking computation is ongoing, we display a simple loading animation, but users can still interact with other visualization panels, change the ongoing computation (e.g, selecting a different subset), or replace the ongoing computation with a new query (e.g., changing an axis). Figure 4 (top) shows an example of a blocking visualization over time.

We implemented the progressive condition by processing the data in chunks of 1,000 data points at a time, with approximate results displayed after each chunk. In total, we refresh the progressive visualization 10 times, independent of the dataset-delay factor. We display the first visualization as quickly as possible, with subsequent updates appropriately delayed so that the final accurate visualization is displayed after the specified dataset-delay Progressive visualizations are augmented with error metrics indicating that the current view is only an approximation of the final result. Even though 95% confidence intervals based on the standard error have been shown to be problematic to comprehend in certain cases [START_REF] Correll | Error bars considered harmful: Exploring alternate encodings for mean and error[END_REF], we still opted to use them for bar charts due to their wide usage and familiarity. We render labels with margins of error (e.g., "±3%") in each bin of a 2D-histogram. Figure 4 (bottom) shows an example of a progressive visualization and how it changes over time.

Note that the confidence intervals in the example are rather small, but their size can change significantly based on the dataset and query.

Our system is implemented in C# / Direct2D. We tested our system and ran all sessions on a quad-core 3.60GHz, 16GB RAM, Microsoft Windows 10 desktop machine with a 16:9 format, 1920x1080 pixel display.

Procedure

We recruited 24 participants from a research university in the US. All participants were students (22 undergraduate, 2 graduate), all of whom had some experience with data exploration or analysis tools (e.g., Excel, R, Pandas). 21 of the participants were currently enrolled in and halfway through an introductory data science course. Our experiment included visualization condition as a within-subject factor and dataset-delay and dataset-order as between-subject factors. Note that the dataset-delay has no direct influence on the instantaneous condition, but it is used as a between-subject factor to test if it affects the other visualization conditions. To control against ordering and learning effects, we fully randomized the sequence in which we presented the different visualization conditions to the user and counterbalanced across dataset-delay conditions. Instead of fully randomizing the ordering of datasets, we opted to create two predefined dataset-orderings and factored them into our analysis. The two possible dataset-ordering values were 123 and 312 (i.e., DS1 followed by DS2 followed by DS3 and DS3 followed by DS1 followed by DS2, respectively), and we again counterbalanced across dataset-ordering. We seed our system's random number generator differently for each session to account for possible effects in the progressive condition caused by the sequence in which data points are processed. While each participant did not experience all possible combinations of visualization, dataset-ordering, and dataset-delay conditions, all users saw all visualization conditions with one specific dataset-delay and dataset-ordering.

For example, participant P14 was assigned dataset-ordering 312, dataset-delay 6s, and visualization condition ordering [blocking, instantaneous, progressive]. In total, this adds up to four trial sub-groups, each with six participants: (datasetdelay = 6s & dataset-order = 123), (dataset-delay = 12s & dataset-order = 123), (dataset-delay = 6s & dataset-order = 312) and (dataset-delay = 12s & dataset-order = 312). Within each subgroup, we fully randomized the order in which we presented the different visualization conditions.

After a 15 minute tutorial of the system, including a summary of the visualization conditions and how to interpret confidence intervals and margins of error, we instructed the participants to perform three exploration sessions. In the case of participant P14, these three sessions included (1) blocking visualizations on DS3, (2) instantaneous visualizations on DS1, and (3) progressive visualizations on DS2 all with 6s delay. Each session was open-ended and we asked the participants to explore the dataset at their own liking and pace. We allotted 12 minutes per session, but participants were free to stop earlier We used a think-aloud protocol [START_REF] Ericsson | Protocol analysis[END_REF] where participants were instructed to report anything they found interesting while exploring the dataset. Throughout each session, we captured both screen and audio recordings and logged low-level interactions (mouse events) as well as higher-level events ("axis changed", "aggregation changed", "textbox brush / filter", "visualization brush," and "visualization updated / stopped / completed"). An experimenter was present throughout the session, and participants were free to ask any technical questions or questions about the meaning of dataset attributes. At the end of the three sessions, we asked participants to give feedback about the tool and whether they had any thoughts regarding the different visualization conditions.

Statistical Analysis

Our study is designed as a full-factorial 3 (visualization conditions) × 2 (dataset-delay conditions) × 2 (dataset-order conditions) experiment. We applied mixed design analysis of variance tests (ANOVA) with visualization condition as the within-subject factor and dataset-delay and datasetorder as the between-subject factors to assess the effects of our factors on the various metrics we computed. Note that the dataset-delay factor should have no influence on trials where the visualization condition is set to instantaneous.

We tested the assumption of sphericity using Mauchly's test, and we report results with corrected degrees of freedom using Greenhouse-Geisser estimates if violated. We report all significant (i.e., p < 0.05) main and interaction effects of these tests. For significant main effects, we conducted further analysis through Bonferroni corrected post hoc tests and for more nuanced interpretation, we opted to include Bayes factors for certain results and report BF 10 factors along with corresponding significance labels [START_REF] Lee | Bayesian cognitive modeling: A practical course[END_REF]. Effect sizes are reported through Pearson's r coefficient, and significance levels are encoded in plots using the following notation: * p < 0.05; * * p < 0.01; * * * p < 0.001.

ANALYSIS OF VERBAL DATA

Inspired by previous insight-based studies [START_REF] Liu | The effects of interactive latency on exploratory visual analysis[END_REF], [START_REF] Guo | A case study using visualization interaction logs and insight[END_REF], [START_REF] Gomez | An insight-and task-based methodology for evaluating spatiotemporal visual analytics[END_REF], we manually coded insights from the audio recordings and screen captures. An insight is a nugget of knowledge extracted from the data, such as "France produces more wines than the US". We followed the verbal data segmentation approach proposed by Liu et. al. [START_REF] Liu | The effects of interactive latency on exploratory visual analysis[END_REF] and adopted their coding process in which the first author did the bulk of the coding, but we iteratively revised finished codes with collaborators to reduce bias. In our case, we decided to count observations that are within the same visualization, have the same semantics, and are on the same level of granularity as one insight. For example: "It looks like country 1 makes the most cars, followed by country 2 and country 3" was coded as one single insight, whereas an observation across two visualizations (through brushing) such as "Country 1 makes the most cars and seems to have the heaviest cars" was counted as two separate insights. We did not categorize insights or assign any quality scores or weights to insights nor did we quantify accuracy or validity of insights. For our analysis, all insights were treated equally.

Number of Insights per Minute

In order to get a quantifiable metric for knowledge discovery, we normalized the total insight count for each session by the duration of the session. Figure 5 shows this resulting "number of insights per minute" metric across different factors. Our results show that this metric was significantly affected by the type of visualization condition (F (1.452, 29.039) = 6.701, p < 0.01).

Post hoc tests revealed that the instantaneous condition showed a slight increase of number of insights per minute over the progressive condition (1.477 ± 0.568 vs. 1.361 ± 0.492, respectively), which was not statistically significant (p = 1.0, r = 0.080). However, the blocking condition reduces the number of insights per minute to 1.069 ± 0.408, which differed significantly from both the progressive (p < 0.05, r = 0.292) and instantaneous (p < 0.001, r = 0.383) conditions. A Bayesian Paired Samples T-Test that tested if measure 1 < measure 2 revealed strong evidence for an increase in insights per minute from the blocking condition to the progressive condition (BF 10 = 13.816), extreme evidence for an increase from blocking to instantaneous (BF 10 = 134.561), and moderate evidence for no change or a decrease between instantaneous and progressive (BF 10 = 0.130). In summary, blocking visualizations produced the fewest number of insights per minute, whereas the instantaneous and progressive visualizations performed equally well.

Insight Originality

Similar to [START_REF] Guo | A case study using visualization interaction logs and insight[END_REF], we grouped insights that encode the same nugget of knowledge together and computed the originality of an insight as the inverse of the number of times that insight was reported by any participant. That is, insights reported more frequently by participants received lower overall originality scores. A participant received an insight originality score of 1 if he or she had only unique insights (i.e., insights found by no other users). The lower the score, the less original the insights were on average. We averaged originality scores across insights for each session and show plots for this insight originality metric across different factors in Figure 6. We did not find any significant effects that influence the insight originality metric, which indicates that the originality score seems unaffected by any of the factors for which we controlled.

ANALYSIS OF INTERACTION LOGS

To analyze how our visualization conditions affect user behavior, we computed several metrics from either the mousemovement events or the high-level system-specific events.

Visualization Coverage

We were interested in analyzing how much of the possible space of visualizations our participants covered. To create a metric for visualization coverage, we computed the set of unique visualizations possible within each dataset. We considered all attributes, both visualization types supported by our system, and all possible aggregation functions. However, we ignored the axis-to-attribute mapping (e.g., a visualization with attribute A on the x-axis and attribute B on the y-axis is considered the same as if the axes were flipped). We then extracted and counted up all visualizations a participant created during a session from the interaction logs. Our final visualization coverage metric is the percentage of possible visualizations a participant created per minute.

Figure 7 plots this metric across different factors. Our analysis shows that the type of visualization condition significantly affects the percentage of the total visualizations that a participant covered per minute (F (2, 40) = 9.847, p < 0.001). Post hoc tests revealed that the instantaneous condition showed a slight increase in percentage over the progressive condition (1.553% ± 0.690% vs. 1.311% ± 0.651%, respectively), which was not statistically significant (p = 0.226, r = 0.251). However, the blocking condition reduces the percentage of total visualizations covered per minute to 0.845% ± 0.449%, which differed significantly from both the progressive (p < 0.01, r = 0.545) and instantaneous (p < 0.001, r = 0.740) conditions. Participants explored more visualizations per minute with the instantaneous and the progressive condition and there is no significant difference between the two.

Note that we did not assign any "importance" or "quality" scores to different visualizations. All visualizations have the same weight, even though some visualizations might be more informative than others or multiple visualizations might convey similar insights. Similarly, our visualization coverage metric is not designed to compare across different users or sessions, such that two users could have the exact same coverage score while looking at completely different parts of the dataset.

Number of Brush Interactions per Minute

While the coverage metric considers the number of possible static visualizations, it does not consider brushing. We measured the brushing interactions by counting the number of brush events from the interaction logs, which we then normalized by the duration of a session. The results, shown in Figure 8, demonstrate that the type of visualization condition significantly affected the number of brush interactions per minute (F (1.335, 26.690) = 17.620, p < 0.0001). Post hoc tests revealed that the instantaneous condition showed a significant increase in number of brush interactions per minute over the progressive condition (4.640 ± 4.095 vs. 2.108 ± 1.744, p < 0.01, r = 0.316), as well as over the blocking condition (1.190 ± 0.737 , p < 0.001, r = 0.413). Furthermore, brushing interactions per minute for the progressive condition were significantly different than for the blocking condition (p < 0.05, r = 0.29).

Additionally, our results showed a significant betweensubject effect for dataset-order (F (1, 20) = 4.568, p < 0.05)). A post hoc test revealed that dataset-order 312 increased the number of brush interactions per minute over datasetorder 123 significantly (3.365 ± 3.284 vs. 1.927 ± 2.429, p < 0.05, r = 0.242). We hypothesize that this effect is due to the structure of DS3 (census dataset), which in turn leads to a learning effect. DS3, to which users where exposed first in dataset-order 321, has considerably fewer qualitative attributes than the other datasets. Users could not use strategies such as looking for trends in 2D histograms or compute averages across attributes and reverted to using the brushing functionality to correlate across different populations of the data. We often observed users manually cycle through one attribute and look for changes in another attribute. For example, users would create two histograms (e.g., one for "marital status" and one for "education") and then manually select different values in the first histogram (e.g., "married", "widowed") to determine whether the second histogram for that subpopulation differed from the overall population.

Visualizations Completed

We extracted the number of times a visualization was completed-the participant waited the full amount of time specified by the dataset delay until a visualization was completely accurate-from our log data. We then divided this count by the number of interactions that forced a visualization to recompute (e.g., axis change, brushing interaction). This gives us the percentage of times an interaction led to a fully accurate visualization. Figure 9 visualizes this metric for different factors. Note that, by definition, this metric is always 100% for the instantaneous condition and we therefore exclude it from post hoc tests. Our results show that the percentage of completed visualizations was significantly affected by the type of visualization condition (F (2, 40) = 169.972, p < 0.0001). Post hoc tests revealed that the blocking condition showed an increase of percentage of completed visualizations over the progressive condition (56.22% ± 15.01% vs. 45.99% ± 14.14%, respectively), which was statistically significant (p < 0.05, r = 0.296). In short, people often moved ahead without waiting for the full result.

Mouse Movement per Minute

Finally, to compute a simple metric of a participant's activity level, we calculated the distance in pixels the mouse was moved per minute and show the results in Figure 10. Our analysis shows that the type of visualization condition significantly affected the mouse movement per minute metric (F (2, 40) = 5.431, p < 0.01). Post hoc tests revealed that the instantaneous condition showed a slight decrease of mouse movement per minute over the progressive condition (303.799 ± 128.243 vs. 313.912 ± 139.1299, respectively), which was not statistically significant (p = 1.0, r = 0.051). However, the blocking condition reduces movement per minute to 258.060 ± 131.536, which differed significantly from both the progressive (p < 0.05, r = 0.315) and instantaneous (p < 0.05, r = 0.245) conditions. Mouse movement is a crude indication of a user's activity level, and our test shows that these levels are lowest with the blocking conditions. Again we find no difference between instantaneous and progressive visualizations.

PERCEPTION OF VISUALIZATION CONDITIONS

During our exit interview, we asked participants to provide feedback about the different visualization conditions they experienced. Most participants liked the instantaneous visualizations best, but preferred the progressive ones over blocking. Below are quotes from our participants that describe these preferences:

"I liked the progressive one better than blocking, but obviously the instantaneous one is best." "I initially felt the loading one A few participants expressed positive remarks towards blocking visualizations or commented that they adapted their strategies because of the wait-time:

"The slow one [blocking] made me feel more confident about what I saw, because there is lots of data behind it." "It [blocking] actually helped me to use the time to think. But it also might limit you from finding really interesting facts, because you're going in with an idea, you're using the loading time to come up with things that you think might be true. Without the loading time you could just randomly mess around with the data and find interesting things." "I used the loading time to do something else."

DISCUSSION

Our analysis provides evidence that overlaps with findings in previous work, such as by Liu et. al. [START_REF] Liu | The effects of interactive latency on exploratory visual analysis[END_REF], which finds that even small latencies have an impact on user performance in exploratory settings. Our results show that knowledge discovery and user activity measures are negatively influenced by the blocking condition when compared to instantaneous visualizations. This intuitively makes sense and it is widely acknowledged that low latency leads to improved user experience and user performance. The more interesting evidence we present in this paper arises when we compare blocking to progressive visualizations or instantaneous to progressive visualizations.

The difference between the progressive and blocking conditions is that users can see approximate results while a query is ongoing rather than a loading animation. However, the overall delay until a final, 100% accurate visualization arrives is exactly the same for both conditions. Yet, our data shows that users generated more insights per minute, had a higher visualization coverage percentage, and displayed higher levels of mouse movements when given progressive rather than blocking visualizations. Additionally, the percentage of completed visualizations (i.e., visualizations where users waited the full number of dataset-delay seconds to get the final answer) is higher in the blocking than the progressive condition. This result suggests that users are efficiently using these in-between and approximate visualizations to either pre-process information, extract insights early, or to decide that the result is not what they were looking for and then move on to the next visualization. A participant expressed it this way: "It's much easier to look at a rough picture of the final data rather than just to see nothing at all. You can start to get an idea of relationships and things like that.". Based on our analysis we accept hypothesis H2 (more insights with progressive than blocking) as well as H4 (higher user activity levels with progressive than blocking).

The instantaneous visualization condition represents a hypothetical system that can provide results to the user almost immediately for any query, minimizing interference to the user's thought process. While increasing amounts of data make systems that provide instantaneous visualizations practically infeasible, the evidence presented in our analysis shows that progressive visualizations might perform almost as well. We did not observe any significant differences across all of our metrics for the instantaneous and the progressive conditions except for the brush interactions per minute metric. We can therefore not accept H1 (more insights with instantaneous than progressive) or H3 (higher user activity levels with instantaneous than progressive).

However, there are some open questions that our study does not address. While our data suggests that approximate visualizations might help users grasp certain characteristics of a dataset rapidly, they simultaneously introduce a set of new challenges. For example, prior work in psychology [START_REF] Cumming | Inference by eye: reading the overlap of independent confidence intervals[END_REF], [START_REF]Understanding the new statistics: Effect sizes, confidence intervals, and meta-analysis[END_REF] has shown that even experts sometimes have trouble interpreting standard accuracy metrics (e.g., confidence intervals, error bars). While we did not observe any cases in our study where participants misinterpreted visualizations based on their uncertainty, these cases are also hard to capture. Our participants often reported high-level trends like "these two categories seem about equal in counts", "there is a slight negative correlation between these two variables". How much of a change would need to occur between the approximate and the final visualization before users would retract these insights? Progressive visualizations expose fundamental limitations when it comes to exploration sessions on datasets where it is important to capture insights that are based on outliers or small subsets of the data. It might take a long time to sample those rare datapoints, and visualizations can therefore take a long time to converge to a view that shows those features prominently.

The second open question is how refresh rates for progressive visualizations affect user interaction and whether or not, in the extreme case, a simple one-sample visualization provides the same benefits. For example, techniques such as the one presented by Stolper et. al. [START_REF] Stolper | Progressive visual analytics: User-driven visual exploration of in-progress analytics[END_REF] that allow users to decide when to update a visualization might be less distracting than our approach of constant regular updates. There are reasons to believe that progressive visualizations have advantages over approximate visualization that are based on a single sample. For example, users can decide individually when a visualization is accurate enough for their liking or their current task. We observed that some users tend to hold off with reporting insights until after several progressions, especially for visualizations where there was a high visual variance between updates. In other cases, participants waited for uncertainty metrics to be in a range small enough for them to make judgments. One of our participants stated: "I didn't want to draw any conclusion right away. Especially in this one visualization where the error bar was across the whole graph, I decided to wait.". While these anecdotal findings make a case for progressive visualizations over simple sampling-based ones, our study does not provide the means for a quantitative comparison between the two.

FUTURE WORK & CONCLUSION

We investigated how progressive visualizations affect users in exploratory data analysis scenarios. Through a controlled experiment, we compared progressive visualizations to blocking and instantaneous visualizations and found significant differences in insight-based metrics and user activity levels across all three approaches. We observed that progressive visualizations outperform blocking visualizations in almost all metrics and that progressive visualizations do not significantly differ from the ideal scenario of instantaneous visualizations in terms of generated insights per minute, insight originality, or visualization coverage percentage. However, this study is just a first step towards understanding progressive visualizations. We studied only simple forms of uncertainty visualizations and we also excluded other factors from the study, such as if users are able to fully grasp the meaning of approximate answers and how distracting the visualization updates are. We plan to carry out several follow-up studies where we compare different update strategies and different types of uncertainty representations. Furthermore, we intend to get a better understanding of how a user's behavior changes between instantaneous and progressive visualizations through indepth sequential interaction log analysis, potentially coupled with eye-tracking data. Gaining such understanding would help optimize visual representations and interactions with progressive visualizations.

Fig. 1 .

 1 Fig. 1. Two coordinated visualizations. Selections in the left filter the data shown in the right.

Fig. 2 .

 2 Fig. 2. Schematic time scale showing when different visualization condition display results.

Fig. 3 .

 3 Fig. 3. Screenshot of our experimental system.

Fig. 4 .

 4 Fig. 4. The blocking and progressive visualization conditions.

  [progressive] was weird because graphs change over time. But after seeing this one [blocking] I can appreciate the value of it [progressive]." "[Blocking] slows the process. But when it's one or the others [instantaneous or progressive], I can see one thing and then that naturally leads to the next thing I want to do. Here [blocking] I have to keep track of the last thing that I loaded while I do something else. It's [blocking] just more stilted and less continuous." "You see a rough picture in the beginning and then you can think about it while it's actually finishing. That's way better than just a loading animation."

Fig. 5 .

 5 Fig. 5. Insights per Minute: Boxplot (showing median and whiskers at 1.5 interquartile range) and overlaid swarmplot for insights per minute (left) overall, (middle) by dataset-delay, and (right) by dataset-order. Higher values indicate better.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Insight Originality: Boxplot (showing median and whiskers at 1.5 interquartile range) and overlaid swarmplot for insight originality (left) overall, (middle) by dataset-delay, and (right) by dataset-order. Higher values indicate more original insights.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Visualizations Completed Percentage: Boxplot (showing median and whiskers at 1.5 interquartile range) and overlaid swarmplot for for completed visualization % (left) overall, (middle) by dataset-delay, and (right) by dataset-order.

ACKNOWLEDGMENTS

The authors wish to thank Angelina Maric from the University of Z ürich for her help with statistical analysis as well as Professor Andries van Dam, Robert Zeleznik, and the reviewers for their comments, feedback, and suggestions.

Emanuel


Giorgio Giorgiani
email: giorgio.giorgiani@univ-amu.fr

Hervé Guillard
email: herve.guillard@inria.fr

Boniface Nkonga
email: boniface.nkonga@inria.fr

SHOCK CAPTURING COMPUTATIONS WITH STABILIZED POWELL-SABIN ELEMENTS

Keywords: Shock capturing, high-order, C1 approximations, splines, CFD, compressible flow

In the last recent years, thanks to the increasing power of the computational machines, the interest in more and more accurate numerical schemes is growing. Methods based on high-order approximations are nowadays the common trend in the computational research community, in particular for CFD applications.

This work is focused on Powell-Sabin (PS) finite elements, a finite element method (FEM) based on PS splines. PS splines are piecewise quadratic polynomials with a global C1 continuity, defined on conforming triangulations. Despite its attractive characteristics, so far this scheme hasn't had the attention it deserves. PS splines are adapted to unstructured meshes and, contrary to classical tensor product B-splines, they are particularly suited for local refinement, a crucial aspect in the analysis of highly convective and anisotropic equations. The additional global smoothness of the C1 space has a beneficial stabilization effect in the treatment of advection-dominated equations and leads to a better capturing of thin layers. Finally, unlike most of other typology of high-order finite elements, the numerical unknowns in PS elements are located in the vertices of the triangulation, leading to an easy treatment of the parallel aspects.

Some geometrical issues related to PS elements are discussed here, in particular, the generation of the control triangles and the imposition of the boundary conditions. The PS FEM method is used to solve the compressible Euler equation in supersonic regime. A classical shock-capturing technique is used to reduce the oscillation around the discontinuity, while a variational multiscale formulation is used to introduce numerical diffusion in the streamwise direction. Some typical numerical examples are used to evaluate the performance of the PS discretization.

INTRODUCTION

Numerical methods with high degrees of regularity are of interest in many physical applications, such as magnetohydrodynamics, shells analysis and vibrations. In particular, finiteelement methods based on spline shape functions have raised a particular interest in the recent years. On one hand, C 1 spline methods allow to discretize higher-order derivatives (that is, derivatives of order superior to the second), and on the other hand, they provide accurate geometrical representation of the computational domain, see for example [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Bazilevs | Isogeometric analysis using T-splines[END_REF] in the context of Isogeometric Analysis. Moreover, it has also been shown recently that the additional global smoothness of the spline interpolant introduces stability to the numerical solution for highly convective equations, see [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF] and [START_REF] Speleers | Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems[END_REF], and also for turbulence computations, see [START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF]. Hence, finite element methods based on spline shape functions can be advantageous also in the context of fluid-dynamic problems.

Nowadays, splines are extensively used in the graphical design industry to create smooth surfaces. The success of splines is basically due to the fact that they have very attractive and unique characteristics: they have a compact representation, they are able to represent curves and surfaces with an arbitrary level of regularity, simply increasing the polynomial degree of the basis, and they are efficient in terms of locality, that is, they can be modified locally by moving a control point, without perturbing the rest of the curve. Spline curves are usually represented as a linear combination of basis functions, called B-splines. The basis functions of a given order p are defined with a recursive relation starting with piecewise constant basis functions (that is, p = 0),

B i,0 (ξ) = 1 if ξ i ≤ ξ < ξ i+1 , 0 otherwise B i,p (ξ) = ξ -ξ i ξ i+p -ξ i B i,p-1 (ξ) + ξ i+p+1 -ξ ξ i+p+1 -ξ i+1 B i+1,p-1 (ξ),
where the points ξ i are the knots of the spline, i = 1, 2, . . . , n + p + 1, and n is the number of basis functions. Hence, the spline curve expressed as a combination of B-splines is

S(ξ) = n i=1 B i,p P i , (1)
where P i are the control points of the spline.

Starting from the one-dimensional case (1), B-splines can be extended to higher dimensions through a tensor product representation. However, tensor product B-splines are restricted to structured rectangular meshes. The refinement procedure with tensor B-splines relies on the insertion of knots, leading hence to a global modification of the domain discretization. Thus, no local refinement is possible. Various solutions are proposed in literature to solve this limitation, for example, the introduction of T-splines in [START_REF] Bazilevs | Isogeometric analysis using T-splines[END_REF].

This work is focused on another approach, based on the definition of bivariate splines on irregular triangulations, that is, Powell-Sabin (PS) splines [START_REF] Powell | Piecewise quadratic approximations on triangles[END_REF][START_REF] Jaxon | Isogeometric analysis on triangulations[END_REF][START_REF] Speleers | Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems[END_REF][START_REF] Dierckx | On calculating normalized Powell-Sabin B-splines[END_REF]. PS splines are piecewise quadratic polynomial with C1 continuity, defined on a unstructured triangulation of the domain. The major advantage with respect to tensor product B-splines is that PS splines allow a straightforward adaptive refinement of the mesh, which is a key ingredient in the simulation of anisotropic equations. Another advantage of PS elements is related to the fact that the unknowns are concentrated on the nodes of the triangulation. This means that there are no unknowns on the element faces or in the interior of the elements. This is not the case for other C 1 interpolation techniques, such as, for example, the Clough-Tocher elements. The consequence is that the stencil is the same for each node of the mesh, which produces a linear system matrix with a regular block shape that facilitates the implementation aspect, enhances the parallel performance and improves the efficiency of the linear system solution.

In this work, a PS finite element method is presented for the solution of the two dimensional (2D) compressible Euler equation in supersonic regime. A SUPG stabilization technique is used to introduce numerical dissipation in the streamwise direction only, and an isotropic diffusion terms is added to avoid oscillation around the discontinuities.

The outline of the paper is as follows. Section 2 introduces the PS splines and their representation. Section 3 deals with the geometrical and mathematical tools to define the PS elements and in particular the generation of the shape functions. A detailed discussion is also given on the imposition of the boundary conditions. Finally, in Section 4 two classical numerical examples are presented.

POWELL-SABIN SPLINES

In this section is introduced the framework and the notation to define spline functions on a triangulation of a polygonal domain. Starting from the definition of bivariate polynomials on a triangle, the goal is to obtain an interpolant of a generic function with C 1 continuity. Considering a single triangle Ω k of vertices V i = (x i , y i ), with i = 1, 2, 3, any bivariate polynomial p(x, y) of degree ≤ 2 in the space Π 2 = { 2 i=0 2-i j=0 a i,j x i y j , a i,j ∈ R} can be written in the Bernstein-Bézier representation

p(x, y) := p(ξ) = i+j+k=2 i,j,k≥0 b i,j,k B 2 i,j,k (ξ(x, y)), (2)
where ξ = (ξ 1 , ξ 2 , ξ 3 ) are the barycentric coordinates of a point (x, y) ∈ R 2 and B 2 i,j,k (ξ) are the Bernstein polynomials of degree 2 in Ω k , that is

B 2 i,j,k = 2! i!j!k! ξ i 1 ξ j 2 ξ k 3 . (3)
The set of Bernstein polynomials B 2 i,j,k (ξ) is a basis for the space of polynomials Π 2 , see [START_REF] Lai | Spline Functions on triangulations[END_REF], hence the six coefficients b i,j,k uniquely define the second order polynomial p(x, y) and they are called the Bézier ordinates with respect to the triangle Ω k 1 . Let Ω ⊂ R 2 be a polygonal domain with boundary ∂Ω. Let T be a conforming triangulation of Ω with vertices V l , l = 1, . . . , N v and elements Ω e , e = 1, . . . , N e . Defining S 2 1 as the linear space of piecewise quadratic polynomials on T , the following interpolation problem is considered: given any set of triples (f l , f xl , f yl ), l = 1, . . . , N v, find s(x, y) ∈ S 2 1 such that,

s(V l ) = f l , ∂s ∂x (V l ) = f xl , ∂s ∂y (V l ) = f yl , l = 1, . . . , N v. (4)
It is clear that such a problem has no solution in general: in fact, problem (4) requires the imposition of nine parameters to define the second order polynomial on each triangle, while only six coefficients are available, see equation [START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF].

Thus, to obtain a solution to the interpolation problem (4), one possibility is to subdivide each triangle in sub-triangles and to define the interpolant in the refined triangulation. One of the solutions proposed by Powell and Sabin in [START_REF] Powell | Piecewise quadratic approximations on triangles[END_REF] is based on the subdivision of each triangle in T into six smaller triangles (PS6-split). Hence, the conditions in (4) are imposed only on the vertices of the original triangulation, while in the other added nodes only C 1 continuity of the interpolating function is imposed. More details can be found in [START_REF] Powell | Piecewise quadratic approximations on triangles[END_REF].

The so called Powell-Sabin refinement of T is denoted with T * and is obtained dividing each triangle in T in six sub-triangles. The procedure to define T * is described in [START_REF] Speleers | Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems[END_REF] and reads as follows.

• Select a split point C k inside each triangle Ω k and connect it to the three vertices of Ω k .

• For each pair of triangles Ω p and Ω q , with a common edge, connect the two points C p and C q . If Ω p is a boundary triangle, connect also C p to an arbitrary point on the boundary edge.

The PS refinement defines also a set of points, called PS points, associated to each vertex.

The PS points associated to the vertex V l are defined as the midpoints of all the edges of the PS refinement containing V l , plus the point V l itself. The PS points are fundamental for the definition of the shape functions as explained in the next section.

Having defined the PS refinement T * , the linear space of piecewise quadratic polynomials with C 1 continuity can be denoted as S 1 2 (T * ). Each element S 1 2 (T * ) is uniquely defined by its values and derivatives at the vertices of T , thus the functional space S 1 2 (T * ) has dimension 3N v .

POWELL-SABIN FINITE ELEMENTS

In this section PS splines are introduced in a finite-element framework to solve partial differential equations, defining shape functions belonging to the space S 1 2 (T * ). Hence, it is necessary to express each element s(x, y) ∈ S 1 2 (T * ) as a linear combination of PS spline basis functions, that is

s(x, y) = Nv l=1 3 r=1 c l,r B (r) l (x, y),
where the functions B

(f i , f xi , f yi ) = 0 except for (f l , f xl , f yl ) = (α, β, γ) = 0. The quantity (α, β, γ
) is called a triple and represents the value of the function and the derivatives with respect the Cartesian axes in a given vertex. Then, it is easy to see that each basis B (r) l , r = 1, 2, 3, vanishes outside the molecule M l of vertex V l defined as the union of all triangles Ω e containing V l . This guarantees the compactness of the support for PS B-splines.

In principle, it is possible to define the PS B-splines simply choosing, for each vertex V l ∈ T , three linearly independent vectors (α l,r , β l,r , γ l,r ), r = 1, 2, 3, and solving the interpolation problem (4). This procedure generates three linearly independent basis functions for each vertex. However, in a finite element context, it is interesting to work with normalized B-splines, that is, basis functions that form a partition of unity on Ω, i.e.

B (r) l (x, y) ≥ 0 ∀x, y ∈ Ω, Nv l=1 3 r=1 B (r) l (x, y) = 1 ∀x, y ∈ Ω,
which leads to the following constrains on the triples

3 r=1 α l,r = 1, 3 r=1 β l,r = 0, 3 r=1 γ l,r = 0.
This property is also crucial to ensure the completeness of the basis.

The procedure proposed in [START_REF] Dierckx | On calculating normalized Powell-Sabin B-splines[END_REF] allows to obtain PS B-splines that form a partition of unity and with good characteristic from a numerical point of view. It is based on the fact that each set of three independent B-splines associated to a vertex V l uniquely define a control triangle

T l = (Q l,1 , Q l,2 , Q l,3
) for the PS spline, with vertices Q l,j = (X l,j , Y l,j ), j = 1, 2, 3. The control triangle of a PS spline as a role similar to the control polygon for univariate splines in terms of control of the shape of the spline. In fact, in [START_REF] Dierckx | On calculating normalized Powell-Sabin B-splines[END_REF] it is shown that the control triangle and the coefficients c l,j define three control points (X l,j , Y l,j , c l,j ) defining a spatial plane that is always tangent to the spline surface at the vertex V l . Optimal shape functions from a numerical point of view are associated to a minimal area control triangle [START_REF] Dierckx | On calculating normalized Powell-Sabin B-splines[END_REF].

The points Q l,j are related to the three triples and the vertex V l via the linear system

  α l,1 α l,2 α l,3 β l,1 β l,2 β l,3 γ l,1 γ l,2 γ l,3     X l,1 Y l,1 1 X l,2 Y l,2 1 X l,3 Y l,3 1   =   x l y l 1 1 0 0 0 1 0   . ( 5
)
The positivity of the B-splines associated to the vertex V l is guaranteed if all the barycentric coordinates with respect the control triangle T l of all the PS points associated to V l are positive. This means geometrically that all the PS points associated to V l are contained in the control triangle T l , see Figure 1. Hence, the procedure of computing the three shape functions for each vertex is:

• find a suitable control triangle for the vertex V l that contains all the PS points associated to V l ;

• compute the triple defining the three shape functions solving the linear system (5).

The control triangles for each vertex are constructed looking for the minimal area triangle that contains all the PS points. This lead to an optimization problem consisting in circumscribing a convex polygon with a minimal area triangle. An optimal algorithm is proposed by O'Rourke in [START_REF] O'rourke | An optimal algorithm for finding minimal enclosing triangles[END_REF] for solving such a problem, with a computational cost O(n), with n the number of PS points. This is the algorithm used here in the numerical tests. In Figure 1 is depicted an example of control triangle, with the relative PS control points, obtained with the O'Rourke algorithm.

In order to understand the role of the control triangle in the definition of the basis functions, a comparison of the condition number of the linear system matrix as been carried out for three structured computational meshes for a square domain [0, 1]×[0, 1]. In the first case, the minimal area triangle have been used to compute the basis. In the second case, the minimal area triangle has been expanded moving (to a distance equal to five times the original one) each vertex along the line joining the vertex with the node to which the triangle belongs. Table 1 shows the results: it is possible to notice how increasing the control triangle area degrades the condition number of mass matrix and the stiffness matrix.

STABILIZED APPROXIMATION OF THE EULER EQUATIONS WITH POWELL-SABIN ELEMENTS

The 2D inviscid Euler equations of gas dynamics are considered,

∂ρ ∂t + ∇•(ρu) = 0 ∂ρu ∂t + ∇•(ρu ⊗ u + pI) = 0 ∂ρE ∂t + ∇•((ρE + p)u) = 0 (6)
where ρ is the density, u is the velocity vector, E is the total energy, I is the identity matrix and p is the pressure, which is defined by the equation of state p = (γ g -1)ρ , γ g = 1.4 is the ratio of specific heats and = E -1/2 u 2 is the internal energy.

Problem ( 6) is recast in a conservative form, ∂U ∂t

+ ∇•F = 0, (7)
having defined the vector of conservative variables U = (ρ, ρu, ρv, ρE) T and the flux tensor F with components, in Cartesian coordinates,

F = F x F y , F x =        ρu ρu 2 + p ρuv (ρE + p)u        , F y =        ρv ρuv ρv 2 + p (ρE + p)v       
, where u and v are the x and y components of the velocity vector.

A stabilized formulation based on PS splines is considered to solve problem [START_REF] Speleers | Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems[END_REF]. To simplify the presentation, in an abuse of notation, the same symbol is used for the numerical approximation, belonging to the finite dimensional space, and the exact solution U . Hence, the stabilized finite element problem corresponding to [START_REF] Speleers | Isogeometric analysis with Powell-Sabin splines for advection-diffusion-reaction problems[END_REF]

reads: find U ∈ [S 1 2 (T * )] 4 such that Ω V • ∂U ∂t dΩ - Ω ∇V . . F dΩ + ∂Ω F V • n dΓ + D SU P G + D SC = 0, (8)
for all V ∈ [S 1 2 (T * )] 4 . The terms D SU P G and D SC are the stabilization terms and they will be discussed later.

A fully implicit first-order time integration scheme is employed to approximate the time derivatives. The non-linear convective term is linearized using a Newton-Raphson procedure of the type

Ω ∇V . . F k dΩ ≈ Ω ∇V . . A k-1 U k dΩ,
where k stands for the current N-R iteration and the third order tensor A is the Jacobian of the flux,

A = ∂F ∂U ,
and can be expressed as a second order tensor for each Cartesian direction, that is

∂F x ∂U =     0 1 0 0 γg-1 2 u 2 -u 2 (3 -γ g )u (1 -γ g )v γ g -1 -uv v u 0 γg-1 2 u u 2 -uH -u 2 (γ g -1) + H (1 -γ g )uv γ g u     , ∂F y ∂U =     0 0 1 0 -uv v u 0 γg-1 2 u 2 -v 2 (1 -γ g )u (3 -γ g )v γ g -1 γg-1 2 v u 2 -vH (1 -γ g )uv -v 2 (γ g -1) + H γ g v     ,
having defined the entalpy, H = E + p/ρ.

Imposition of the boundary conditions

Solid walls and free-stream boundaries are considered for defining the boundary conditions of system [START_REF] Jaxon | Isogeometric analysis on triangulations[END_REF]. Walls are considered as slip boundary conditions, that is, the normal component of the flux is set to zero on a wall. This kind of condition does not entail particular problems, and it is weakly imposed by the boundary integral in [START_REF] Dierckx | On calculating normalized Powell-Sabin B-splines[END_REF]. On the contrary, inlet and outlet boundary conditions are imposed setting different components of the unknown vector at the boundary nodes, depending on the subsonic/supersonic regime at the boundary, see [START_REF] Donea | Finite element methods for flow problems[END_REF]. The problem of imposing boundary values deserves particular attention in the framework of PS finite elements. Since these kind of conditions are typically set on straight boundaries, the interest is focused on imposing boundary values with PS elements on polygonal boundaries.

A careful choice of the control triangle associated to the boundary node V l simplify the task of imposing boundary values. Two different cases must be considered: the boundary node belongs to two edges forming an angle different to π, or equal to π. From equation ( 5), the following relations are derived

α l,1 Q l,1 + α l,2 Q l,2 + α l,3 Q l,3 = V l , (9) β l,1 Q l,1 + β l,2 Q l,2 + β l,3 Q l,3 = e x , (10)
γ l,1 Q l,1 + γ l,2 Q l,2 + γ l,3 Q l,3 = e y , (11)
where e x and e y are respectively the unity vector in the x and y directions. Equation [START_REF] O'rourke | An optimal algorithm for finding minimal enclosing triangles[END_REF] states that (α l,1 , α l,2 , α l,3 ) are the barycentric coordinates of the vertex V l with respect to the control triangle, while [START_REF] Speleers | Powell-Sabin splines with boundary conditions for polygonal and non-polygonal domains[END_REF] and [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] mean that (β l,1 , β l,3 , β l,3 ) and (γ l,1 , γ l,2 , γ l,3 ) are respectively the barycentric coordinates of the vectors e x and e y with respect to the control triangle.

In the case V l forms an angle different from π, if the control triangle has two sides aligned with the two boundary edges as in Figure 2 (left), from equation ( 9) derives that only one shape function has value different from zero in V l . This allows to directly set the boundary value to the unknown associated to the non-zero shape function. Moreover, from [START_REF] Speleers | Powell-Sabin splines with boundary conditions for polygonal and non-polygonal domains[END_REF] and [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] derives that the other two shape functions have zero derivative in the directions the two edges concurring in V l . In particular, using the notation depicted in Figure 2 (left) where triangle vertices are now called A, B and C, it results

α l,B = α l,C = 0, (⇒ α l,A = 1), β l,B γ l,B • t = 0, β l,C γ l,C • r = 0.
Denoting with U the component of U to be imposed, and û the L 2 projection of the prescribed values on the space of univariate quadratic splines on ∂Ω, the resulting equations relative to the unknowns of the vertex V l are where U l,A , U l,B , U l,C are the unknowns related to the shape functions with respectively triples (α l,A , β l,A , γ l,A ), (α l,B , β l,B , γ l,B ) and (α l,C , β l,C , γ l,C ).

U l,A = û(V l ), U l,A β l,A γ l,A • r + U l,B β l,B γ l,B • r = ∇û(V l ) • r, U l,A β l,A γ l,A • t + U l,C β l,C γ l,C • t = ∇û(V l ) • t, that is U l,A = û(V l ), U l,B = ∇û(V l ) • r -û(V l ) β l,A γ l,A • r β l,B γ l,B • r , U l,C = ∇û(V l ) • t -û(V l ) β l,A γ l,A • t β l,C γ l,C • t , r t V l A C B r V l A C B
In the case V l belongs to a plane angle, if the control triangle is aligned with the boundary as in Figure 2 (right), there is one shape function with zero value on V l , see equation [START_REF] O'rourke | An optimal algorithm for finding minimal enclosing triangles[END_REF]. Hence, the unknown related to this shape function is not affected by the boundary condition. Moreover, again equations ( 10) and ( 11) guarantee that the same shape function has zero derivative in the direction parallel to the boundary. Thus, with reference to Figure 2 (right), results

α l,A = 0, α l,B = 0, α l,C = 0, β l,A γ l,A
• r = 0, leading to the following relations for the unknowns in the vertex

V l U l,B α l,B + U l,C α l,C = û(V l ), U l,B β l,B γ l,B • r + U l,C β l,C γ l,C • r = ∇û(V l ) • r, that is U l,B = û(V l ) β l,C γ l,C -α l,C ∇û(V l ) • r α l,B β l,C γ l,C -α l,C β l,B γ l,B • r , U l,C = û(V l ) β l,C γ l,C -α l,C ∇û(V l ) • r α l,B β l,C γ l,C -α l,B β l,C γ l,C • r .

Stabilization by artificial viscosity

Stabilized formulations are obtained adding extra terms in the Galerkin weak form, with the goal of reducing the instabilities with the introduction of artificial viscosity. Here, two stabilization terms are considered.

The first term is the streamline upwind Petrov-Galerkin (SUPG) term, used to produce a stable upwind discretization without introducing excessive numerical dissipation. The SUPG method was introduced by T. Hughes in [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] and it is one of the most established stabilized formulations in finite element flow computations. The SUPG method introduces a certain amount of artificial viscosity in the streamline direction only.

The second stabilization term is a shock capturing term, and it is used to prevent oscillations around the discontinuities arising in the solution in supersonic regimes. These spurious oscillations might lead to severe accuracy loss or stability problems. Differently from the SUPG stabilization term, the shock-capturing stabilization term introduce an isotropic artificial diffusion, but only in a sharp zone surrounding the discontinuities.

Since a first order time discretization is employed, a simplified version of the SUPG stabilization term is used

D SU P G = N e e=1 Ωe A T ∇V . . τ SUPG A∇U dΩ e , (12)
which does not take into account the complete residual of the original equation, but only the convective term. The shock capturing stabilization term is

D SC = N e e=1 Ωe ∇V . . τ SC ∇U dΩ e . (13)
The terms τ SUPG and τ SC are stabilization matrices. Various options to compute these terms can be found in literature. In this work, for the SUPG stabilization term, two options have been retained: the first one is a simple diagonal matrix

τ SUPG ≡ τ SP = τ I,
and a constant value τ is used in the whole mesh, τ = dt/2, being dt the time step. The second choice, denoted τ T z , is inspired by the work of Tezduyar [START_REF] Tezduyar | Stabilization and shock-capturing parameters in SUPG formulation of compressible flows[END_REF] τ

T z = (1/2dt) (-2) + 3 r=1 3 q=1 u • ∇B (q) r (-1/2) .
For the shock-capturing operator, two choices described in [START_REF] Tezduyar | Stabilization and shock-capturing parameters in SUPG formulation of compressible flows[END_REF][START_REF] Wervaecke | A fully coupled RANS Spalart-Allmaras SUPG formulation for turbulent compressible flows on stretched-unstructured grids[END_REF] are retained, and are described using four sensors. The first sensor is based on the relative density variation and it is defined, in each element, as

[ ∇ρ] e = 3 r=1 3 q=1 ∇ρ • ∇B (q) r |∇ρ| .
The second sensor is based on the relative gradient according to the principle axes and is written

[ ∇U ] e = 2 d=1 ∂U ∂x d 2 1 2
, while the third sensor scales the with the values of the unknowns

[ U ] e = U .
Finally the fourth sensor takes into account the local relative convective flux

[A ∇U ] e = A∇U .
Hence, the first shock capturing artificial viscosity is defined as a combination of the sensors with a managing parameter β as follows

τ SC1 (β) = [ ∇ρ] e [ ∇U ] -(2-β) e [A ∇U ] -β e ,
and the second

τ SC2 (β) = [ ∇ρ] e [ ∇U ] -(2-β) e [ U ] (1-β) e [A ∇U ] -β e .
In [START_REF] Wervaecke | A fully coupled RANS Spalart-Allmaras SUPG formulation for turbulent compressible flows on stretched-unstructured grids[END_REF], it is suggested to use β = 1 is for smooth shocks and β = 2 for strong shocks. A compromise between the two definitions is proposed,

τ SC = τ SC (β = 1) + τ SC (β = 2) 2 ,
and it is the choice used in this work.

NUMERICAL TESTS

5.1 Numerical simulation of a supersonic flow past a forward facing step This is classical test case for 2D numerical codes for supersonic flows and a detailed description can be found in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. It consists in a supersonic flow entering an infinite long wind tunnel with a step. The interaction of the supersonic flow with the step and the tunnel walls creates a typical patter of shock reflections.

The numerical set up of the simulation is as follows (see Figure 3): the computational domain has dimensions 3 × 1 length units, and the step is 0.2 length units high and it is located at 0.6 length units from the left boundary. The flow is entering from the left boundary with a uniform horizontal Mach 3 velocity and pressure p = 1. The tunnel walls are non-slip boundaries, while the right boundary is a free exit. The simulation is started setting the variables at the inlet boundary values. No particular treatment is done for the corner of the step, which represent a singular point and it is known for introducing numerical errors, and no positivity correction was needed for the test cases here illustrated. However, with strong refinement in the corner, a negative pressure was developed in a point few elements after the corner and would require a special treatment. This point is still under investigation.

For the SUPG stabilization, the simple τ SM is used. For the shock capturing term, the form τSC1 is retained. The time step dt is linearly incremented at the beginning of the simulation, until reaching a final CFL=1. In order to compare the results with [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF], four computational meshes are employed with decreasing element size h = 1/40, 1/80, 1/160, 1/320. The results for a simulation time t = 4 are depicted in Figure 4 for the density variable. A correct shock reflection pattern is obtained and the results are globally similar to the reference [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF]. While the Mach steam is almost completely disappeared in the finest mesh, as expected, there is no Kelvin-Helmholtz instabilities near the top of the domain. This could be due to the extra diffusion introduced by the triangular discretization, with respect the quadrangular used in [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF], and it is still not fully understood.

In Figure 5 is also shown the shock profile on a horizontal line at y = 0.1, for the four meshes, with a comparison with the element size. It is possible to notice that there are no overshoots in the solution across the discontinuity and the shock is spread across about four element widths, which is acceptable.

The Double Mach reflection test

The second test case considered is the Double Mach reflection test, also introduced in [START_REF] Woodward | The numerical simulation of two-dimensional fluid flow with strong shocks[END_REF]. It represents the impinging of a strong shock wave, at Mach 10, against an oblique wall, which produce a typical reflected shock with a curved shape, and a bubble of denser gas. The set up of the test is the following, see Figure 6 : the computational domain is a box of dimensions 3 × 1, and a wall boundary condition is set in the lower boundary starting from x = 1/6. The right boundary and the remaining of the lower boundary (from x = 0 to x = 1/6) is set to an outflow boundary, while the left and the top boundaries are inlet boundaries. In the upper boundary, a time variable condition is set, simulating the traveling of the shock wave to the right. In the initial condition, the shock wave creates an angle of 60 • from the x axis and impinges the wall at x = 1/6. Post-shock conditions are set to the right of the wave, pre-shock to the left.

In order to avoid spurious oscillation of the solution in the first steps of the simulation, the fields are initialized with a smooth shock profile in the domain, with the discontinuity smeared over four elements. Moreover, the time step is increased with a cubic ramp until reaching a CFL=1 in 100 time steps.

For the SUPG stabilization, similar results are obtained for τ SM and τ T z . The shock-capturing stabilization however resulted more critical in this case: the parameter τSC1 used in the previous case resulted too diffusive, hence the form τSC2 is used instead. As in [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V[END_REF], four meshes are considered for comparison purpose, with element size h = 1/60, 1/120, 1/240, 1/480. Again, a stable solution is obtained with a sharp shock resolution. However, the shock patterns inside the bubble is not obtained, and further investigation are needed. A positivity correction was used in this case: the value of the internal energy is artificially set to 10 -5 in those Gauss points where a negative value is obtained.

CONCLUSIONS

In this paper is presented a finite element method based on Powell-Sabin (PS) splines, for the solution of the 2D compressible Euler equations in supersonic regime. The mathematical and geometrical tools to define the PS shape functions are presented, in particular, the subdivision of the elements that allows to define the PS splines as piecewise second order polynomials in each sub-triangle, the definition of the PS points and the control triangles.

A stabilized formulation is adopted, with SUPG and shock-capturing term. A simple isotropic artificial diffusion technique is adopted for the shock capturing term.

  are called Powell-Sabin B-splines and c l,r are the coefficients of the representation. Each B-spline B (r) l (x, y) can be seen as the solution of the interpolation problem (4) with all

Figure 1 :

 1 Figure 1: Powell-Sabin control triangle and relative Powell-Sabin points for the vertex V l .

Figure 2 :

 2 Figure 2: Control triangles for imposing boundary values: on an angle different from π (left) and equal to π (right).

FlowFigure 3 :Figure 4 :Figure 5 :

 345 Figure 3: Forward facing step case: computational domain and boundary conditions.

Figure 6 :

 6 Figure 6: Double Mach reflection case: computational domain and boundary conditions.

  Two classical examples are shown: the Woodward and Colella forward-facing step with Mach 3 flow and the Double Mach reflection test. Satisfactory results are obtained in terms of stabilization of the solution and resolution of the discontinuities. Further investigations are needed to understand the differences with the results with the reference solution and to set up an efficient positivity correction procedure. h = 1/60 h = 1/120 h = 1/240 h = 1/480

Figure 7 :

 7 Figure 7: Double Mach reflection problem: density ρ, 30 equally spaced contour lines at t = 0.2.

Table 1 :

 1 Condition number comparison for the stiffness matrix and the mass matrix obtained with optimal and non-optimal control triangles in a square domain: N e is the number of triangular elements of the mesh, h is the element size.Cond. numb. of the stiffness matrix Cond. numb. of the mass matrix

	Ne	h	optimal	non optimal	optimal	non optimal
	200 70.71E-03	83	389	171	4466
	800 35.36E-03	326	453	163	4466
	3200 17.68E-03	1305	1585	154	4469

Note that (2) can be expanded as p(ξ) = b

2,0,0 ξ 2 1 +b 0,2,0 ξ 2 2 +b 0,0,2 ξ 2 3 +2b 1,1,0 ξ 1 ξ 2 +2b 1,0,1 ξ 1 ξ 3 +2b 0,1,1 ξ 2 ξ 3 .