This article discusses the vtk stream line filter and how it is implemented.
General idea of the streamline algorithm
Assuming there is vector field in space, and we want to visualize it, we might put a point (fluid element) in it (seed), then that point will move among the space because of the existance of the vector field (The flow field is steady here, and the data is generated by one timestep). We might trace the trajectory of that particle to see how it moves with the increase of the timestep. This line is called the streamline.
The formal defination of streamline: refer to this video: https://www.youtube.com/watch?v=8MUPQbazkLQ.
Streamlines are lines that are tangent to the velocity field.
More formal and detailed description can be found here. The caculation of the streamline is the integration of the vector field based with the increase of the time step (this can be a small segment). If we combines multiple segments together, we can get a final streamline which across multiple step.
Other relates concepts: vector field reconstruction. Refer to this work (Vector field reconstruction from sparse samples with applications). The input is a set of sample point with the vector, the output is an image with different color that shows the vector field in the whole domain.
There are some realted works that use the ML techniques to do this (https://www.youtube.com/watch?v=sz_MbjftQTQ), such as SVM or deep learning things. Given the sparse set of particles in the vector field, the goal is to create a function that approximate the vector field over the domain When we say the reconstruction, we want to reconstruct that function from the sample particles. The video mentioned above provides a good explanation.
The streamline example in paraview
There are two examples on vtk-example about the streamline, one is this (https://kitware.github.io/vtk-examples/site/Python/Visualization/StreamLines/) another is this (https://kitware.github.io/vtk-examples/site/Cxx/Visualization/Kitchen/).
The associated data can be found here (https://github.com/lorensen/VTKExamples/tree/master/src/Testing/Data/)
We first show how to load these data into the paraview and show the results of using the stream line in it.
The first example is the combxyz.bin and the combq.bin case. Thi data format is called the PLOT3D. In the Paraview, we first load the combxyz.xyz using the PLOT3D file fomat, and then click the Q file name at the properties section to load all the data. (A typical plot3d data contains the q file and the xyz file). Then the file can be loaded successfully. For the kitchen.vtk example, the data is in vtk format, and we just need to load it direactly into the vtk.
After that, just click the streamline trace filter and we can know the streamline. At the information section of the paraview, we can see the type of the streamline is the polygonal mesh. The color of the streamline in paraview can also change with when we use different variables.
By clicking the gear button at the properties section, we can show all the settings of the associated filter. We can adjust the control the seed, the number of steps and the granularity of step (length of step) etc.
This is another good tutorial (https://www.youtube.com/watch?v=1jml3j3BmUg) about how to use the streamline for paraview, we can use the tube filter and glyph graph to show the direaction of the streamline.
Using the stream line in VTK code
This is one example code about the streamline based on VTK in cpp code
It use the plane as a source, then setting necessary variables such as propagation steps and the integration step. After that, it writes the resutls into the disk
This is the another example that use the kitchen.vtk data (https://github.com/wangzhezhe/5MCST/blob/master/vtk_example/streamline/streamline_exampleKitchen.cpp) The seed source is a line in this case.
We can load these generated data into the paraview and do the further analysis. Then we dive into more details to see how the stream trace filter is implemented in VTK.
Understanding the how VTK implement the stream line filter
It is not easy to understand all the details about how the VTK implement the stream tracer. We only take some key steps here.
The assocaited filter vtkStreamTracer.h or vtkStreamTracer.cxx is located at VTK/Filters/FlowPaths in the vtk repo, there is also vtkGenericStreamTracer in the VTK/Generic repo and it looks that the code are similar.
The start point of the vtkStreamTracer filter is vtkStreamTracer::Integrate after all kinds of data pre processing. The core operation is integration operation use the Runge-Kutta method. integrator->ComputeNextStep. The point2 is the next step, then update the point1 and use the new computed point2 to update the code gradually. The readme of that function can show lots of details, here are some common parameters: seeds is the start point, we could assume we put an particle at this position, the steps indicates how many step we want to go, and then we move that imaginary particle along the current field to see how it moves. That is the key idead of the stream line. We can specify more specific parameters such as the step size or the granularity of the step. In paraview, click the gear button at the left side, it shows all parameters and their default settings used by associated filter.
Some idea about the runge-kutta method
These are several useful materials:
https://www.youtube.com/watch?v=IYyECJmxfG0
The code example that solves the real problem
https://www.youtube.com/watch?v=1FYrnwqWQNY
This is a good theoretical explanation:
https://lpsa.swarthmore.edu/NumInt/NumIntSecond.html
On the online tutorial such as wiki page, they always provide a basic form which is the 1d situation. The actual code ususlly focuse on 2d or 3d cases. Instead of using x, y, z … The actucal code ususlly use a vector of x to represents the input indepedent variables. If we understanding naively for runge-kutta method, it just tries to estimate a proper k at the tangential direction to compute the next point value. The higher order method has a more accurate results.
Another key point is that, in the tutorial, we always have a concrete function at the right side as an example, we know the function form. But for the numerical case, such as streamline case, we do not have that function and we just have a set of points and each points has a 3d coordinates and associated scalar or vector values. In the VTK, there is FunctionValues associated with each function, it compute the value of x for the function f. x and f are two parameters of the FunctionValues. In the streamline case, it uses the vtkInterpolatedVelocityField, which can compute the function value based on the data sets. According to the VTK text book, the change rate of stream line trajactory within specific delt t is the vector field.
For example, VTK/Common/Math/vtkRuneKutta2.h provides necessary function for the RungeKutta2.
xprev is the previous data
xnext is the next step data
t is the current time
delT is the time slice
vtkFunctionSet specifies an abstract interface for set of functions
of the form F_i = F_i(x_j) where F (with i=1..m) are the functions
and x (with j=1..n) are the independent variables.
The this->Derivs is updated gradually, it first represents the k1 then represents the k2. FunctionValues compute the value at the next point (t+h/2)
Some notes about the flow filter in vtkm
For the particle advection and the streamline filter, the FilterParticleAdvectionSteadyState is adopted to execute actual integration task. for the PathParticle and Pathline filter, the FilterParticleAdvectionUnSteadyState is adopted to execute the actual integration operation.
Following the typical way of the vtkm filter, the start position of the flow filter is flow/FilterParticleAdvectionSteadyState.cxx or flow/FilterParticleAdvectionUnSteadyState.cxx (v2.0), one important step here is to transfer the input dataset into the integrator, during the execution of the filter, what we manupleated is integrator instead of the raw data blocks. This is entity associated with data block for executing integration opertation. Which is created at the beginning of the filter by function CreateDataSetIntegrator
The particle adopted by flow visulization is inherited from the general particle under vtkm/Particle.h. The self defined particle used for flow visualization is under flow/worklet/Particles.h file. If we search SetTerminate under the filter/flow folder, we can get all necessary information such as when particles are actually terminated.
Some optimizations
In VTK-m it adopts multiple strategies to optimize the particle advection operations.
One idea is called the block duplication (or prefetching), the idea is simple, assuming there are two ranks and two blocks, instead of putting one block at each rank, there are two blocks for each rank. So when particles fly through one block (such as block0) to another block(such as block1), it can decide to send particle to block 1 either in rank0 or in rank1. It is potential to save some communication overhead. and make work more balanced distributed among all ranks (avoid the communication overhead). If there is some monitoring service could help to prefetch associated block to the same ranke, these communication overhead can be saved.
The view of particles, blocks, ranks. Depedning of different scenarios, store information on ranks or each particles or blocks to improve the performance. Or distributed the associated among different resources to improve the overhead (such as parallel over seeds or parallel over data blcoks)
The mixture of RK4 method and Euler method for small steps when there is boundry acrossing. RK4 method is accurate compared with Euler method. When a particle flies close to the boundry, the RK4 method might fail when the actual distance towards the boundry is less than current step size. Then we can dynamically change the step size, such as using the binday search to dynamiclly shrink the step size to let the particle advected a little bit more. By this way, we can let particle flies close to the boundry as possible. Then we the particle flies out of the current boundry to the next data block, we can then use Euler method once for the convenience of computation. If we do not use extra small steps, one step based on Euler method will lead the particle to fly out of the boundry which can introduce the inaccuracy.
Another optimization is to terminate the particle when it is trapped into a particular block, which is mentioned in the “A Study of Parallel Particle Tracing for Steady-State and Time-Varying Flow Fields”. Assuming the syncrounous version of the MPI is adopted, there are explicit splition of round during the advection process. After one round, particle in the same block will be terminated rather than advancing to the next round.
Key matadata for streamline
BoundsMap let’s use vtkm’s bounds map as an example. There are two interface, the first interface treat each block as different one. Another interface used for bounds map where there are duplication. User can set the specific block id for these blocks. The idea of bounds map is to map the block id into the id of the rank.
References
https://kitware.github.io/vtk-examples/site/VTKBook/06Chapter6/#63-vector-algorithms
pipeline of the scientific vis
https://medium.com/researchsummer/visualization-of-a-vector-field-9402615c780a
