Flow visualization and field line advection in computational fluid dynamics: application to magnetic fields and turbulent flows

Accurately interpreting three dimensional (3D) vector quantities output as solutions to high-resolution computational fluid dynamics (CFD) simulations can be an arduous, time-consuming task. Scientific visualization of these fields can be a powerful aid in their understanding. However, numerous pitfalls present themselves ranging from computational performance to the challenge of generating insightful visual representations of the data. In this paper, we briefly survey current practices for visualizing 3D vector fields, placing particular emphasis on those data arising from CFD simulations of turbulence. We describe the capabilities of a vector field visualization system that we have implemented as part of an open source visual data analysis environment. We also describe a novel algorithm we have developed for illustrating the advection of one vector field by a second flow field. We demonstrate these techniques in the exploration of two sets of runs. The first comprises an ideal and a resistive magnetohydrodynamic (MHD) simulation. This set is used to test the validity of the advection scheme. The second corresponds to a simulation of MHD turbulence. We show the formation of structures in the flows, the evolution of magnetic field lines, and how field line advection can be used effectively to track structures therein.


Introduction
The visual display of three-dimensional (3D) steady and unsteady vector fields output from numerical simulations is an active research area among the visualization community. While the display of scalar quantities is to a large degree a solved problem, with mature algorithms readily available and frequently employed by scientists in a variety of disciplines, including computational fluid dynamics (CFD), numerous challenges remain for vector data. Techniques that work reasonably well for 2D vector fields fall far short with the addition of a third spatial dimension: ambiguous interpretation of geometric shapes projected onto a 2D display, visual clutter, and algorithmic performance are but some of the problems encountered with 3D vector field display.
In this paper, we briefly survey current approaches to vector field visualization with a focus on those methods best suited for high-resolution computational physics simulations. Comprehensive, more general surveys may be found in [1]- [3]. We next describe in detail methods that we have implemented as part of an open source visual data analysis tool, VAPOR, whose general capabilities were previously reported in this journal [4]. We also present a novel technique, termed field line advection, that we have developed for portraying the dynamics of magnetic fields (or the vorticity field in the case of neutral fluids) under the influence of advection in a velocity field. Lastly, we report our experiences using these techniques in the investigation of two sets of simulations. We use two simulations of the viscous and ideal magnetohydrodynamic (MHD) equations with the same initial conditions to validate the field line advection scheme. Then, we study magnetic field lines and their advection by the fluid in a 1536 3 simulation of MHD turbulence.

Survey of vector field visualization methods
A variety of diverse methods exist for discrete vector field visualization. We follow the taxonomy put forth by Post et al [3], classifying these techniques into four areas described below. 4 turbulent flows the number of critical points can be enormous and some form of reduction, such as the hierarchical methods proposed by de Leeuw [20] or Bremer et al [21] are required. To the authors' knowledge these techniques have not been applied to unsteady 3D fields.
In this paper, we describe several geometric techniques we have implemented in the VAPOR package [4]. In particular, streamlines are used to convey instantaneous information about magnetic or velocity fields, while pathlines are available for investigation of the evolution of unsteady flow fields. We augmented these techniques with a novel, combined steady/unsteady field line advection technique. Our reasons for selecting geometric methods are primarily due to their familiarity and acceptance by the scientific community. Nevertheless, as noted briefly above, a number of issues arise with the application of these methods to the simulation outputs of the type we are interested in. Here we expand further upon these issues, and later we describe our implementation in detail as it relates to the problems noted. We also discuss the solutions and workarounds we have employed to mitigate these difficulties.

