Visualization of Feature Separation in Advected Scalar Fields

TL;DR

This paper introduces a novel spatio-temporal visualization method for feature separation in fluid flows, using particle tracking and segmentation to better understand the evolution of scalar features over time.

Contribution

The paper presents a new visualization technique that segments feature volumes based on their contribution to separated features, improving analysis of topological changes in fluid flows.

Findings

Effective visualization of feature separation in multiphase flows

Application to liquid-gas and liquid-liquid flows demonstrates utility

Enhanced understanding of scalar feature evolution over time

Abstract

Scalar features in time-dependent fluid flow are traditionally visualized using 3D representation, and their topology changes over time are often conveyed with abstract graphs. Using such techniques, however, the structural details of feature separation and the temporal evolution of features undergoing topological changes are difficult to analyze. In this paper, we propose a novel approach for the spatio-temporal visualization of feature separation that segments feature volumes into regions with respect to their contribution to distinct features after separation. To this end, we employ particle-based feature tracking to find volumetric correspondences between features at two different instants of time. We visualize this segmentation by constructing mesh boundaries around each volume segment of a feature at the initial time that correspond to the separated features at the later time. To convey temporal evolution of the partitioning within the investigated time interval, we complement our approach with spatio-temporal separation surfaces. For the application of our approach to multiphase flow, we additionally present a feature-based corrector method to ensure phase-consistent particle trajectories. The utility of our technique is demonstrated by application to two-phase (liquid-gas) and multi-component (liquid-liquid) flows where the scalar field represents the fraction of one of the phases.

Figures

Tables

Table 1. Table I: Computation times for the three datasets. For Drop Collision (Drops) and Oil Inclusions (Oil), the measurement was performed on a desktop computer, and on the Cray system for the Non-Newtonian Jet (Jet) dataset. “Time” is the total computation time, “Adv. Step+ S 𝑆 S ” is the time for one step (done for each simulation step) of advection and extraction of S 𝑆 S . “ B 𝐵 B ” is time for extraction of boundaries (done once).

Dataset	#Steps	#Particles	Adv. Step+$S$	$B$	Total
Drops	$400$	$1.1\text{e}6$	$54.7\text{s}$	$5.8\text{s}$	$6.8\text{h}$
Oil	$500$	$1.1\text{e}6$	$54.8\text{s}$	$6.4\text{s}$	$8.3\text{h}$
Jet	$237$	$1.7\text{e}6$	$29.5\text{s}$	$6.1\text{s}$	$3.9\text{h}$

Table 2. Table II: Computation times for different value of parameter r 𝑟 r for Drop Collision and Non-Newtonian Jet datasets. Same notation as in Table  I .

Dataset	$r$	#Particles	Adv. Step+$S$	$B$	Total
Drops	$0$	$1.8\text{e}4$	$2.1\text{s}$	$0.3\text{s}$	$0.9\text{h}$
Drops	$1$	$1.4\text{e}5$	$7.9\text{s}$	$0.9\text{s}$	$1.5\text{h}$
Drops	$2$	$1.1\text{e}6$	$54.7\text{s}$	$5.8\text{s}$	$6.8\text{h}$
Jet	$0$	$1.7\text{e}6$	$29.5\text{s}$	$6.1\text{s}$	$3.9\text{h}$
Jet	$1$	$13.6\text{e}6$	$234.0\text{s}$	$46.2\text{s}$	$18.9\text{h}$

Equations

$$V_{t_0}^{t_F}(i,j)=\{\mathbf{x}:\mathbf{x}\in M_i(t_0)\wedge {\varphi }_{t_0}^{t_F}(\mathbf{x})\in M_j(t_F)\},$$

$$B_{t_0}^{t_F}(i,j)=\partial V_{t_0}^{t_F}(i,j).$$

$${ϵ}_p=\sum _{i\in I}(|{\mathbf{x}}_{p,i}^{\prime }-{\mathbf{x}}_{p,i}|+|{\mathbf{x}}_{p,i}^{\prime \prime }-{\mathbf{x}}_{p,i}^{\prime }|+|{\mathbf{x}}_{p,i}^{\prime \prime \prime }-{\mathbf{x}}_{p,i}^{\prime \prime }|),$$

Full text

1 Introduction

Flow visualization is widely used for the analysis of natural and engineering processes related to fluid motion. As simulation data continuously grow in size, mainly due to increasing computational power of both supercomputers and consumer desktops, feature extraction and tracking has gained particular attention. In feature-based visualization, visual data representation is reduced to important characteristics, typically physical quantities, flow topology, and quantities derived from generic scalar fields [30]. These features are then tracked in time to determine topological events, such as feature birth or split. The topology can be intuitively conveyed with a graph, where edges correspond to the features and nodes to the events. Such an approach reduces the amount of information to the essential part, relevant to the research at hand, and therefore allows for effective analysis of potentially enormous data.

While such feature-level visualization reveals the overall feature topology dynamics, detailed spatial information on volumetric partitioning of features is difficult to obtain with existing techniques. Thus, we propose a visualization approach that can reveal the separation dynamics within features at some initial time with reference to a later target time. To this end, we extract boundaries around regions within a feature that correspond to other features in the course of time, and additionally convey the temporal information of the separation by means of separation surfaces.

In the investigation of topology in general time-dependent flow, ridge extraction in the finite-time Lyapunov exponent (FTLE) is a standard method for the determination of regions that exhibit different flow behavior. The FTLE measures the separation of neighboring trajectories, and ridges in the FTLE separate those regions. With our visualization method, features representing physical quantities can be analyzed from a similar view point: given a feature at some point in time, we want to know how it will evolve over the course of time. However, whereas the ridges in the FTLE expose maximum separation of trajectories after a given time interval, our boundaries reveal separation of quantities after a time interval.

