PYTHON

velocity field gradient

def Lambda2(macroscopic_velocities, negative_cutoff = 0):
  """Use vortex detection algorithm "Lambda2" to get which nodes in your velocity field are inside a vortex.
  """
  # First define the gradients of the velocity field: Δu
  _gradients = np.gradient(macroscopic_velocities, 
                           axis = (0, 1, 2))
  gradients = np.einsum('dD...-> ...dD', 
                        np.array(_gradients))

  # Get the strain rate tensor: S = (Δu + Δuᵀ)/2
  transposed_gradients = np.einsum('...dD -> ...Dd', 
                                   gradients)
  strain_rate = (gradients + transposed_gradients)/2

  # Get the vorticity tensor: Ω = (Δu - Δuᵀ)/2
  vorticity = (gradients - transposed_gradients)/2 

  # Get the eigenvalues (λ) of S² + Ω². Use np.eigh to get ordered eigenvalues. 
  S2_O2 = strain_rate**2 + vorticity**2
  eigen_values, errors = np.linalg.eigh(S2_O2)

  # Select the locations where λ₂ < the negative cutoff. 
  # Negative cutoff must be 0 or less.
  in_vortex = (eigen_values < negative_cutoff)[..., 1]
  return in_vortex