Geometric field visualization challenges
Three principal challenges arise when employing geometric methods for 3D time varying volumes, particularly those with highly complex dynamics and computed on high resolution computational grids: interactive performance, avoiding visual clutter, and effective seed placement.
Unless a costly, highly parallelized visual supercomputing environment is available, maintaining interactivity is difficult. IO, in particular, presents a substantial bottleneck; three component fields must be read from disk for each time step. A 1024 3 grid of 32-bit floating point values requires approximately 10 s to read a single time step assuming a storage system capable of delivering data at 100 MB s −1 . Once the data are resident in memory (assuming sufficient primary memory is available) the integration algorithm itself may be computationally expensive depending on the length of integration and number of integral objects displayed. The position of points along a curve must be found by solving the ordinary differential equation: dp where p, v, and t are respectively position, velocity, and time. Integrating equation (1) yields Equation (2) must be evaluated numerically. Simple Eulerian methods have been shown to lack the accuracy needed [22], and typically higher order, more expensive methods such as Runge Kutte are employed. Equally as important as the numerical integration method is the choice of step size. If too large, successive points will skip grid cells, missing features in areas of strong gradients. If too small, needless computation will be performed. Finding cells containing the current integration step, particularly for less regular grids, may also be expensive, possibly requiring an exhaustive search of the domain. The algorithm must also interpolate vector components to arbitrary spatial and temporal coordinates. Lastly, coordinates or field values may need to be transformed back and forth between physical and computational space. All of these operations must be performed at each integration step and for each curve.
In addition to their computational expense, the generation of the large numbers of curves necessary to convey global field information is likely to result in an unusably cluttered image. Reducing the number of curves will of course decrease the clutter, but may make it impossible to portray a global sense of the field or may lack coverage necessary to ensure the capture of significant phenomena. Recent research has led to algorithms that adaptively control the integration length and density of curves [23]- [25] in an effort to reduce screen clutter while preserving information content, but these methods are only applicable to steady fields (streamlines).
Closely related to the issue of visual clutter is the effective placement of seed points-the initial starting position of an integral object-within the field. Random or regular distributions of seed points of sufficient density to capture all important features result in excessive clutter as noted earlier. If dense placement strategies are combined with adaptive curve length and density control methods, reducing the homogeneous distribution of display primitives, important areas of convergence or divergence may be obscured. A number of feature-sensitive and densityguided seed placement algorithms have been reported in the literature aimed at ensuring good coverage while reducing the number of seed locations generated through dense distributions. Verma first introduced a 2D strategy for placing seed points in the vicinity of critical points using templates depending on the critical point-type [26]. This method has been extended to 3D by Ye et al [27]. Random seed placement may also be biased by physical properties of the data, such as regions of high local vorticity as reported by Sanna et al [28]. But with large numbers of critical points or other phenomena, as expected with highly turbulent flow, these approaches will not guarantee a sparse curve distribution.

Vector field visualization in VAPOR
With the aforementioned challenges in mind we describe the capabilities and features of the geometric vector field visualization system that we have implemented in VAPOR. The principal means for ensuring interactive performance is by way of a wavelet-based multi-resolution data representation as reported by us previously [4]. Variables and fields in the data can be accessed at full resolution, or the resolution can be reduced by powers of two in all three spatial dimensions. Further, the data model supports rapid data subsetting allowing fast access to small subregions of the volume at high resolution (or the full domain at coarse resolution).
The computational time of flow integration is usually dominated by the time required to load the vector field from disk. Only when many thousands of flow lines are displayed at very high accuracy does the integration and/or rendering time dominate. Consequently, it is desirable to explore flow lines at lowered resolution or in small spatial domains before performing massive high-resolution flow calculations. Users are able to accelerate the scientific discovery process by formulating hypotheses in an interactive mode, based on low-resolution or spatially limited flow integration. Such hypotheses can be validated using increased resolution at possibly noninteractive computation rates.
Interactivity is also related to integration efficiency. In VAPOR, flow integration is performed using fourth order Runge-Kutta integration. During integration, the vector field is determined at arbitrary points in the 3D volume by performing tri-linear interpolation of field values within the grid cells. For unsteady (time-varying) flow integration, the vector field is determined at arbitrary temporal points by linearly interpolating the field between its spatially interpolated values at the nearest available time values. An adaptive integration algorithm (similar to the algorithm used by Kenwright and Lane [29]) is utilized to enable the integration to proceed rapidly in regions where the vector field is nearly constant, but to achieve high accuracy Top: field lines resulting from seed points randomly distributed throughout a volume containing a magnetic flux tube (left), and resulting from a distribution biased toward increased magnetic field strength (right). Bottom: a volume rendering of a current sheet (left), and field lines resulting from seed points selected with a 2D probe depicting current intensity.
when the vector field rapidly changes direction from one cell to the next. The integration is performed by successively selecting points along a curve determined by field magnitude and direction at the most recent sample point. The step size (distance between successive sample points) is adapted (repeatedly doubled or halved) to ensure that the change in angle between successive field samples is between 3 • and 15 • . Users are provided with an accuracy parameter that determines the minimum and maximum step sizes, which are by default smaller than the grid cell size. In very turbulent datasets, the field magnitude can increase by several orders of magnitude in passing from one cell to the next. To maintain accuracy with such data, the adaptive step size is halved multiple times until the resulting step size satisfies the required accuracy.
The problem of seed placement is addressed in VAPOR with two techniques. First, the distribution of random seed points can be biased by the magnitude of a field. Consider, for example, the magnetic flux tube in figure 1 (top, left). The flux tube occupies a very small subset of the entire volume, so that random seed placement with a uniform seed distribution pattern would poorly sample or entirely miss the flux tube. In many situations like this, flow 7 seed placement is improved by biasing the random distribution so as to position most seeds where a field magnitude is very large or very small. Figure 1 (top, right) shows the magnetic field lines that result from biasing the seed distribution toward increased magnetic field intensity.
A second seed point placement strategy supported by VAPOR is the explicit manual specification of nonrandom, individual seed point coordinates, which may be calculated outside of the VAPOR environment and supplied to VAPOR as a list of points, or may be selected through an interactive data probe. The data probe, which can be manipulated by a GUI, provides an arbitrarily oriented 2D slice of the data, visualized using simple pseudo coloring. Seed points can thus be placed precisely in the proximity of interesting structures, identified with the aid of visual inspection. For example, in figure 1 (bottom) it is desired to position flow seeds slightly above and slightly below the current sheets, a task easily accomplished via the interactive probe depicted.