Our visualization technique expands upon traditional feature visualization in several ways. First, it allows for static visualization of dynamic processes, and therefore reduces visual clutter. Second, it combines advantages of standard feature tracking methods and FTLE-based methods, in that it allows for a detailed inspection of the separation of features. Third, it is applicable to a broad class of multiphase flow, including two-phase flow (liquid in gas surrounding) and multi-component flow (liquid-liquid or multi-component liquid in gaseous surrounding).

An important type of flow simulation is two-phase flow. Example cases are simulations of droplet collisions or liquid jets. The break-up of liquid jets into droplets appears in many technical applications as well as in nature. Examples range from a simple faucet over spray painting and spray drying in food processing to fuel injection in combustion engines. Especially in cases where the disruptive forces are much stronger than the cohesive forces, called atomization, the jets break up into many small droplets shortly after injection. Another important line of research is liquid-liquid flow, where immiscible fluids (for instance, oil and water) interact, for example, in emulsification processes. Such flow has often considerably different characteristics compared to two-phase flow, e.g., due to high viscosity of both components. Analysis of these processes, which occur on both very small temporal and spatial scales, is still very difficult.

Our visualization approach helps with the analysis of features in multiphase flow by visualizing the changing topology with respect to the feature at a selected time, and hence providing a static representation of the dynamics of features. To summarize, the main contribution of our work is a spatio-temporal visualization of feature separation dynamics that reveals the feature volume distribution among developed features as well as temporal information of the separation. As a necessary prerequisite, we propose a corrector scheme for particle integration in multiphase flow that ensures phase-consistent trajectories.

2 Related Work

The field of research most related to our work is feature tracking. Post et al. [30] provide a survey on the topic. Graph representation is a commonly employed approach to convey temporal evolution of features, e.g., in the analysis of combustion simulations [42] or combustion experiments [32], and the features can be abstracted by glyphs [31]. Often, features are defined by a threshold. Bremer et al. [8] propose a hierarchy-based approach that alleviates the dependency on predefined thresholds. For the inspection of features at a given time, a linked-view approach is employed in these methods. Gu and Wang [16] compute time-dependent state transition probabilities for volumetric data and visualize a 3D view of the volume together with a 2D graph representation of the transitions. Grottel et al. [15] visualize the evolution of molecular clusters on a timeline as an addition to a 3D representation of the molecules to monitor the quality of the clustering. Laney et al. [22] track bubbles and visualize their split and merge behavior as a graph. Further research on feature tracking concentrates on clustering methods, including the work by Ozer et al. [29], where user-defined feature characteristics are used to determine feature groups. A survey on graph representations is provided by Wang [41]. Our method differs from these approaches, in that we focus on the spatio-temporal aspect of feature representation, where the topology is directly conveyed in the feature volume.

For the analysis of the combustion process in engine simulation, Garth et al. [14] proposed a set of visualization methods that operate on time-varying unstructured grids. A recent work by Sauer et al. [35] utilizes particle and volume data from the simulation runs to track features over longer time intervals. In our method, we insert particles after the simulation run. On the one hand, this gives a better control over the particle density and therefore the detail level, on the other hand, however, it necessitates some corrector schemes for phase-consistent advection in multiphase flow, as presented in this paper.

Delocalized quantities are employed in unsteady flows to statically investigate the dynamics of scalar fields [36], where quantities are averaged over time along integral curves. Our concept can be viewed as a special case of this, where only the values at trajectory end points are stored at the seed point. The visualization of delocalized quantities can be improved using uncertainty information from the FTLE field [33]. For the topology analysis in general unsteady flow, the FTLE has become a standard visualization method [17]. Sanderson [34] proposed an alternative inspired by Lyapunov exponents that is based on traveled particle distance instead of particle divergence.

The dynamics of fluid flow is often represented by surfaces or volumes. In the work by van Wijk [40], stream surfaces are obtained by extracting isosurfaces from a scalar field. This in turn is generated by advecting streamlines backward up to the boundary with predefined scalar values and resampling those values along the streamlines. Becker et al. [4] use flow volumes in unsteady flow to reveal the flow dynamics around regions of interest. An implicit version of the flow volumes by Xue et al. [43] allows for more detailed inspection of the flow.

In our work, the labels corresponding to the volumetric contributions can be interpreted as different materials. Material interface reconstruction is a challenging topic in visualization, especially with respect to volume preservation in multi-material configurations [2, 26, 7]. For the extraction of volume boundaries, we adhere to the representation by the marching cubes algorithm [25], since in our case the material distribution within the cells is explicitly given by points with assigned labels. Moreover, as the number of materials (i.e., labels) can be arbitrarily large for a given investigation, it would complicate the computation of the interface based only on the fraction information. Further research in material interface analysis includes the work by Obermaier et al. [27], where the stability of reconstructed interfaces is investigated by comparing it with time surfaces.

For the visualization in multiphase flow, we use a corrector method to ensure that the advected particles remain in the given phase during tracing. Related to this problem are approaches for surface tracking. Stam [38] proposed a method for the computation of interface velocities to properly translate surfaces. Bojsen-Hansen et al.[5] developed a method for tracking surfaces undergoing topology changes, without prior information on the underlying physics. In the work by Bonneel et al. [6], Lagrangian transport is used for correct displacement interpolation. Solomon et al. [37] optimize the transportation in terms of Wasserstein distances, which allows for efficient shape transformation. In our work, we analyze the evolution of volumes, and computing Wasserstein distances in this case would be computationally prohibitive and also would not guarantee physically correct correspondences.

