Eigenvalue computation

Hello all,

I am implementing Lambda-2 criterion in Paraview 5.6.0 mounted on Linux. What it does is to provide a view of the presence of a vortex in a 3D flow by knowing when at least two of the eigenvalues of (||S||^2 + ||Omega||^2) are negative. Here, S and Omega are the symmetric and antisymmetric parts of the Velocity gradient tensor.

By following the instructions of some other posts where they dealt with this implementation, I have made the following:

  • Gradient of unstructured data set. Inputting Velocity as Scalar Array and Compute Gradients
  • Python calculator. Expression “eigenvalue(strain(Velocity)**2 + (Gradients - strain(Velocity))**2)”. Array name “Lambda”. This computes the three eigenvalues of the expression ||S||^2 + ||Omega||^2
  • Calculator. “Lambda_Y”. Result Array Name “Lambda2”. Get the second eigenvalue of the previous result. Since the three eigenvalues are ordered such that Lambda1 >= Lambda2 >= Lambda3, if Lambda2 is negative, Lambda3 is negative too, so a vortex is found.

1- I am having some issues since the results that I obtain are not the expected ones. So I took one step back and I have been looking terms by terms what is happening there. I have realized that if I compute the Gradient of Unstructured Data Set to get the velocity gradient and I compute my means of a Python calculator “strain(Velocity)”, the tensor I am obtaining in both cases does not share the trace. Both of them should share the trace since strain(Velocity) = 1/2*(Gradient(Velocity) + transpose(Gradient(Velocity)), and it is not the case.

2- Also, I am evaluating if the issues comes from the eigenvalue calculation itself. To do so, I tried to simple compute the eigenvalues of the velocity gradient by means of the function eigenvalue(Gradients) as well as by means of numpy.linalg.eigvals(Gradients). Think this second version to compute should work better, however, it draws the following error (it does not recognize the Gradients tensor as a tensor):
Traceback (most recent call last):
File “”, line 4, in
File “/home/softs/softs/ParaView-5.6.0-MPI-Linux-64bit/lib/python2.7/site-packages/paraview/calculator.py”, line 192, in execute
retVal = compute(inputs, expression, ns=variables)
File “/home/softs/softs/ParaView-5.6.0-MPI-Linux-64bit/lib/python2.7/site-packages/paraview/calculator.py”, line 140, in compute
retVal = eval(expression, globals(), mylocals)
File “”, line 1, in
File “/home/softs/softs/ParaView-5.6.0-MPI-Linux-64bit/lib/python2.7/site-packages/numpy/linalg/linalg.py”, line 1021, in eigvals
_assertRankAtLeast2(a)
File “/home/softs/softs/ParaView-5.6.0-MPI-Linux-64bit/lib/python2.7/site-packages/numpy/linalg/linalg.py”, line 204, in _assertRankAtLeast2
‘at least two-dimensional’ % a.ndim)
numpy.linalg.linalg.LinAlgError: 0-dimensional array given. Array must be at least two-dimensional

Thank you in advance,

Andrea Arroyo

Hello,

Here you have a code snippet to compute Lambda 2 in a Python Programmable Filter:

import numpy as np
from vtk.numpy_interface import algorithms as algs

INPUT=inputs[0]
Vx = INPUT.PointData[“X_VELOCITY”]
Vy = INPUT.PointData[“Y_VELOCITY”]
Vz = INPUT.PointData[“Z_VELOCITY”]

Velocity_Vector = algs.make_vector(Vx,Vy,Vz)

STRAIN = strain(Velocity_Vector)
V_gradient = gradient(Velocity_Vector)
AAA=STRAIN**2+(V_gradient - STRAIN)**2
Lambda_2 = np.linalg.eigvals(AAA[:])
Lambda_2 = np.real(Lambda_2)
output.PointData.append(Lambda_2, “Lambda_2”)

Regards

Miguel

3 Likes

Hello,

I’d like to calculate the principal strains so I modified this code into

import numpy as np
from vtk.numpy_interface import algorithms as algs

INPUT=inputs[0]
Dx = INPUT.PointData["MS_Displacement_X"]
Dy = INPUT.PointData["MS_Displacement_Y"]
Dz = INPUT.PointData["MS_Displacement_Z"]

Displacement_Vector = algs.make_vector(Dx,Dy,Dz)

STRAIN = strain(Displacement_Vector)
principal_strain = np.linalg.eigvals(STRAIN[:])

outarray0.Arrays[ii] = principal_strain[:,0]
outarray1.Arrays[ii] = principal_strain[:,1]
outarray2.Arrays[ii] = principal_strain[:,2]

output.PointData.append(outarray0, "principal_strain_0")
output.PointData.append(outarray0, "principal_strain_1")
output.PointData.append(outarray0, "principal_strain_2")

However, there’s an error shows:

AttributeError: 'VTKNoneArray' object has no attribute 'shape'
Traceback (most recent call last):
  File "<string>", line 22, in <module>
  File "<string>", line 10, in RequestData
  File "C:\Program Files\ParaView 5.10.1-MPI-Windows-Python3.9-msvc2017-AMD64\bin\Lib\site-packages\vtkmodules\numpy_interface\algorithms.py", line 861, in make_vector
    return algs.make_vector(arrayx, arrayy, arrayz)
  File "C:\Program Files\ParaView 5.10.1-MPI-Windows-Python3.9-msvc2017-AMD64\bin\Lib\site-packages\vtkmodules\numpy_interface\internal_algorithms.py", line 558, in make_vector
    if len(ax.shape) != 1 or len(ay.shape) != 1 or (az is not None and len(az.shape) != 1):
AttributeError: 'VTKNoneArray' object has no attribute 'shape'

Any ideas? Thanks in advance!

Best,
Gening

I think that you have a few errors in your code:

  1. STRAIN**2 is probably an elementwise square. What you actually want is the matrix product that matches the signature (n,k),(k,m)→(n,m), which matmul should provide
  2. np.linalg.eigvals does not order the eigenvalues, but you want lambda_2 to correspond to the second largest eigenvalue.

This filter should give the proper result, but I haven’t validated it:

import numpy as np
from vtk.numpy_interface import algorithms as algs

# Compute the lambda_2 Vortex criterion, see pp.76-77 in
# J. Jeong and F. Hussain, “On the identification of a vortex,” 
#  Journal of Fluid Mechanics, vol. 285, pp. 69–94, 1995, doi: 10.1017/S0022112095000462.

vvector = inputs[0].PointData['velocity']

vstrain = strain(vvector)
vskew = gradient(vvector) - vstrain

aaa = matmul(vstrain, vstrain) + matmul(vskew, vskew)

# since aaa is symmetric, it only has real eigenvalues
lambdas = eigenvalue(aaa)
lambdas = real(lambdas)
lambda2 = sort(lambdas)[:,1]

output.PointData.append(lambda2, 'lambda2')