[gaolbert]

Dear All,

I have been working on the following problem for a while now, and neither search machines nor any of the other threads I found provided me with a satisfying answer:

Based on the Paper of Haimes&Sujudi (1995) I am trying to implement a routine that provides the position of vortex cores in turbulent flows with free vortices:
There are three main steps:
1. In the vicinity of vortex cores, the gradient tensor of the velocity field will have 1 real and two complex conjugate eigenvalues.
2. For each cell, the eigenvector corresponding to the real eigenvalue is evaluated and the projection of the actual velocity onto the direction of the real eigenvector (n=normalized eigenvector) is subtracted to obtain a 'reduced' velocity w:
w=u - ( u*n)*n. This velocity will be exactly 0 in the vortex core.
3. Using linear interpolation, the places in which w has a value of 0 are evaluated on the cell faces. In the original implementation (in pV3) the field values (e.g. vorticity) were then interpolated to those positions and exported. For my purpose, it would suffice to simply set a volScalarField value v[i]=1 if a value of w[i]=0 (or similarly small) is detected for the corresponding cell and leave it at v[i]=0 otherwise. Using e.g. paraview for postProcessing, I could then write out the position and cell values of all cells with a value of v=1.

So far so good. I checked the routines I have written and they seem to work for non-complex eigenvalues. But as soon as I try with the complex sort, I get the following error (after a considerable amount of warnings due to complex eigenvalues:

--> FOAM Warning : 
    From function eigenValues(const tensor&)
    in file primitives/Tensor/tensor/tensor.C at line 157
     complex eigenvalues detected for tensor: (1011.6531 196.83104 189.08022  -33987.54 -6265.4809 -6080.2742 29041.968 5457.2612 5242.8922)
--> FOAM Warning : 
    From function eigenValues(const tensor&)
    in file primitives/Tensor/tensor/tensor.C at line 157
     complex eigenvalues detected for tensor: (-158.45982 -5.1601569  -16.312047 -34371.929 -4067.9722 -5129.8848 28292.557 3404.9461  4233.1667)
etc.....
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib64/libc.so.6"
#3  Foam::eigenValues(Foam::Tensor<double> const&) at ??:?
#4   Foam::eigenValues(Foam::Field<Foam::Vector<double> >&,  Foam::UList<Foam::Tensor<double> > const&) at ??:?
#5  ? at ??:?
#6  ? at ??:?
#7  __libc_start_main in "/lib64/libc.so.6"
#8  ? at /home/abuild/rpmbuild/BUILD/glibc-2.22/csu/../sysdeps/x86_64/start.S:121

Here is the code I use. Sorry for the ugly workaround using the for-loop. I read that the error might be due to a zero division and tried to make sure every value inside the tensor was unequal to zero. I already tried setting VSMALLTens up using createFields.H, but wasn't sure how to assign default values for a volTensorField. If anyone has any recommendations for improvement, please shoot:

const volTensorField gradU(fvc::grad(U));
    volTensorField VSMALLTens = gradU;
    forAll(VSMALLTens, i)
    {

        VSMALLTens[i]=ROOTVSMALL;

    }

    volScalarField EigValsOfgradU = eigenValues(gradU+VSMALLTens)().component(vector::Y);

I only use the second eigenvalue as OF stores the real eigenvalue in this position when complex eigenvalues are detected.
Do you have any recommendations on how I may avoid the Floating point exception? Thanks in advance!

- using OF 1606+