We exemplify the utility of our technique using multiphase flow simulation datasets. Investigation of droplets is an active area of research [44, 39, 11], and so is the analysis of liquid jets [13, 24, 10]. We refer the reader to Fuster et al. [12] for a detailed introduction to multiphase flow simulation and to Lefebvre [23] for a thorough description of liquid atomization and sprays.

3 Simulation Data

The data used in this paper are computational fluid dynamics (CFD) simulations of two-phase (liquid-gas) and two-component (liquid-liquid) flow, based on the Navier-Stokes equations and discretized on rectilinear grids [9]. The velocity is represented by a cell-based vector field $\mathbf{u}(\mathbf{x},t)$.

In case of scientific multiphase flow simulations, the volume of fluid (VOF) method [19] is typically used for interface tracking, where an additional volume fraction field $f(\mathbf{x},t)$ is maintained for each cell. A cell contains only the gas phase when $f=0$, only the liquid phase when $f=1$, and the interface is located within cells with $0<f<1$. For liquid-liquid flow, the gas phase is replaced by the surrounding liquid phase. To maintain sharp interfaces, piecewise linear interface calculation (PLIC) [45] is used which approximates the interface by planar patches. In this paper, we use linear interpolation in time for tracking, whereas for sampling of the fraction field, we use the PLIC interface to determine the sampled value. Only features that undergo advective transport are considered, since, e.g., diffusive features would be difficult to capture with the particle-based approach, and are therefore beyond the scope of our work.

Our technique employs particle advection in order to determine the volumetric correlation of features between a reference time step $t_{0}$ and the target time step $t_{F}$. Usually, only a certain fraction of timesteps of a simulation are stored for later analysis. Using these data with reduced time resolution can lead to errors in the estimation of trajectories, especially in multiphase flow, where particles potentially leave the initially assigned phase for this reason. This problem could be avoided if particles were advected during the simulation, which is increasingly popular, as reported by Sauer et al. [35]. Additionally, particles could be densely populated using the method proposed by Agranovsky et al. [1] in order to capture more details of the topology of feature dynamics. In our datasets, however, particle data was not provided. Nevertheless, to still ensure robust particle advection in terms of phase consistency in multiphase flows, we introduce a feature-based corrector method, as described in Section 5.6, which utilizes the $f$-field in a corrector step during particle advection to ensure that particles stay in the respective phase throughout integration.

4 Visualization Method

We define a feature $M$ as a region with the value of the quantity $f$ exceeding a corresponding threshold $\tau$: $M=\{{{\mathbf{x}}}:f({{\mathbf{x}}},t)>\tau\}$. In our experiments, $f$ represents the fraction of fluid in a cell in multiphase flow and $\tau$ is set to [math], so that the liquid phase (or one of the phases in liquid-liquid flow) represents the features. For other scalar quantities, $\tau$ could be set to a specific value of interest. Feature separation is illustrated in Figure 1, where an initial feature $M_{1}$ at time $t_{0}$ splits within the time interval $]t_{0},t_{F}[$, resulting in three features $M_{\{1,2,3\}}$ at time $t_{F}$.

4.1 Separation-Based Feature Segmentation

With our visualization technique, we want to capture how features develop topologically in the course of time. For instance, if a feature splits into two new ones, we want to learn how the volume of the initial feature is divided among them, or, in other words, what the volumetric contributions of the new features in the initial one are. To accomplish this goal, we have to determine the spatio-temporal correspondence between the feature $M_{i}(t_{0})$ at time $t_{0}$, and features $M_{j}(t_{F})$ at some other time $t_{F}$. That is, for each $M_{j}(t_{F})$, we define its volumetric contribution within $M_{i}(t_{0})$ as

[TABLE]

where $\phi_{t_{0}}^{t_{F}}({{\mathbf{x}}})$ maps the initial position ${{\mathbf{x}}}$ at time $t_{0}$ to its position at time $t_{F}$, as it is advected by the flow. We refer to $\Delta t$ = $t_{F}-t_{0}$ as the computation time interval. For our visualization, we extract the closed boundary of the volume $V$:

[TABLE]

Note that for $\Delta t=0$, $B$ bounds the whole feature volume if $i=j$, and $B$ does not exist if $i\neq j$. Also, note that our method allows $t_{F}$ to be earlier in time than $t_{0}$ if we set $\Delta t<0$. In Figure 1, the method is illustrated for a simple case where a feature splits into three.

Volumetric contributions $V_{t_{0}}^{t_{F}}(i,j)$ can be composed of disconnected segments, either due to a merge followed by a split of the initial features, or due to disjoined volume segments inside the initial feature that together form a separate feature. As will be shown in Section 5.4, our visualization method allows one to readily identify and analyze these cases.

Separation Surfaces

To support temporal analysis of feature separation, we supplement our technique with the extraction of temporal separation surfaces, denoted $S$. While the above method finds the volumetric correspondences in the initial time $t_{0}$ to the features at target time $t_{F}$, here, the separation surfaces are constructed within the whole interval $]t_{0},t_{F}]$. They divide the volumetric contributions in the initial feature as the corresponding features split into new features.

This is illustrated in Figure 2, where the initial feature $M_{1}(t_{0})$ has split into two at time $t_{1}$, and the resulting feature $M_{2}(t_{1})$ has further separated into new features $M_{2}(t_{2})$ and $M_{3}(t_{2})$ at time $t_{2}$. In Figure 2, both split events are illustrated by the separation surfaces $S_{t_{0}}^{t_{1}}$ and $S_{t_{0}}^{t_{2}}$ that divide the volume segment according to the volumetric contributions of the feature resulting from the separation.

