EnSight's Particle Trace Integration (RK) Algorithm

What type of integration algorithm does EnSight use to compute its particle traces?

EnSight uses an RK4 integration algorithm for massless particles (and an RK45 algorithm for massed particles) with adaptive time stepping.  In essence, the particle trace algorithm follows.

Particle Trace Algorithm

  1. Start with (and find) the emitter's xyz location in an element.
  2. Determine the starting dt (time delta) based on the characteristic length of the element's size and velocity.
  3. Feed the xyz, dt, and velocity into the RK4 integrator (see RK Details and Variable Interpolation below) and approximate new dxyz, dtheta values.
  4. Check on adaptive Dt Step-size Conditions (see below) and/or Exit Conditions (see below), and either:
  • Accept and apply the approximated values to the next iteration point (i.e. xyz[i+1]=xyz[i]+dxyz[i], time[i+1]=time[i]+dt[i], theta[i+1]=theta[i]+dtheta[i]), and then repeat steps 3-4 using with the new iteration point xyz[i+1] and either the same or an adapted dt; or
  • Reject the approximated values, adjust dt, and repeat steps 3-4 using the last iteration point xyz[i] and the adjusted dt.

RK Details
The streamline algorithm essentially starts integration using an Runge-Kutta 4th order (RK4) algorithm for massless particles [i.e. see], and an RK45 algorithm for massed particles [see any Runge-Kutta-Fehlberg Method). Essentially, the RK algorithms are based on an initial dt that is computed from the characteristic length of the element and the average velocity over the element where the emitter point for a trace is located.

The RK algoriths basically solve V=dX/dt for X[i+1] where dX=V*dt with dX=X[i+1]-X[i] which gives X[i+1]=X[i]+V*dt (see Variable Interpolations below for how we interpolate the velocity).  Substituting the RK4(X,V) function for V yields: X[i+1]=X[i]+RK4(X,V)*dt and so on ... (see The RK4(X,V) is a 4th order blend using 4 Euler steps and a the same dt over those 4 Euler steps. Note: for streamlines, 'dt' is a measure of the residence time; whereas, for pathlines 'dt' is the actual time increment.

RK Performance Consideration
Performance and accuracy of this type of algorithm is really based on mesh resolution. The coarser the mesh, the more chance for the integration to skip over local features. Some of our step-size adaptations depend on the angle of the trace with its previous segment, as well as the angle of the twist of the segment with its previous segment. For a more refined solution, it might be necessary to increase the resolution of the mesh spatially, and if transient - sometimes temporally as well to resolve localized flow features, i.e. sinks, sources, &c.

Variable Interpolation (i.e. velocity, &c.)
Within elements:
For variables based at the nodes EnSight uses element primitive (i.e. bar, triangle, quad, tet, pyramid, prism, and hex) shape, or blending, functions to interpolate variable values (i.e. velocity, &c.) at each end point of a trace segment within an element.

For variables based at the elements EnSight assumes the value of the variable is constant throughout the element.

Between time samples: (i.e. pathlines)
EnSight linearly interpolate the values between transient samples. For example, in the situation where the integrated time is found between two time steps, we first

  • Bracket the integrated time value within the closest min and max time steps.
  • Evaluate the location of the trace with both time steps (getting the appropriate shape-function values at the nodes of the element where-in the trace point lie).
  • Linearly interpolate the values between these time limits based on their time ratio contributions for each time step.

Dt Step-size Conditions
Dt step conditions are modified based on either the trace segment's rate of rotation (twist), rate of angle change (flatness and/or more curvature), maximum number of segments in an element, or below Exit Conditions.

Exit Conditions
The trace stops when it encounters either: zero velocity, maximum (residence) time, maximum number of segments, or exits the fluid domain.

Was this article helpful?
0 out of 0 found this helpful
Have more questions? Submit a request


  • Avatar
    Pauline ROCA

    Dear Mel Spencer,

    Thanks a lot for these details. Could you be more precise on how the "starting dt" computation ?

    In case of the computation of streamlines inside a vascular segmentation, i.e. a volume representing a vessel with a given diameter D within an estimated velocity vector field defined on a 3D grid at a given spatial resolution (sx,sy,sz), what would be the formula combining "the characteristic length of the element's size and velocity" ?


    Thank you in advance for your help,

    Best regards,


    Pauline Roca

  • Avatar
    Mel Spencer

    dt = || [ MAX(Cell(X)) - MIN(Cell(X)) ] * Vel(X) || / || Vel(X) || / 2


    X = <x,y,z> tupple,

    || Vel(X) || = SQRT( Vel(x)*Vel(x) + Vel(y)-Vel(y) + Vel(z)*Vel(z) ).

    Cell(X) = the <x,y,z> node set defining a cell

    || [ MAX(Cell(X)) - MIN(Cell(X)) ] * Vel(X) || =

      SQRT( ((MAX(Cell(x))-MIN(Cell(x))*Vel(x))^2 + ((MAX(Cell(y))-MIN(Cell(y))*Vel(y))^2 + ((MAX(Cell(z))-MIN(Cell(z))*Vel(z))^2 )

  • Avatar
    Pauline ROCA

    Thanks a lot for your answer. I have some questions more. Thank you in advance for your help :

    1) Could you give any reference or assumptions used to set this starting dt ?

    2) In the exemple I mentioned above :

        a) could you precise what are "Cell(X)", and "MAX(Cell(X))" please ?

        b) could you mention the interpolation method used ? (interpolation between the elements, in case on a 3D grid)

    Thanks a lot in advance,

    Best regards


  • Avatar
    Mel Spencer

    1) Could you give any reference or assumptions used to set this starting dt ?

    Starting dt criteria is based on establishing an initial dt that is a sized according to the characteristic length of the cell, the average velocity in the cell, and half that distance (allows 2 steps).

    2a) could you precise what are "Cell(X)", and "MAX(Cell(X))" please ?

    As stated above, Cell(X) = the <x,y,z> node set defining a cell, i.e. means all the coordinates of all the nodes representing the cell

    MAX(Cell(X)) = <MAX(Cell(x)), MAX(Cell(y)), MAX(Cell(z)) > = vector of the MAX <x,y,z> of the cell.  Ditto for MIN.

    2b) could you mention the interpolation method used ? (interpolation between the elements, in case on a 3D grid)

    Standard FE interpolation functions.  (see attached file)