Programmable filter - eigenvalues/lambda2

[flotus1]

I followed this guide in an attempt to get Eigenvalues / Lambda2 calculated by ParaView:

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. …

My script looks like this:

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

input=inputs[0]
Velocity_Vector = input.PointData['velocity']
STRAIN = strain(Velocity_Vector)
V_gradient = gradient(Velocity_Vector)
AAA=STRAIN**2+(V_gradient - STRAIN)**2
print(shape(AAA))
Lambda_2 = np.linalg.eigvals(AAA[:])
Lambda_2 = np.real(Lambda_2)
output.PointData.append(Lambda_2, 'Lambda_2')

Which seems to work up to the line “Lambda_2 = np.linalg.eigvals(AAA[:])”
Here I get the error “numpy.linalg.LinAlgError: 0-dimensional array given. Array must be at least two-dimensional”
The shape of AAA is printed as (12058557,3,3)
The input for this filter is point data with a velocity field in vector form.
PV version 5.11.0

I am new to programmable filters and Python in general. A nudge into the right direction would be greatly appreciated.

[flotus1]

Btw, in case this is an issue of me not working with “numpy arrays”:
I tried various methods with code snippets I found on the internet, to retain the “Velocity_Vector = algs.make_vector(Vx,Vy,Vz)” line from the original. Maybe that makes a difference compared to just importing the velocity vector as-is. So far without success.

For example, I added some calculators earlier in the chain of filters to have the velocity components available as scalars vx, vy, vz.
Changing my script to
input=inputs[0]

Vx = input.PointData["vx"]
Vy = input.PointData["vy"]
Vz = input.PointData["vz"]
print(shape(Vx))
Velocity_Vector = algs.make_vector(Vx,Vy,Vz)

I get the shape of Vx printed as (12058557,)
But the last line throws the error “AttributeError: ‘VTKCompositeDataArray’ object has no attribute ‘shape’”
???

[Francois_Mazen (François Mazen (Kitware))]

Hi @flotus1

have you tried without the copy operation [:], e.g.

Lambda_2 = np.linalg.eigvals(AAA)

Not sure to understand why the copy is needed here.

By the way you may want to read the numpy/VTK doc here: 6. Using NumPy for processing data — ParaView Documentation 5.11.0 documentation

Best,

[flotus1]

Thanks for the reply.
Yes, I tried without “[:]”, among a few other things. To no avail, the error remains the same.
I will try to understand the documentation you linked when I have a few hours of downtime.

I was hoping for some kind of a shortcut here, that doesn’t involve me learning python/numpy.
After all, this should not be horribly difficult to pull off. Yet there are plenty of instances of people asking to do the same thing online. Most of these discussions seem to go nowhere. Some of the guides you find flat out don’t work. And I found a plugin that is no longer maintained for current versions of PV: VCG - Vortex Criteria - ParaView Plugin - Visual Computing Group - Heidelberg University

Worst of all: ParaView has computation of lambda2 built-in, in the “Vortex Cores” filter. But that is not particularly useful, since the output of that filter is lines.

I also tried

...
print(isinstance(AAA, np.ndarray))
np.asarray(AAA)
print(isinstance(AAA, np.ndarray))

Which yields “False” both times

[flotus1]

I can’t seem to make this work. None of the commands I find in the documentation, or in similar threads on the internet, seem to help.

What I did now is this: I downloaded the “official” ParaView testing data files. And opened the “timestep_0_15.vts” file in PV. After applying the “Clean to grid” filter, it is shown as an “Unstructured Grid” data type. Lo and behold, using the programmable filter from my initial post, adjusted to the “Momentum” vector present in this data, works without a hitch.

After importing the data I actually want to work with, PV shows the data type as “Multi-Block Dataset”.
It is originally cell-centered data in an Ensight Gold file format. So far, I used the “Cell data to point data” filter in order to work with point data. Applying the “Clean to grid” filter here does not change the data type. It is still “Multi-Block Dataset”.

So I guess this is my new question:

1. what do I need to change in my programmable filter, in order to work with Multi-Block Datasets
2. -or- which filters do I need to apply to a Multi-Block Dataset, for PV to change the data type to unstructured grid. Or any other data type where the usual way of writing programmable filters works. EDIT: the “Merge Blocks” filter seems to do the trick

[mwestphal (Mathieu Westphal (Kitware))]

@flotus1

which filters do I need to apply to a Multi-Block Dataset, for PV to change the data type to unstructured grid. Or any other data type where the usual way of writing programmable filters works.

MergeBlocks.

[flotus1]

Yeah, just after I finished typing, I found that myself.

[flotus1]

With that out of the way, I tried comparing the results of my filter to the ones obtained with this plugin: VCG - Vortex Criteria - ParaView Plugin - Visual Computing Group - Heidelberg University
Had to go back to PV version 5.6 to use it though.

The results don’t look too promising. With Lambda2 as a criterion to identify vortices in turbulent flows, the results on the left (from the plugin) are much more plausible.
Bottom row of images is an iso-surface with a value of -5e9.

[bastian]

@flotus1

My script looks like this:

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

input=inputs[0]
Velocity_Vector = input.PointData['velocity']
STRAIN = strain(Velocity_Vector)
V_gradient = gradient(Velocity_Vector)
AAA=STRAIN**2+(V_gradient - STRAIN)**2
print(shape(AAA))
Lambda_2 = np.linalg.eigvals(AAA[:])
Lambda_2 = np.real(Lambda_2)
output.PointData.append(Lambda_2, 'Lambda_2')

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.

I wrote a filter that should give you 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')