To create the separation surfaces and determine the time at which they occur, changes in contributions are detected within small time intervals. Specifically, for each time increment $\delta t=t_{2}-t_{1}$, where $\delta t\ll\Delta t$, contributions $V_{t_{0}}^{t_{2}}(i,k)$ in each feature $M_{i}(t_{0})$ are computed and compared with the previous contributions $V_{t_{0}}^{t_{1}}(i,j)$. If the number of contributions $V_{t_{0}}^{t_{2}}(i,k)$ for which $V_{t_{0}}^{t_{1}}(i,j)\cap V_{t_{0}}^{t_{2}}(i,k)\neq\emptyset$ is greater than one, a separation has occurred in the given interval inside the volume contribution $V_{t_{0}}^{t_{1}}(i,j)$, and the separation surface is generated in $M_{i}(t_{0})$ where different $V_{t_{0}}^{t_{2}}(i,k)$ adjoin.

Theoretically, the temporal evolution of feature separation could be accomplished by repeatedly extracting the closed boundaries $B_{t_{0}}^{t_{k}}(i,j)$ for varying $k$. This would, however, lead to repeated construction of overlapping boundaries that would be difficult to analyze. The separation surfaces $S$, on the other hand, provide an open structure that complements the extracted boundaries $B$.

5 Numerical Approach

In our method, we operate on both the Eulerian frame, in which the simulation data is defined, and the Lagrangian frame, where the inserted particles represent the feature volume and allow us to determine the volumetric correspondences.

Our visualization framework works as follows. The input to our approach are time series of a vector field ${\mathbf{u}}$ for particle advection and a scalar field for feature definition. The initial time $t_{0}$ is chosen, at which particles are seeded within $M_{i}(t_{0})$. These particles are then advected up to time $t_{F}$, whereby integration is performed between consecutive simulation steps, and so is the construction of separation surfaces. At time $t_{F}$, the connected components that define the features $M_{j}(t_{F})$ are determined, and for each advected particle, the surrounding feature is found. The feature labels are then mapped to the seed points, and finally, boundaries are extracted for each label within $M_{i}(t_{0})$. In Figure 3, the algorithm is presented as a diagram. In the following sections, each step, marked on the diagram with the corresponding section number, is described. Please see also the supplementary material for a detailed pseudo-code of the algorithm.

5.1 Particle Seeding and Advection

At time $t_{0}$, particles are seeded at ${\mathbf{s}}_{p}$ if $f({\mathbf{s}}_{p},t_{0})>\tau$. This is illustrated in Figure 4, where $\tau$ is represented by the PLIC interface that divides the cell into liquid and gaseous phase, and seed points are located only at positions ${\mathbf{s}}_{p}$ enclosed by the PLIC patch.

To determine if ${\mathbf{s}}_{p}$ lies within the phase of interest, we compute the seed position relative to the PLIC patch. Therefore, we have to compute the patch orientation and position. We first compute the gradient of the $f$-field at the cell center ${{\mathbf{x}}}_{c}$, and then normalize it to obtain the PLIC normal ${\mathbf{n}}$. Next, we determine the attachment point ${\mathbf{a}}$ as the most distant cell node from the PLIC patch in the direction of ${\mathbf{n}}$. The PLIC patch position $l$, i.e., the distance from ${\mathbf{a}}$, is obtained iteratively by ensuring that the volume enclosed by the patch equals the value $f$ in the given cell. Then, we compare $l$ with the distance $d$ from ${\mathbf{s}}$ projected onto the normal. If $d<l$, a particle is positioned at ${\mathbf{s}}_{p}$. In the figure, the distance $d>l$ for ${\mathbf{s}}_{m+3}$ implies wrong phase, and hence, no particle will be seeded there. Please note that for anticipated application to different data models, trilinear interpolation can be used instead of PLIC reconstruction to determine the threshold position.

The seed position ${\mathbf{s}}_{p}$ and the maximum number of seeds per cell $n_{c}$ are controlled by a global refinement parameter $r$ in the form $n_{c}=(2^{d})^{r}$, where $d$ is the data dimension ($d=3$ in our datasets). For $r=0$, the seeds are positioned at the centers of the simulation cells. For $r>0$, the cells are recursively divided $r$ times into equally sized subcells, and the particles are positioned at the centers of these subcells.

To compute the flow map $\phi_{t_{0}}^{t_{F}}({\mathbf{s}}_{p})$, we set the particles $p$ at the seed points and advect them up to time $t_{F}$ (Figure 1) using the standard fourth-order Runge-Kutta scheme performed between consecutive simulation time steps in the interval $[t_{0},t_{F}]$. The seed points ${\mathbf{s}}_{p}$ are stored for later processing to determine the correspondence between features at times $t_{0}$ and $t_{F}$.

5.2 Feature Labeling

To extract the features at time $t_{F}$, we first create a grid with the same structure as the input simulation grid. We equip this grid with a bit mask, where cells $c$ with $f({{\mathbf{x}}}_{c},t)>\tau$ are set to $1$ and the others to [math]. With this bit mask, we can find connected components that correspond to the features $M_{j}(t_{F})$. To this end, we employ region growing. Here, two feature cells $c_{m}$ and $c_{n}$ are considered connected if they are face neighbors. Each feature is then identified by a label $j\in\mathbb{N}$. The connected components are stored in a label field, which again structurally corresponds to the simulation grid. This greatly simplifies particle assignment, as described in the next section. All grid cells without connected components are marked with an invalid label (i.e., $-1$). In Figure 1, the labeled features are marked red, green, and yellow. Feature labeling can be performed independently of the advection step, since determination of feature labels only requires the scalar field $f({{\mathbf{x}}},t_{F})$.

5.3 Particle Label Assignment