Field line advection
The need for understanding the dynamics of one vector field under the influence of a second one arises in several areas of CFD. As examples, we will discuss hydrodynamic (barotropic) flows, where in the inviscid case the vorticity is advected by the velocity field, and ideal MHD flows, where the magnetic field is advected by the velocity in the absence of Joule dissipation (which, in general, applies to any flow for which an equivalent of Kelvin's or Alfvén's theorem exists [30,31]). Here, we describe a technique for visualizing the time-varying behavior of a vector field, assuming its dynamics result from advection imposed by a second field. We refer to this technique as 'field line advection' because it amounts to advecting entire field lines (not simply points) in a time-varying velocity field.
A seemingly straightforward formula for advecting a field line in a velocity field from time T to time T + 1 would be the following: (i) perform a steady integration, determining a field line, L, at time T ; (ii) select a sequence of points along L; (iii) advect the sequence of points at time T in the velocity field to time T + 1, resulting in a sequence of points at time T + 1; (iv) connect the sequence of points at time T + 1 to form a new line, L .
Under the conditions stated above (existence of a form of Kelvin or Alfvén's theorems, and ideal flows), the preceding recipe gives a procedure to advect a sequence of points. Unfortunately the resulting line, L , may not trace an actual field line if, for example, the system has nonzero dissipation. To overcome this deficiency we utilize a method proposed by Nordlund [32] as follows: (i) perform steady integration at time T , as above; (ii) select only one representative point P from L, based on a selection criterion described below; (iii) advect the point P in the velocity field to time T + 1, resulting in a point Q at time T + 1; (iv) perform steady integration at time T + 1, using the point Q as the seed. 8 For MHD flows, Nordlund proposed that the representative point at time T (in step (ii)) be obtained by maximizing the magnetic field intensity along the field line. The significance of this choice for MHD is explained in the next section. VAPOR provides several modifications of this algorithm for special situations as described below.
The advected point Q may be outside the region where the field values are known (an event whose likelihood is increased when small subregions are explored as facilitated by VAPOR). An alternative field line advection mode for this situation is to consider a set of points {P 1 , P 2 , P 3 , . . . , P N } along the field line and advecting the entire set to time step T + 1, resulting in points {Q 1 , . . . , Q M }, where Q i are the advected points remaining inside the domain bounds. A seed point Q is selected from the points Q i as the point having the largest field magnitude at time T + 1. We term this scheme prioritization after advection.
Another modification was proposed by Nordlund to deal with situations where field magnitude oscillates along field lines, resulting in two or more local maxima. In this case it may be desirable to only advect a small portion of the field line containing a local maximum of the field intensity, rather than using the global maximum of field magnitude along its full length. The user can select a lower cutoff value C, so that in step (ii) above, the search for the representative P is stopped if the field magnitude is less than C times the field magnitude at the seed point. This has the effect of shortening the advected field lines so as to reflect the local behavior of the field.
VAPOR also extends the scheme by allowing users to choose arbitrary fields for steady integration (steps (i) and (iv)), unsteady advection (step (iii)), and prioritization (step (ii)). As is shown in the next section, this latter option is useful to advect field lines in viscous (dissipative) flows choosing points for the advection for which the effect of dissipation is a minimum.

The MHD equations
In this section, we show the field line visualization capabilities of VAPOR, and apply the field line advection scheme discussed in the previous sections. To this end, we consider two sets of numerical simulations of the incompressible MHD equations. These equations describe the evolution of the velocity and the magnetic field in a conducting flow. In particular, the equations describe the evolution of the fields in planetary cores [33] and in experiments using liquid sodium or gallium [34]. To a good approximation, they also describe the large-scale dynamics of plasmas in stars, the interplanetary and interstellar medium, and the Earth's magnetopause and magnetosphere.
The MHD equations derive from Maxwell's equations for sub-relativistic velocities (as a result, the displacement current can be neglected). They consist of the momentum equation for the velocity field, with the addition of the Lorentz force to take into account the back-reaction of the magnetic field in the flow, and the induction equation for the magnetic field. In dimensionless alfvénic units the equations read [31] ∂v ∂t Equation (4) is accompanied by ∇ · b = 0, indicating the lack of magnetic monopoles; and equation (3) is accompanied by the incompressibility condition ∇ · v = 0. For an incompressible flow, the pressure P in equation (3) is obtained by solving a Poisson equation that follows from the condition ∇ · v = 0. The fields v and b are respectively the velocity field and the magnetic induction; j = ∇ × b is the current density, ν is the kinematic viscosity, and η is the magnetic diffusivity.
In the ideal case (ν = η = 0) the total energy is conserved by the equations. The ideal MHD equations have nonlinear exact solutions u = ±b, for which all nonlinear terms in equations (3) and (4) vanish. These solutions correspond to Alfvén waves. The 'alfvénicity' of the flow (alignment between the velocity and the magnetic field) can be measured by the following quantity, known as the cross helicity, This quantity is also conserved by equations (3) and (4) in the ideal case. In the following sections it will be convenient to measure the relative cross helicity ρ C = 2H C ( |v| |b| ) −1 , which takes values between −1 and 1, and can be interpreted as the cosine of the angle between the two fields. The induction equation (4) expresses simple physics: magnetic field lines in a conducting fluid behave as material lines, which are advected by the velocity field. From equation (4) the following relation for the time evolution of the magnetic flux through an orientable surface S can be proven (see e.g. [31] This result is known as Alfvén's theorem. In the ideal case (η = 0) the magnetic flux is conserved. If we consider a magnetic flux tube as the aggregation of magnetic field lines through a small closed curve, every curve embracing the flux tube moves with the fluid. As a result, the magnetic field lines are advected by the velocity field. We note, however, that the amplitude of the field can change as a result of the deformation of the surface S: if the fluid motions stretch the surface S and change its area, the magnetic field amplitude has to change to maintain the magnetic flux constant. From the fact that the topological structure of the magnetic field cannot change with time because the field lines are frozen to the fluid, it follows that in ideal MHD flows a third quantity is conserved, the magnetic helicity: where A is the vector potential, defined as b = ∇ × A.
The dynamics of magnetic fields have many similarities with hydrodynamic turbulence [31,35]: as an example, we can recall the Batchelor analogy between vorticity and induction. Indeed, in the hydrodynamic case Kelvin's circulation theorem for barotropic flows is equivalent to Alfvén's theorem for MHD flows. As in the MHD case, when Kelvin's theorem holds, individual vorticity (ω = ∇ × v) field lines behave as material lines and vorticity field lines move with the fluid advected by the velocity (see e.g. [30]). The conserved quantity associated with this property of vorticity field lines is the kinetic helicity, For simplicity, in the following we will focus on MHD examples of field line visualization and advection.
In the presence of dissipation (nonzero viscosity in the hydrodynamic case, or nonzero magnetic diffusivity in the MHD case) the advection of the field lines is broken: the presence of dissipation allows field lines to diffuse away [31]. As a result, the topological structure of the field changes as the flow evolves, and the connectivity of the field lines can change: the field lines can 'reconnect'. There is no unambiguous way to follow a field line in time, and advecting all points will give a line that will not correspond to the field line at a later time, as was discussed in section 2. Attempts have been made to define a 'viscous velocity' or a 'velocity of slip' in such cases, but if such a velocity could be defined conservation of the helicities and of the connectivity of the field lines would follow for this new velocity, as pointed out by Moffatt [31].
In this context, it is worth pointing out that VAPOR does not make any assumption about how the velocity is defined, and that the user can make different choices, e.g. by computing a posteriori a new field to use as the velocity for the advection scheme. In some cases, although non-unique, the introduction of a velocity of slip can be useful. An obvious choice for MHD would be to use as the local velocity for the advection the velocity at which the local frame of reference electric field vanishes, given by u = E × B/B 2 , where E is the electric field and ∇ × E = − ∂B ∂t [32]. As a result, in all of our experiments described below we use the prioritization after advection scheme presented in section 2: using the velocity v we advect five uniformly distributed points at time T and prioritize the points remaining in the domain at time T + 1, selecting a single point for the generation of a new line at time T + 1. The obtained line is a field line (within numerical errors) at each time. Choosing a good criteria (small local dissipation) to select the point P at time T + 1 ensures that advection dominates over diffusion, and that the field line can be followed for a reasonable amount of time. The method also allows for changes in the topology of the field lines, since only one advected point results in a new field line computed after each advection step.
The MHD equations develop, from smooth initial conditions, thin current sheets. These are the regions in which magnetic reconnection takes place and where the connectivity of magnetic field lines changes, while magnetic energy is rapidly dissipated. The local dissipation in MHD goes as ηj 2 , and current sheets are characterized by large local values of the current density. In 2D [36,37] regions of strong current density and reconnection are also associated with neutral x-points of the magnetic field, and quadrupolar structures in the vorticity. In 3D the geometry of these regions is more complex and sheets can be formed without the need for magnetic null points [38]. This observation motivated a possible prioritization criterium: the point P is chosen as the maximum of the magnetic energy along the field line. In general, this gives a point that is far from the current sheet, and where advection dominates over the dissipation. Another possible criteria for prioritization is to choose P at the minimum of the current along the magnetic field line. Similar considerations can be done for hydrodynamic flows.
The advection scheme was validated in simple MHD flows at low resolution, for which the evolution of magnetic field lines is well known (e.g. the Orszag-Tang vortex [38]). In the following two sections, we present a comparison of field line advection in ideal and viscous simulations (section 3.2), and advection of magnetic field lines in a turbulent flow (section 3.3). In both cases, the MHD equations were integrated using a pseudospectral method [39,40] to compute spatial derivatives, and Runge-Kutta of second or fourth order to evolve the equations in time. The pseudospectral methods are of very high order, have exponential convergence, and it is worth noting here that these methods are conservative and non-dispersive [41], so the ideal simulations have no numerical dissipation up to the time when the smallest scales resolved by the grid are excited. All the analysis was done interactively, even in cases for which the amount of data exceeded tens of terabytes.

Comparison of advection in ideal and viscous flows
In this section, we compare two runs with the same initial conditions, the first ideal (η = ν = 0), and the second with nonzero viscosity and magnetic diffusivity. The initial conditions are given by a generalized MHD Taylor-Green flow [42] v = v 0 (sin x cos y cos zx − cos x sin y cos zŷ), b = b 0 (cos x sin y sin zx + sin x cos y sin zŷ − 2 sin x sin y cos zẑ).
The initial velocity field corresponds to two counter-rotating eddies with a shear layer between them. The initial magnetic field is proportional to the curl of the velocity. This flow has several symmetries that are preserved by the MHD equations as time evolves, and that can be imposed in a code to save memory and computing time [43,44]. Simulations need only be done in a box of size [0, π/2] 3 [43], while the dynamics in the rest of the [0, 2π ] 3 box can be reconstructed from the symmetries. The simulations used here for the comparison correspond to a total of 2048 3 grid points [42], although smaller domains need only be visualized because of the symmetries. Details of the simulations and of the time evolution of this flow can be found in [42].
Here, we focus on the comparison of how field lines are visualized and advected by VAPOR in the ideal and the viscous case. To this end we consider one simulation with η = ν = 0 (run A), and another simulation with η = ν = 1.3E − 4 (run B). As previously mentioned, the MHD equations tend to develop from smooth initial conditions sharp structures (sheets) in the current. The magnetic field changes direction fast around these structures, where most of the Joule dissipation takes place. As a result, the ideal run can only be continued until the smallest scale resolved by the grid is reached: as soon as the thickness of the current sheet reaches this scale, the simulation is stopped. This happens for t ≈ 2.6. We compare the two runs for t between 2.45 and 2.55. These times are of special interest because at t ≈ 2.45 structures in runs A and B are still thick and the flows are similar, while for later times thinner structures develop quickly and the resistive effects in run B start to become relevant.
We start by analyzing field lines and studying the behavior of the field line advection scheme in VAPOR in run A. For this run, the advection scheme is exact (except for numerical errors) for the magnetic field. Figure 2 shows the current density intensity and magnetic field lines in a small subvolume at four different times: t ≈ 2.45, 2.48, 2.52 and 2.55. A thin sheet-like structure can be observed in the current intensity. Seeds for the magnetic field lines at t ≈ 2.45 were densely distributed along a line perpendicular to the sheet. At t ≈ 2.45, the magnetic field is strong above and below the region with strong currents (indicated by the gray field lines), and is weak in the sheet (white lines). Above the structure the magnetic field runs in a certain direction and as it crosses the structure it rotates and changes direction in 90 • . Note the anticorrelation between regions of strong current (where dissipation dominates in viscous runs) and regions of strong magnetic field. The magnetic field lines for t > 2.45 were obtained by advecting the field lines at t ≈ 2.45 using the prioritization after advection method described in section 2, prioritizing points with maximum magnetic intensity.
As time evolves, magnetic field lines above and below the structure are advected by the velocity field towards the current sheet. This is clear at t ≈ 2.55; note how the field lines above the sheet (initially distributed along a vertical line), now get closer to the structure and displaced from the vertical line. The field lines in the current sheet change with time more slowly and after t ≈ 2.48 do not seem to move. This can be understood by inspecting the velocity field in the same region (figure 3). Above and below the structure in the current intensity, the velocity field points towards the structure. However, near the structure the velocity field points along the structure, thus explaining the evolution of the advected magnetic field lines observed in figure 2. However, in the middle of the current sheet the velocity field is nonzero and large, giving rise to the question, addressed below, of how the magnetic field lines in this region can remain almost static.  Left: visualization of current intensity and magnetic field lines in run A at t ≈ 2.45 from the corner of the subvolume. Note the double current sheet, and the abrupt change in the direction of the magnetic field lines above and below the structure. Right: relative cross helicity, ρ C , in the same region. Magenta corresponds to ρ C > 0.9; only regions with ρ C > 0.75 are shown, regions not aligned are transparent. Figure 4 shows another visualization of the current intensity and magnetic field lines in the same structure, with a rendering of the relative cross helicity ρ C from the same viewpoint. This quantity measures the alignment between the velocity and the magnetic field, and takes values between −1 (when the fields are anti-aligned) to 1 (when the fields are aligned). The region in the middle of the current structure has maximum alignment, with ρ C > 0.9. There are two other 'alfvénic' regions above and below the structure. The alignment between the velocity and the magnetic field in the middle of the structure explains the slow evolution by advection of the magnetic field lines in the current sheet. When the velocity field goes along the magnetic field line, the magnetic field line is not displaced. These results show that the advection scheme in VAPOR is consistent with the physics of the ideal MHD equations. The behavior observed in the evolution of the magnetic field lines in run A (see figure 2), and indicated by the advection scheme in VAPOR, corresponds to the expected dynamics of magnetic field lines near a current sheet in MHD flows [36]- [38], [45]: magnetic field lines point in different directions above and below the structure, and the magnetic field intensity is in general small in the current sheet. The velocity field pushes magnetic field lines towards the sheet, where in the dissipative case magnetic energy can be dissipated and the magnetic field lines can reconnect.
A comparison of the field lines at t ≈ 2.55 in run A (ideal) and B (viscous) is shown in figure 5. The field lines at this time were obtained (as in figure 2) by advecting the same seeds starting from t ≈ 2.45 for both runs. Only small differences can be observed between the ideal and viscous runs: most of the differences in the advection of the magnetic field lines are concentrated in the vicinity of the current sheet. As previously mentioned, for a viscous flow there is no unambiguous way to track a field line, and results of the field line advection should be interpreted with care. However, the similarities between both runs and the deviations in regions of strong dissipation are remarkable. This is the region where the dynamics of the magnetic field is affected by Joule dissipation, and the prioritization after advection method with prioritizing by minimizing the local dissipation ensures that even after the advection the magnetic field line corresponds to an actual field line. The mean quadratic distance between the field lines for runs A and B as time evolves is shown in figure 6. The distance is computed as where N is the total number of points in all the advected field lines, i is the index that labels each point, and x (A) and x (B) are respectively, the Cartesian coordinates of the field lines in runs A and B at a given time. The average quadratic distance increases slightly, as resistive effects become more relevant in run B and the evolution of the two runs depart. However, the average departure is small. If the departure is computed for each field line in both runs (characterized as the field lines computed from the same seed), larger departures are observed in the region with large current intensity (and dissipation), while departures become smaller in the rest of the box. Figure 6 also shows a histogram of the absolute value of the differences between the magnetic field intensity in the visualized subvolume at t ≈ 2.45 and 2.55 (the values are normalized by the magnetic field intensity in the ideal run). Although at early times differences are small, at t ≈ 2.55 the development of a tail in the histogram with large deviations (of order one) is observed. These differences develop in the region of strong dissipation. Differences in the velocity field at the same time are smaller (by a factor of ten), while differences in the current density are larger (also by a factor of ten).
Another example of current visualization and the result obtained after doing field line advection from the same seeds in runs A and B is shown in figure 7. In this figure, we started at t ≈ 2.45 from a set of seeds placed just above and just below the region with strong current, and some seeds placed in the region of maximum current. The same seeds were advected in VAPOR using data from runs A and B. As with the previous set, the region with strong current was observed to become thinner as the field lines were pushed together towards that region. At t ≈ 2.55 the region with strong current is thicker in run B (the effect of dissipation), and the largest deviations between the field lines are again observed for field lines in the current sheet.

Advection in a turbulent flow: tracking structures
The results presented so far indicate that field line advection can be a powerful technique to visualize the evolution of magnetic fields in MHD (or vorticity in hydrodynamics) and to understand the dynamical evolution of a flow. The flows studied so far have a simple geometry, with a simple (quasi-2D) current sheet configuration and magnetic field lines pushed towards the structure, where dissipation occurs if the flow is viscous/resistive. However, in astrophysics and geophysics Reynolds numbers (the ratio of nonlinear to dissipative terms) are typically large and the flows are unstable. Magnetic fields are then embedded in a medium which is turbulent. As a few examples, planetary magnetospheres, the solar wind, and the convective region in the sun can all be viewed as laboratories for turbulence, with numerous satellites that observe its properties [46]- [48]. In the solar wind, the magnetic energy spectrum is close to the Kolmogorov spectrum E M (k)∼k 5/3 [49], whereas in the magnetosphere of Jupiter, the spectrum (determined using Galileo spacecraft data) is anisotropic [50] and follows a scaling law ∼k −2 ⊥ (where k ⊥ refers to wavevectors perpendicular to the local mean field).
As discussed in section 1, the visualization and analysis of turbulent flows present several extra difficulties. Visualization of field lines in 3D turbulence suffers badly from visual occlusion, because flows are complex and the structures are intermittent and tend to form dense clusters [51]- [53]. Methods that detect singular or 'interesting' points in the flow also result in cluttered images if no form of hierarchical reduction is used, and to the best of our knowledge these methods have not been used to track structures in highly complex flows such as those that result from numerical simulations of fully developed turbulence. In a previous work [4], we discussed several methods available in VAPOR to interactively pick regions in a flow and place seeds for field line integration at a fixed time. These methods proved to be useful to identify structures in turbulent flows (see e.g. [45]). In this section, we consider how the field line advection scheme can be used to extend these capabilities by giving a way to visualize the evolution of a turbulent flow and to track its structures.
For the analysis, we use data stemming from a 1536 3 simulation of MHD turbulence [45,54]. As in the simulations discussed in the previous section, equations (3) and (4) are integrated using a pseudospectral method in a periodic cubic box of side 2π . However, in this case there are no symmetries imposed in the initial conditions nor in the equations, and as a result there are no savings in storage or in computing time. The initial conditions are where it is assumed complex conjugate terms are added to have real fields, v ABC denotes a helical Arn'old-Beltrami-Childress (ABC) flow (see e.g. [55]) and φ v,k and φ b,k are random phases with Gaussian distribution, chosen to have divergencefree velocity and magnetic fields. The initial kinetic and magnetic energy spectra are ∼ k −3 exp[−2(k/k 0 )] 2 . The system decays freely in time. The kinematic viscosity and magnetic diffusivity are ν = η = 2 × 10 −4 . More details of the initial conditions and a study of the evolution of the system can be found in [45,54].
Here, we will focus our analysis on a small structure that appears in the flow as turbulence develops. The structure was identified in [45] and used for studies of visualization of turbulent flows in [4]. It corresponds to a 'current roll': a cylindrical structure with strong current intensity formed as a current sheet folds and rolls as the result of an instability akin to the hydrodynamic Kelvin-Helmholtz instability. The structure is shown in figure 8. Seeds for the magnetic field lines were placed in the vicinity of the current sheets. As in the examples discussed in the previous section, the magnetic field changes direction sharply above and below the sheets. In the roll (and in a second rolled structure on the top left corner of the subvolume) magnetic field lines go along the roll. Velocity field lines integrated from the same seeds are also shown in figure 8. The velocity field points on the average in a direction perpendicular to the sheets (except for local alignments of the fields in the sheets). As a result, we expect these structures to be advected and deformed by the velocity field as time evolves.  ≈ 1.6). The roll can be seen in the center of the subvolume (260 × 160 × 200 grid points). Right: current intensity and velocity field lines, integrated from the same seeds. Figure 9 shows the result of advecting the magnetic field lines up to t ≈ 2.8 using twelve snapshots of the fields. As time evolves, the roll is advected and deformed, and the current sheets in the region studied are stretched. The current sheet initially in the bottom of the subvolume is stretched until it almost occupies the entire region. It is remarkable that the advection scheme tracks these structures, and that magnetic field lines that were initially placed just above and below the current sheets move together with the sheets. In particular, the current roll initially in the center of the region is deformed as time evolves. In the first three visualizations, the structure and the field lines can be seen moving together. However, at late times (see the last visualization) the deformation and the decay in the intensity of the current is such that the cylindrical structure cannot be identified anymore. However, the bundle of magnetic field lines associated with it can still be seen, pointing to the position where the deformed structure is located.
It is worth noting that the evolution of the current sheets and the evolution of the magnetic field lines are obtained from two methods that are different. The current intensity is visualized directly using the solution data that results from integrating numerically the MHD equations, while the magnetic field lines are advected using the prioritization after advection scheme applied to the solutions of the MHD equation. Moreover, even in ideal MHD the current intensity is not frozen to the velocity field. The current is the curl of the magnetic field, and regions of strong current are associated with sharp gradients in this field (sharp changes in its amplitude and/or direction). As the magnetic field lines are advected, the gradients displace and the current sheets move. The overall behavior of both visualized quantities is consistent (although the evolution of each quantity is obtained by different methods), and compatible with the physics of the MHD equations. The example shown here is more complex than the previous example, not only because the flow is turbulent, but also because the evolution of the flow is faster, more dissipative, and current sheets change their position and geometry rapidly with time. Figure 9. Visualization of current intensity and advected field lines at four different times (from left to right and top to bottom) starting from t ≈ 1.6 up to 2.8. Twelve snapshots of the fields were used in the integration. The associated movie (available from stacks.iop.org/NJP/10/125007/mmedia) shows the temporal evolution. Figure 10 shows the same snapshots, but visualizing only the current intensity. The complexity of the turbulent flow can be clearly seen in the last images. The structures can be identified in the first two or three images, but are too deformed to be identified in the last image. Keeping in mind limitations associated to viscous/dissipative effects, field line advection complemented by volume rendering of scalar quantities proves to be a valuable tool to visualize the dynamics of a turbulent flow and track the evolution of its embedded structures.

Conclusions
In this work, we presented a visualization technique for the dynamics of vector fields inspired by the physics of the hydrodynamic and MHD equations. For both cases, the existence of some form of Kelvin or Alfvén's theorems imply that field lines associated with a certain field are (in the ideal case) frozen to the velocity field and advected as material curves. This fact can be used to visualize the dynamics of one vector field under the influence of a second by what we called 'field line advection': the advection of entire fields (instead of simply points) in a time-varying velocity field.
In viscous flows, the presence of dissipation introduces a limitation to the idea of pure advection because field lines cannot be uniquely identified: advecting all points in a field line at a time T will not result in a field line at a time T + 1. Based on the observed fact that dissipation in hydrodynamic and MHD flows tends to take place in localized regions, we advect only one representative point and use the advected point as a seed to compute a new field line at T + 1. Different prioritization criteria to select the representative point at time T (or a set of points at time T from which a single point is selected at time T + 1) are allowed to give flexibility to the method, although the general rule is to use a criteria such that the selected point is in a region where the effect of dissipation is small compared with the advection.
We used two datasets to illustrate the method. In the first, we considered an ideal and a viscous simulation of an MHD flow, and we focused on a quasi-2D structure: a current sheet in the vicinity of which magnetic field lines change direction rapidly. Then, we used the method to study the evolution of structures in MHD turbulence. Field line visualization for all these examples tends to suffer badly from visual cluttering. We showed how the combination of interactively picked regions of interest, interactively placed seeds for field line integration at a fixed time, and field line advection, can give a way to identify and track structures and to visualize the evolution of the flow.