I recently found myself in the situation where I wanted to get the trajectories of streamlines starting from a specific point. Lots of visualization software offer the possibility to plot streamlines. Two of these are Paraview and Tecplot. Paraview is open-source and generally good, although on dense meshes, streamlines tend to stop randomly. Tecplot is commercial software and I never had any problems plotting streamlines with it. However, with neither do I know how to extract streamline trajectories. If you know how to, please let me know!

In the meantime, some quick code to get streamline trajectories in Python. I am using results calculated with oomph-lib, see this earlier tutorial. The output format is optimized for Tecplot and in this example, contains 5 values: x-position, y-position, x-velocity, y-velocity and pressure. The following function will read those values from a file:

from __future__ import absolute_import, division, print_function import numpy as np def _read_file(path): '''Read output file with 5 values: x, y, u, v, p ''' file_handle = open(path, 'r') #open file lines_list = file_handle.readlines() #read all lines x=[]; y=[]; u=[]; v=[]; p=[]; for line in lines_list: #go through each line ent = line.split() # split into separate entries #ignore lines with text and those that have a different number of entries #the tecplot file contains zone information and information on the #connectivity of each element if len(ent) == 5 and ent[0][0] != 'Z': x.append(float(ent[0])) y.append(float(ent[1])) u.append(float(ent[2])) v.append(float(ent[3])) p.append(float(ent[4])) return np.array(x), np.array(y), np.array(u), np.array(v), np.array(p)

The data is unstructured. This causes a bit of complication. To use the in-build python functions to e.g. plot the flow field, regularly spaced data is needed. The function scipy.interpolate.griddata thankfully does that for us.

</pre> from __future__ import absolute_import, division, print_function import numpy as np from scipy.interpolate import griddata def interpolate_data_to_uniform_grid(data, n_pts=100, plot=False): ''' Approximates unstructured data with a uniform grid Input: data - a 2D numpy array with 5 columns: x, y, u, v, p n_pts - int: number of points in the uniform grid. Going above 1'000 quickly eats GB's of RAM... ToDo: accept tuple for x and y dimension plot - bool: show results, mainly used for debugging Returns: The data interpolated on a uniform grid: xi, yi, ui, vi, pi. Note that xi and yi are arrays not necessarily of equal length and a value for every combination of xi and yi is returned in ui, vi and pi. ''' x = data[0,:]; y = data[1,:]; u = data[2,:]; v = data[3,:]; p = data[4,:]; # define grid. xi = np.linspace(np.min(x),np.max(x),n_pts*2) yi = np.linspace(np.min(y),np.max(y),n_pts) # grid the data. ui = griddata((x, y), u, (xi[None,:], yi[:,None]), method='cubic') vi = griddata((x, y), v, (xi[None,:], yi[:,None]), method='cubic') pi = griddata((x, y), p, (xi[None,:], yi[:,None]), method='cubic') # griddata fills in NAN if it cannot fit, as is the case when outside the # area of the mesh. For our purposes, zero velocity there should ensure that # nothing every goes there ui = np.nan_to_num(ui) vi = np.nan_to_num(vi) pi = np.nan_to_num(pi) # contour the gridded data, plotting dots at the randomly spaced data points. if plot: fig = plt.figure(figsize=(8, 6));ax = fig.add_subplot(111) plt.quiver(x,y,u,v,zorder=1) plt.axes().set_aspect('equal', 'datalim') plt.title("Original Data") fig = plt.figure(figsize=(8, 6));ax = fig.add_subplot(111) plt.quiver(xi, yi, ui, vi,zorder=1) plt.title("Gridded Data") plt.axes().set_aspect('equal', 'datalim') return np.array(xi), np.array(yi), np.array(ui), np.array(vi), np.array(pi)

Now that we have uniformly gridded data, it is easy to interpolate the velocities, using scipy.interpolate.interp2d. We can then simply call the interpolation to get the velocity at any point in the domain. The streamlines are calculated using a simple Euler method $x_{n+1} = x_n + u(x_n) * \delta_t$, where $x_n$ is the position at the the $n$th time-step, $u(x_n)$ the velocity at the position $x_n$ and $\delta_t$ the time-step.

from __future__ import absolute_import, division, print_function import numpy as np from scipy import interpolate def get_single_streamline(xi, yi, ui, vi, start_pos, DELTAT_T=1e-3, MAX_COUNTER=1e6): '''Get trajectory of a streamline starting at start_pos Input: xi, yi, ui, vi - Position and velocity on a uniform grid start_pos - Tuple, x and y position of starting point DELTAT_T=1e-3 - Time step used MAX_COUNTER=1e10 Max number of steps Returns: traj_x, traj_y - x and y position of streamline ''' #normalize speed, so that the time-step means the same thing max_speed = np.max(np.sqrt(ui*ui + vi*vi)) ui=ui/max_speed; vi=vi/max_speed #Interpolate the grid #This is somewhat doggy as we are interpolating from interpolated data. #This allows us to get the value of the velocities at any point. # A more elegant way would be to interpolate the unstructured data but this works okay. fui = interpolate.interp2d(xi, yi, ui, kind='cubic') fvi = interpolate.interp2d(xi, yi, vi, kind='cubic') #Define the dimension of our data MAXX=np.max(xi);MINX=np.min(xi) MAXY=np.max(yi);MINY=np.min(yi) # Now use Euler to predict the next position based on the velocity traj_x=[start_pos[0]] traj_y=[start_pos[1]] counter=0; while(counter < MAX_COUNTER): # Stop when we reach the edge of the domain if (traj_x[-1] > MAXX or traj_x[-1] < MINX or traj_y[-1] > MAXY or traj_y[-1] < MINY): break; #Output to ensure the user that we are still doing something if counter%(MAX_COUNTER/10.0) ==0: print("%dth iteration" %counter) #find next position by simple euler xn = traj_x[-1] + fui(traj_x[-1], traj_y[-1]) * DELTAT_T yn = traj_y[-1] + fvi(traj_x[-1], traj_y[-1]) * DELTAT_T traj_x.append(xn); traj_y.append(yn) counter +=1 return np.array(traj_x), np.array(traj_y)

All we need to do now is call this function to get the streamlines. An (messy) example is provided in this folder for the case considered in this tutorial. I added some functions to plot the boundaries of the domain. My aim was to get the relationship of the the $y$ position above the obstacle and the resulting angle with respect to the centreline in the expansion. In the graphs below you can see that the streamlines plotted with Tecplot show the same qualitative behaviour as those plotted with this python code.