In the next step, we have to determine those particles $p$ for which $\phi_{t_{0}}^{t_{F}}({\mathbf{s}}_{p})\in M_{j}(t_{F})$. For each advected particle $p$ at time $t_{F}$, we determine by which feature $M_{j}$ it is bounded (Figure 1) and assign the corresponding label $j$ to the particle. If it does not lie inside any feature, a non-valid label (i.e., $-1$) is assigned to it. Since we operate on rectilinear grids, we can find the bounding cell by iterating over the node coordinates $x_{i}$ of the grid until $x_{i}<x_{p}<x_{i+1}$, with $x_{i}\in\{x,y,z\}$. The feature label $j$ is then mapped to the stored seed points ${\mathbf{s}}_{p}$ corresponding to the advected particles (shown by dashed arrows in Figure 1). Additionally, to visualize the temporal evolution of the features, the intermediate particle positions are stored during the computation of the visualization, and the computed labels are assigned to them at this stage. For scalar quantities beyond multiphase flow, the value $f$ at the particle position can be trilinearly interpolated. If the particle is inside a cell without a label but the trilinear interpolation results in $f>\tau$, we can compute the gradient of $f$ and assign the label of the nearest cell with $f>\tau$ in gradient direction.

5.4 Boundary Extraction

The labeled particles allow us to extract the volume boundaries $B_{t_{0}}^{t_{F}}$ within the feature $M_{i}(t_{0})$ that correspond to the features $M_{j}(t_{F})$, according to Equation 2. For each label $j$, we find all seed points ${\mathbf{s}}_{p}$ with this label and compute a bounding box around them. Subsequently, a rectilinear grid is generated within the bounding box such that the seed points coincide with the grid nodes. At each grid node, we set a value $1$ if a seed point is actually located at this position, and [math] otherwise. This way, we can employ the standard marching cubes algorithm to extract the boundaries. Since we require that the grid node positions correspond with the seed points, the size and number of cells is implicitly controlled by the refinement parameter $r$. In Figure 1, the boundaries $B$ are marked by darker curves that bound the volumes $V$.

The disconnected feature segments (Section 4) can be detected either directly from the color of the boundaries (each feature maps to a different color label of the volume regions, and disconnected volume regions belonging to the same feature have the same color) or using the particles from intermediate time steps that also carry the feature label information.

5.5 Extraction of Separation Surfaces

To extract the separation surfaces $S$ described in Section 4.1, after each advection substep, feature labeling (Section 5.2) and particle label assignment (Section 5.3) are performed to obtain the current feature segmentation. The resulting seed labels are stored for two subsequent time steps $t_{k}$ and $t_{k+1}$. For each label at $t_{k}$, the corresponding seeds and their labels at $t_{k+1}$ are analyzed. More than one unique label at $t_{k+1}$ indicates that the segment has split in the current time interval, and therefore, separation surfaces are extracted by performing a modified marching cubes algorithm for each $2$-combination of the labels. In the algorithm, node values are set to “$+$” and “$-$” which correspond to either of the labels, and the surface passes through the edge centers between “$+$” and “$-$” nodes. To ensure open surfaces, the outer nodes (i.e., not belonging to any of the two labels) are marked with an invalid negative value. In the final step of marching cubes, the edges of each triangle are tested for whether they lie on edge with the outer node, and discarded in this case. Finally, the surfaces and the corresponding time stamps $t_{k+1}$ are stored for visualization.

5.6 Phase-Consistent Trajectories in Multiphase Flow

Since typically only a fraction of simulation time steps is saved for analysis and visualization (due to potentially high storage demands of large simulations), and because multiphase flow is highly nonlinear (due to, among others, surface tension forces), there is usually not enough information to ensure phase consistency of the advected particles, i.e., the particles may stray off to the other phase during advection. However, assuming that physically-based phase transitions do not occur in the simulation (which is the case in our applications), the particles must stay in the same phase throughout advection. To accommodate this requirement, and to provide plausible results with the available information, we introduce a three-stage approach that utilizes the scalar field representing the phase of interest to correct positions of stray particles during advection.

Particle Corrector

After each integration step between $t_{k}$ and $t_{k+1}$, we test each advected particle if it remained in the assigned phase. If this is not the case, a three-stage procedure is applied to translate the stray particle back to its original phase. In the first stage, the neighborhood of the particle is searched for valid particles (i.e., particles that remain in the assigned phase at $t_{k+1}$). That is, in the $3\times 3\times 3$ cell neighborhood at time $t_{k}$ (i.e., before the current advection step) the nearest neighbor (if any) is chosen. Its displacement vector $\mathrm{d}{\mathbf{x}}_{p}={\mathbf{x}}_{p}(t_{k+1})-{\mathbf{x}}_{p}(t_{k})$ is applied to the stray particle, i.e., ${\mathbf{x}}^{\prime}(t_{k+1})={\mathbf{x}}(t_{k})+\mathrm{d}{\mathbf{x}}_{p}$ (Figure 4). Since this does not guarantee that the particle will be back in the original phase, in the second stage, the nearest cell with $f>0$ is searched (using a loop over neighboring cells), and the particle is translated to the boundary of this cell along the line connecting the particle with the cell center, as shown in Figure 4, where the new particle positions are denoted ${{\mathbf{x}}_{p}}^{\prime\prime}$. Afterwards, in the third stage, for each particle in the interface cell, the position of the PLIC patch is calculated at the current simulation time step. Then, it is tested if the particle position lies within the volume enclosed by the PLIC patch. If not, the particle is translated along the line ${{\mathbf{x}}_{p}}-{\mathbf{a}}$ to the patch (Figure 4).

6 Implementation

We have implemented our visualization method as a plugin in the ParaView visualization framework [3]. To improve the visual representation of the boundaries $B$ and separation surfaces $S$, we smooth the meshes extracted by the marching cubes algorithm. Although some error is introduced in the process, a smoothed representation improves perception, and hence we consider it a reasonable trade-off. For the visualization of the labeled features in multiphase flow datasets, we extract the geometric representation of the PLIC interfaces using the method by Karch et al. [20] and resample the labels on the PLIC patches. Apart from the boundaries we also store intermediate particle positions after every $n$th simulation time step ($n=8$ in our applications), in order to visualize temporal evolution of the features, see for instance Figure 5.

For the second stage of the particle corrector scheme (translation to a cell boundary), we adopted the line-box intersection algorithm proposed by Kay and Kajiya [21].

6.1 Parallelization

For large datasets, parallelization is necessary due to large memory requirements of the simulation data. For instance, in the Non-Newtonian Jet simulation (Section 7.3), each simulation time step consists of over $20\,\textrm{GiB}$, which could not be processed on a regular desktop computer. Hence, we employ data parallelism (i.e., the approach adopted by ParaView), where the data is split among several processes (possibly running on different machines), and each process works on a preassigned subdomain. In this approach, explicit communication is necessary for data exchange.

Due to the global nature of particle advection, those particles that leave a subdomain must be sent to the processes managing the subdomains the particles have entered. To avoid cases where the Runge-Kutta substeps are computed across neighboring subdomains, we store ghost cells around each subdomain. Another global algorithm we employ in our technique is connected component labeling. Our implementation is based on the method proposed by Harrison et al. [18]. Additionally, the transfer of particle labels to their respective seed points must be handled explicitly. That is, for each particle we store its id within the original subdomain, and the process id responsible for that subdomain. This information is then used to transfer the labels to the correct processes and to save the label at the correct position in the label array.

In addition to the process-level parallelization, we also employ thread-level parallelism for particle advection, i.e., we employ the OpenMP [28] parallel-for loops to speed up iteration over particles being advected.

7 Results

We demonstrate the utility of our method on three datasets from direct numerical simulations of incompressible flow. The Drop Collision and Non-Newtonian Jet datasets are two-phase flow simulations, where liquid and gas phases occur simultaneously. The Oil Inclusions dataset contains oil inclusions in water surrounding. For Drop Collision and Oil Inclusions, $r=2$, whereas for Non-Newtonian Jet, $r=0$.

7.1 Drop Collision

The first dataset is a simulation of a two-phase flow with two water droplets that collide peripherally. This dataset consists of $256^{3}$ cells, which cover a domain size of $1\,\text{cm}^{3}$, and $461$ time steps. In the top row of Figure 5, selected simulation time steps (from $t_{0}$ to $t_{F}$) are shown. We analyze the features from time $t_{0}=0.478\,\text{ms}$ (i.e., $60$th time step, just before the collision) up to time $t_{F}=2.85\,\text{ms}$ (i.e., $461$st time step, with small droplets resulting from the collision). The two merged drops form a flat shape that splits into an inner disk and an outer ring. Finally, both structures separate into small droplets.

Our visualization in the bottom row reveals the volumetric correspondences between the volumes of the two initial drops and the drops at the final time step. It can be seen that many small droplets are formed from narrow sections that extend radially from the collision center. These small drops, however, form only on the sides, whereas the volume sections at the bottom and top contribute to relatively large drops.

In Figure 6, a section of the right drop from Figure 5 has been cut out to reveal the inner structure of the volume contributions. Interestingly, almost all droplets that develop after the collision are formed from volumes that are close to the collision center and propagate outward to the back side. This structure bears some similarity to cracks in a solid material that propagate from the impact point. In Figure 6, separation surfaces $S$ are shown. Here, the breakup of inner disk and outer ring is apparent in the form of cone like structures in both droplets. The surfaces within the cones have similar colors, indicating simultaneous separation of multiple droplets.

For this dataset we have additionally applied our method in reverse time direction, i.e., we computed the correspondence of the last simulation time step to the first one to reveal how exactly the initial two droplets contributed to the later smaller droplets. The result is shown in Figure 7. In Figure 7, the two droplets at time $t_{F}=0.478\,\text{ms}$ are visualized by the PLIC interface. They are colored by the connected component labels, and their contributions in the droplets at time $t_{0}=2.85\,\text{ms}$ are visualized in Figure 7. As can be seen, the volume from the red drop dominates on the upper right part of the droplets at later physical time. Figure 7 shows a selected droplet that rotates after separation from the disk-shaped drop. The boundary from the blue droplet is transparent to reveal the inner structure of the drop. The interface between the two boundaries $B$ forms an S-shape, and there is distinguishable symmetry of the two volumes. The investigation of the interplay between volume distribution and rotational motion of the feature could shed light on the details of drop dynamics, this is, however, out of scope of this paper.

7.2 Oil Inclusions

The next dataset is a simulation of colliding oil inclusions in water surrounding. For this dataset, $f$ represents volume fraction of oil. The simulation domain consists of $256^{3}$ cells covering $1\,\text{cm}^{3}$, and $501$ time steps that span a time interval of $12.5\,\textrm{ms}$. In Figure 8, selected simulation time steps are shown. The deformation of the oil inclusions differs strongly from the water inclusions in the first dataset due to the presence of liquid surrounding. Even before the collision, the initially spherical inclusions deform strongly and they lose a considerable amount of momentum due to friction. In the final time steps after the collision, the inclusions quickly form spherical shapes and their velocity drops to almost zero.

Figure 9 shows the boundaries $B_{t_{0}}^{t_{F}}$ with the upper left part of the front drop cut out to reveal the inner segmentation. The structure is quite distinct from the one in the previous dataset, as it is less ordered and contains many small segments in the central part.

The separation surfaces for the same configuration are shown in Figure 9. Interestingly, in the upper segment, there is a closed surface (marked with red box) that indicates split followed by merge event. With our technique, we can further investigate this part by selecting the particles comprising the enclosing segment (and therefore, a developing inclusion), as visualized in Figure 10. Here, the first time step was selected just before the split, the second shows the small part that separated (red box) and the last three show the re-merging process. In Figure 10 we identified the particles that belong to the splitting part and shown them for all stored time steps (i.e., every $5$th simulation step), with time color-coded. The bend parts clearly indicate the split and merge (red boxes).

This dataset demonstrates the ability of our technique to perform a detailed inspection of the separation processes.

7.3 Non-Newtonian Jet

The last dataset is a simulation of the injection of a jet made up of a non-Newtonian shear thinning aqueous solution of Praestol 2500 (0.3% weight) into air at $p_{\text{atm}}=1\,\text{bar}$. The jet is introduced into a rectangular computational domain with a diameter of $\textrm{D}=1.2\,\text{cm}$ and with the velocity profile from a short nozzle with a mean velocity of $u=75\,\text{m/s}$, and therefore a Reynolds number $Re=19\,000$. The domain has a size of $42\,\textrm{D}$ in the direction of the injection and of $10\,\textrm{D}$ in each of the directions of the injection plane. The domain is discretized over $2688\times 512\times 512$ cells. The temporal development of the jet is sketched in Figure 11. At $t_{0}=1.60\,\textrm{ms}$, the jet has just entered the computational domain and has started to interact with the surrounding gas. The resulting effects include the bending back of the jet tip and the development of surface waves on the jet core. Due to the high Reynolds number, the jet becomes unstable very quickly, and at $t_{F}=3.97\,\textrm{ms}$, the core has expanded substantially in radial direction and has mostly broken up into ligaments and droplets—the atomization has begun.

Figure 12 shows the boundaries $B$ and separation surfaces $S$ at $t_{0}$. Small polygons have been removed for clear overview. Apparently, in the investigated time interval, the front of the jet does not undergo strong disintegration (except for the tip), whereas in the rear, many segments imply the origins of elongated droplets that develop later. The separation surfaces at the bottom of the figure are uniformly distributed, although small droplets detach earlier, as indicated by small dark blue surfaces.

To better understand the atomization process, we use the proposed visualization method to investigate the origin of ligaments within the liquid core. In Figure 13, we can observe the spatial origin at $t_{0}$ for a selected number of ligaments, and their position at $t_{0}+\Delta t$ is displayed in Figure 13. In Figure 13, the three features are in close proximity in stream-wise $x$-direction. This is surprising, as their original liquid mass in Figure 13 was distributed over a rather large range in $x$-direction. In Figure 13, the liquid mass of the ligaments is scattered within the jet core. This is especially noticeable for the droplet colored in green, which can be seen in the zoomed cutout. This indicates a strong influence of the flow field inside the jet core, as the liquid mass has to combine before separating from the core.

In Figure 14, the temporal evolution of one selected feature (colored) is shown with the entire jet (gray) in the background. The temporal evolution is from $t_{0}$ at the top to $t_{0}+\Delta t$ at the bottom. The development in the first three time steps confirms our assumption from Figure 13. The liquid mass, which will later form the ligament, is initially scattered over one third of the jet length. With time, the liquid mass near the jet surface is being decelerated due to the interaction with the surrounding atmosphere. This allows the rest of the mass, which started out more central in the jet, to catch up, and by the third time step, the bulk of the liquid mass is located within the first surface wave, while a few parts are located in the next wave. From the underlying jet, we can observe that the waves are already starting to deform strongly from the shape of a periodic wave. At the next time step, nearly all of the liquid mass has merged and moved even closer toward the surface. The liquid mass has reached its most compact form. In the next time step, the surface wave has disintegrated and the liquid mass has become a complex branching ligament. With time, the ligament moves further in radial direction, interacting more strongly with the atmosphere and therefore is stretched in $x$-direction, until we can observe the first droplet breaking off the ligament in the final time step.

With the proposed visualization method, we were able to make interesting observations. Contrary to the impression we get from Figure 11 that the jet is mainly expanding in radial direction, we observed in Figure 13 that features which are in close proximity in $x$-direction at $t_{0}+\Delta t$ originate from rather distant positions at $t_{0}$. From this, we can infer that the complex nature of the velocity field leads to features moving at higher velocity in $x$-direction. We could also observe that this is even true for the velocity field within the jet core. Contrary to the rather solid seeming core in Figure 11, both Figure 13 and especially Figure 14 present how the flow field inside the core is creating the analyzed features from liquid mass scattered throughout the core. Furthermore, Figure 14 is an excellent representation of the different steps of primary jet break-up showing the formation of surface waves, the deformation of the waves, the rupturing of the waves into ligaments, and the stretching and break-up of the ligaments into droplets.

The presented method allows for a simple and intuitive way to analyze the spatial and temporal development of the primary break-up of liquid jets. Thanks to this method an insight into the movement patters of select features of the break-up process could be gained which would otherwise have been hard and unlikely to obtain.

7.4 Performance

The computation of the visualization for the Oil Inclusions and Drop Collision datasets was performed on a commodity desktop with Intel i7 3.6 GHz processor (4 cores) and $32\,\textrm{GiB}$ RAM. For the Non-Newtonian Jet dataset, we have used $64$ nodes on a Cray XC40 System, each with two Intel Xeon E5 processors and $128\,\textrm{GiB}$ of memory, and one process running per node. Table I shows the configurations (i.e., number of simulation time steps and the number of seeded particles) and timings (i.e., total time, average time for one advection step and extraction of $S$, and time for extraction of $B$) for the three datasets. Although all datasets have comparable numbers of particles, the computation of the Non-Newtonian Jet dataset is not much faster, despite the parallelization. This is due to the fact that with the data parallelism in the multi-process configuration, the particles are not evenly distributed among the processes. Moreover, the exchange of particles across the subdomains introduces some overhead after each integration step.

We have also investigated the performance dependency on parameter $r$. As Table II shows, in case of Drop Collision, the difference between $r=0$ and $r=1$ is relatively small arguably due to better usage of thread-level parallelism. For $r=2$, on the other hand, the number of particles increases the computation time substantially. For the Non-Newtonian Jet, it can be seen that with $8$-fold increase in particle number, there is $8$-fold increase in time for advection step and the total time increases about $5$ times, suggesting that the interprocess communication constitutes considerable time which is less dependent on the number of particles.

For the visualization we have used a commodity desktop for all datasets, since the data produced by our technique (i.e., the particle data and mesh for boundaries $B$) does not exceed $1\,\textrm{GiB}$, even for the Non-Newtonian Jet dataset.

7.5 Analysis of the Corrector in Multiphase Flow

To quantify the amount of applied correction that is necessary to keep the particles within the initial phase in multiphase flow (see Section 5.6), we accumulate the introduced displacement for each particle over all time steps $i$:

[TABLE]

where the first term is the translation with the displacement vector of the nearest phase-consistent particle (Figure 4), the second term in the sum corresponds to the translation to the nearest cell with $f>0$ (Figure 4), and the third term is the translation to the PLIC patch in the interface cells (Figure 4). For the analysis, we visualize the value $\epsilon_{p}$ at the seed positions ${\mathbf{s}}_{p}$. In Figures 15 and 15 the displacement $\epsilon_{p}$ is shown for the Drop Collision and Non-Newtonian Jet datasets, respectively. The particles that were translated most correspond with the boundary positions. This is easy to interpret, since the particles that most probably leave the liquid phase are close to the feature interface, i.e., at unstable regions. Although the maximum value of the accumulated correction is relatively large ($0.07$ of the domain size for Drop Collision and $0.06$ of the domain size for the Non-Newtonian Jet), it occurs only for small number of particles, and the correction for the majority of particles is considerably smaller (as indicated by pale blue color in the figures).

Robustness

In Figure 16 we also investigate the sensitivity of the particle advection to temporal resolution of data and different variants of the particle corrector. The selected droplet region from the Drop Collision dataset is particularly prone to error due to the collision which is difficult to capture with integration. For the analysis, we use the temporal resolution of $50$ (left column) and $400$ (right column) simulation time steps spanning the same time interval. The top row shows the resulting boundary reconstruction for particles with no correction applied, in the middle row, the particles were moved to the closest cell and consequently to the PLIC patch without the adjustment by the displacement vector. In the bottom row, full correction is applied. With no correction, even with high temporal resolution, there are regions containing particles that left the initial phase, as indicated by dark blue surfaces. Interestingly, for two-stage correction, the stray particles are gone, but in Figure 16(c), some particles at the front are incorrectly assigned. With full correction, little difference can be seen between both resolutions. This shows that our method can provide reliable results also for lower temporal resolution, although in our experiments we used higher temporal resolution to further increase reliability of the results.

8 Discussion

As already mentioned in Sections 3 and 5.6, there are two problems with tracking particles in multiphase flow. First, reduced temporal resolution of the simulation data does not allow for accurate integration due to the dynamics of phases. Particle advection is particularly challenging where phases collide, as their momentum changes rapidly. Second, the involved surface tension forces lead to very complex behavior, e.g., in regions where inclusions disintegrate. Several methods have been proposed to find mesh and volume correspondences. However, those methods either work only on 2D manifolds, or provide correspondence based on optimization problems that typically cannot guarantee physical correctness, e.g., by ignoring the underlying velocity field. For instance, the Wasserstein distances may result in physically inaccurate distant correspondences. Moreover, applied to volume, the correspondences from one cell usually scatter among many cells, which additionally complicates the choice of the correct particle path.

Therefore, we devised a new correction method which, admittedly, still has some limitations. Due to the translation to the nearest cell, this method introduces asymmetry in the advection, i.e., the particles have different trajectories depending on the integration type (forward or backward in time). In fact, we have observed that for the reverse advection in the Drop Collision dataset, the drop in the bottom right in Figure 7 is identified as coming from only the red drop in Figure 7, whereas in the forward advection, we have observed that both drops contributed to it. We have investigated several approaches to solve this problem, which, however, due to its complexity, remains beyond the scope of this work. In the future, we would like to further pursue possible solutions to the asymmetry.

9 Conclusion and Future Work

In this paper, we have proposed a novel visualization method for the analysis of the separation of features. Whereas the boundaries $B$ determine feature volumes, the separation surfaces $S$ provide the information on the time point at which a given feature separation has occurred. Hence, both methods complement each other to provide a comprehensive space-time visualization of feature separation.

For the analysis of multiphase flow, we have introduced a corrector method that utilizes the scalar field representing one of the fluid phases in order to ensure phase-consistent particle trajectories. We have demonstrated the utility of our method for liquid-gas and liquid-liquid flow datasets. However, it can be as well applied to more complex multiphase simulations (e.g, with gas and two liquid components), since each phase is usually described by a separate scalar field, and therefore allows for direct application of our method. Similarly, transitions between liquid and solid states could be investigated, where both states are stored as separate field data and therefore can be used to distinguish between the phases.