[Svensen]

Hi, I'm trying to implement a precise viscosity model from experimental measurements.

I've created a folder "Spline" and two files "Spline.C" and "Spline.H".
I've redefined calcNu() method in the following was:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
    volScalarField strain = strainRate();
     volScalarField nu = strain;
     
//for each element of a field
     for (unsigned i=0; i<strain.size(); i++){
         double n = strain[i]; //get a strain-rate
    
        const double dn = 250/sizeof(visc);
        n /= dn; //scale it to match an index of array
    
/*visc array includes 2000 values of measured viscosity*/
        unsigned x1 = floor(n);
        if (x1 < 1) x1 = 0; //check minimum value
        if (x1 >= sizeof(visc)) x1 = sizeof(visc)-2; //check maximum value
        
        unsigned x2 = x1 + 1;
    
        if (x1 != x2){
            double y1 = visc[x1];
            double y2 = visc[x2];
        
            double k = y2 - y1; //linear interpolation
            double b = x2*y1 - x1*y2;
        
            nu[i] = k*n + b;
        }
        else{
            nu[i] = visc[x1];
        }
     }
    
    return nu;
}

It compiles with no error, however after start of pimpleFoam it exits with the following log:

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::viscosityModels::Spline::calcNu() const at ??:?
#4  Foam::viscosityModels::Spline::Spline(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#5  Foam::viscosityModel::adddictionaryConstructorToTable<Foam::viscosityModels::Spline>::New(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#6  Foam::viscosityModel::New(Foam::word const&, Foam::dictionary const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#7  Foam::singlePhaseTransportModel::singlePhaseTransportModel(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#8  ? at ??:?
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  ? at ??:?

What can be wrong with my code ?

[wyldckat]

Greetings Svensen,

I've come to this thread, after looking for it due to the following bug report: http://bugs.openfoam.org/view.php?id=2316
In that bug report you have a slightly different source code:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
        volScalarField strain(strainRate());
     
        //for all elements in the field
     for (unsigned i=0; i<strain.size(); i++){
         double n = strain[i];
    
        const double dn = 250/sizeof(visc);
        n /= dn; //calc an index in array
    
        unsigned x1 = floor(n);
        if (x1 < 1) x1 = 0; //correct minimum value
        if (x1 >= sizeof(visc)) x1 = sizeof(visc)-2; //correct maximum value
        
        unsigned x2 = x1 + 1;
    
        if (x1 != x2){
            double y1 = visc[x1]; //visc array contains the measured viscosity values 
            double y2 = visc[x2];
        
            double k = y2 - y1; //because x2-x1 = 1
            double b = x2*y1 - x1*y2; //because x2-x1 = 1
        
            strain[i] = k*n + b;
        }
        else{
            strain[i] = visc[x1];
        }
     }
        
    return dimensionedScalar("one", dimLength*dimLength, 1.0)*strain;
}

Furthermore you commented that:

After some experiments I've just found that problem is in line "strain[i] = k*n + b;". However I still don't understand why the problem happens here. OpenFOAM documentation states that volScalarFiels is based on List, therefore a value can be assigned to individual element. Is it some issue in OpenFOAM sources ?

It's very possible that you've isolated a symptom, but not the cause of the problem.
If we take a look at the printed stack trace, the first tell-tale sign is this:

#1  Foam::sigFpe::sigHandler(int) at ??:?

This means that there was a https://en.wikipedia.org/wiki/SIGFPE - which quoting from there:

The SIGFPE signal is sent to a process when it executes an erroneous arithmetic operation, such as division by zero (the name "FPE", standing for floating-point exception, is a misnomer as the signal covers integer-arithmetic errors as well).

Other possible reasons is if there is a bad logarithm calculation or a division by infinite and so on.

When we take into account this signal for a bad mathematical operation and that the stack trace indicates that this occurs when the transport model is being created for the first time, I'm guessing that something wrong is going on with the initialization of the other fields:

• It is not clear what "visc" is and how it was defined.
• These two operations are the most likely suspects for a bad mathematical operation:
  

           const double dn = 250/sizeof(visc);
          n /= dn; //calc an index in array

• And as I mentioned in the bug report, I suspect that "sizeof(visc)" is simply being used incorrectly, because "sizeof" gives the number of bytes the variable needs on RAM, not the number of items on a list. Which means that if "visc" is in fact a list, then this means that any accesses with "visc[x1]" will result in illegal memory accesses.

When in doubt, use old-school debugging, i.e. use Info to output to screen the values that are being used, e.g.:

Info << "size visc: " << sizeof(visc) << endl;
Info << "dn: " << dn << endl;

Beyond this, the way that OpenFOAM goes through lists is by using forAll, e.g.:

forAll(strain, itemi)
{
   Info << strain[itemi] << endl;

}

Best regards,
Bruno

[Svensen]

Thanks, wyldcat.

The first error was with sizeof. I forgot that sizeof returns a number of bytes instead of number of elements.

However, the main problem was in the way of how "visc" array was initialized. I defined it as a private member of a class and initialized it in the body of constructor.

BUT OpenFOAM uses a constructor with initialization list, where nu_ object initializes before body of constructor is executed. That's why calcNu() method executes first and only after this my constructor will initialize a "visc" array. This produces "strange" values in "visc" array during execution of calcNu() method.

I've tried to move an initialization of nu_ object from initialization list to body of constructor, but compilation failed.
To overcome this, I've simply moved an initialization of "visc" array to the code of calcNu() method. It compiles without errors and executes first step without error, but then crashes.

I need some time to check my code to find out why it crashes.

But nevertheless, can I initialize a "visc" array outside from calcNu() method ? Because it seems weird to initialize it every time the calcNu() method executes.
I've attached the sources below.

[Svensen]

About crash...

First of all, I've started with Newtonian viscosity, i.e. I've set all elements of visc array to 4.1 mPa*s, which corresponds to kinematic viscosity of 3.9e-6 m2/s. I've checked that every element of "strain" object is initialized by this value during calcNu() method. However the simulation crashes again.

The return command looks like this:

return dimensionedScalar("one", pow(dimLength, 2), 1.0)*strain;

Error log is the following:

Courant Number mean: 0.970093 max: 303.708
Time = 0.002

PIMPLE: iteration 1
GAMG:  Solving for Ux, Initial residual = 0.931379, Final residual = 5.80411e-08, No Iterations 1000
GAMG:  Solving for Uy, Initial residual = 0.95176, Final residual = 2.56661e-08, No Iterations 1000
GAMG:  Solving for Uz, Initial residual = 0.924978, Final residual = 1.00825e-07, No Iterations 1000
GAMG:  Solving for p, Initial residual = 0.663798, Final residual = 8.97244e-08, No Iterations 48
time step continuity errors : sum local = 1.64355e-05, global = 1.05359e-05, cumulative = 1.05436e-05
GAMG:  Solving for p, Initial residual = 0.00216456, Final residual = 5.99289e-08, No Iterations 11
time step continuity errors : sum local = 0.00179684, global = -0.00101075, cumulative = -0.00100021
PIMPLE: iteration 2
[9] #0  Foam::error::printStack(Foam::Ostream&) at ??:?
[9] #1  Foam::sigFpe::sigHandler(int) at ??:?
[9] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
[9] #3  double Foam::sumProd<double>(Foam::UList<double> const&, Foam::UList<double> const&) at ??:?
[9] #4  Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[9] #5  Foam::GAMGSolver::solveCoarsestLevel(Foam::Field<double>&, Foam::Field<double> const&) const at ??:?
[9] #6  Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMatrix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
[9] #7  Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
[9] #8  ? at ??:?
[9] #9  ? at ??:?
[9] #10  ? at ??:?
[9] #11  ? at ??:?
[9] #12  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
[9] #13  ? at ??:?
[sergey-HP:109750] *** Process received signal ***
[sergey-HP:109750] Signal: Floating point exception (8)
[sergey-HP:109750] Signal code:  (-6)
[sergey-HP:109750] Failing at address: 0x3e80001acb6
[sergey-HP:109750] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x36cb0)[0x7fc9651ffcb0]
[sergey-HP:109750] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0x37)[0x7fc9651ffc37]
[sergey-HP:109750] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x36cb0)[0x7fc9651ffcb0]
[sergey-HP:109750] [ 3] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam7sumProdIdEEdRKNS_5UListIT_EES5_+0x39)[0x7fc966566689]
[sergey-HP:109750] [ 4] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam5PBiCG5solveERNS_5FieldIdEERKS2_h+0x7d5)[0x7fc9663b93d5]
[sergey-HP:109750] [ 5] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver18solveCoarsestLevelERNS_5FieldIdEERKS2_+0x207)[0x7fc9663d46a7]
[sergey-HP:109750] [ 6] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver6VcycleERKNS_7PtrListINS_9lduMatrix8smootherEEERNS_5FieldIdEERKS8_S9_S9_S9_S9_S9_RNS1_IS8_EESD_h+0xbd8)[0x7fc9663d63f8]
[sergey-HP:109750] [ 7] /home/sergey/OpenFOAM-dev/OpenFOAM-dev/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam10GAMGSolver5solveERNS_5FieldIdEERKS2_h+0x49b)[0x7fc9663d818b]
[sergey-HP:109750] [ 8] pimpleFoam[0x459a28]
[sergey-HP:109750] [ 9] pimpleFoam[0x462b1b]
[sergey-HP:109750] [10] pimpleFoam[0x462dfb]
[sergey-HP:109750] [11] pimpleFoam[0x424978]
[sergey-HP:109750] [12] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xf5)[0x7fc9651eaf45]
[sergey-HP:109750] [13] pimpleFoam[0x4253e8]
[sergey-HP:109750] *** End of error message ***
--------------------------------------------------------------------------
mpirun noticed that process rank 9 with PID 109750 on node sergey-HP exited on signal 8 (Floating point exception).

It seems to me that something wrong happens when I return the "strain" object from calcNu() method. Maybe some pointer violation or type mismatch...

[Zeppo]

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
...
    volScalarField strain(strainRate());
...
    return dimensionedScalar("one", dimLength*dimLength, 1.0)*strain;
}

->

Foam::viscosityModels::Spline::Spline
(
    const word& name,
    const dictionary& viscosityProperties,
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    ...
    nu_(calcNu())
{}

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
    const volScalarField& strainRate = strainRate()();
    tmp<volScalarField> tNu
    (
        new volScalarField
        (
            IOobject
            (
                "nu",
                strainRate.time().timeName(),
                strainRate.mesh(),
                strainRate.readOpt(),
                strainRate.writeOpt()
            ),
            strainRate.mesh(),
            sqr(dimLength)*strainRate.dimensions()
            strainRate.boundaryField().types()
        )
    );

    volScalarField& nu = tNu.ref();
    //volScalarField& nu = tNu(); //- for OF < 4.0
    forAll(nu, i)
    {
        ...
        nu[i] = ....
        ....
    }

    //- optionally you can also alter your boundary values:
    volScalarField::GeometricBoundaryField& bNu = nu.boundaryField()
    //volScalarField::Boundary& bNu = nu.boundaryFieldRef() //- for OF < 4.0
    forAll(bNu, i)
    {
        scalarField& pNu = bNu[i]; // patch field
        forAll(pNu, j)
        {
            pNu[j] = ... // patch cell face
        }
    }

    return tNu;
}

[Svensen]

Dear Zeppo,

thank you very much for help !

I've inserted your code in Spline.C and fixed some syntax issues in it. It compiles fine, however during execution the error with reference occurred. I've tried to comment every line of code from end to beginning and it seems that something wrong with this line:

tmp<volScalarField> tNu
    (
        new volScalarField
        (
            IOobject
            (
                "nu",
                strain_Rate.time().timeName(),
                strain_Rate.mesh(),
                strain_Rate.readOpt(),
                strain_Rate.writeOpt()
            ),
            strain_Rate.mesh(),
            sqr(dimLength)*strain_Rate.dimensions(),
            strain_Rate.boundaryField().types()
        )
    );

error log states that:

Selecting incompressible transport model Spline
[3] #0  Foam::error::printStack(Foam::Ostream&)[4] 
[4] 
[4] --> FOAM FATAL ERROR: 
[4] hanging pointer at index 0 (size 6), cannot dereference
[4] 
[4]     From function const T& Foam::UPtrList<T>::operator[](Foam::label) const [with T = Foam::fvPatchField<double>; Foam::label = int]
[4]     in file /home/sergey/OpenFOAM-dev/OpenFOAM-dev/src/OpenFOAM/lnInclude/UPtrListI.H at line 107.
[4] 
FOAM parallel run aborting

I've attached modified sources as well as full log below

[Zeppo]

Ok. Then replace your calcNu() with the code following below:

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::Spline::calcNu() const
{
    const volScalarField& strainRate = strainRate()();
    tmp<volScalarField> tNu
    (
        new volScalarField
        (
            IOobject
            (
                "nu",
                strainRate.time().timeName(),
                strainRate.mesh(),
                strainRate.readOpt(),
                strainRate.writeOpt()
            ),
            strainRate.mesh(),
            dimensionedScalar
            (   
                "oneNu", 
                sqr(dimLength)*strainRate.dimensions(), 
                1.0
            )
        )
    );
    return tNu;
}

[Svensen]

Dear Zeppo, thank you very much for help !

You suggestion is right and I was able to implement just a simple Newtonian model:

forAll(nu, i){
        nu[i] = visc[0];
    }

This code works fine when nu is not depend on strainRate. However, if I need a viscosity as a function of strain rate:

forAll(nu, i){
        nu[i] = function(strain_Rate[i]);
    }

it crashes.

It seems to me that i-index for nu is not the same as i-index for strain_Rate.

Maybe you know how the correct index for strain_Rate can be obtained ?
The sources as well as log file are attached below

[Zeppo]

Both strain_Rate and nu have the same number of elements, so the code like this

forAll(nu, i)
{
    Info << nu[i] << " " << strain_Rate[i] <<endl;   
}

should be unconditionally valid. Commenting out the code you should isolate the problem to the single line which is the cause of the error

[Svensen]

Finally I've figured out a problem. I' found that this code works properly on a normal execution like "pimpleFoam", but code crashes on parallel run "pimpleFoam -parallel".

Maybe you know if there are some additional changes are required to source code to support parallel execution ?

Thanks !

[Zeppo]

Which line triggers the error? This one?

double n = strain_Rate[i] / dn;

[Svensen]

Yes, you are right. Something wrong here.

[Zeppo]

run this

forAll(nu, i)
{
    Info << i << ": " << nu[i] << " " << strain_Rate[i] <<endl;   
}

[Svensen]

In the first step for every nu it prints i: 1 1 like

10380: 1 1

and so on.

At the beginning of second step it crashes before any output starts.
The full log is attached

[Zeppo]

I am not an expert but what I see is that the error message never mentions your library file (.so-file). You built your custom viscosity class into a user library, right? It is libOpenFOAM.so that is mensioned in the error message as the sourse where the error comes from.
You said that the very same code runs flawlessly in serial, correct?

[Svensen]

Yes, serial execution works fine.
I will post this issue on issue tracker.

Thanks a lot for help !

[wyldckat]

Quick answer:

Originally Posted by Svensen

Yes, serial execution works fine.
I will post this issue on issue tracker.

Not a bug in OpenFOAM, it's a limitation of the implementation you've done

I need to understand what is the "visc" vector meant to be. Is it a lookup table?

[Svensen]

Yes, it is just a viscosity data. visc[0] contains the kinematic viscosity for shearRate=0; visc[1] contains the kinematic viscosity for shearRate=0.125; visc[2] contains viscosity for shearRate=0.250 and so on.

for shearRate between points, for example for shearRate=0.2 the linear interpolation will be used to compute nu for this point.

The sources are attached above.

[Svensen]

Can this problem be treated ?

[wyldckat]

Greetings Svensen,

There are a few boundary conditions that demonstrate how lookup tables can be used. Using lookup tables as OpenFOAM already does, would ensure that the bug is not coming from a possibly flawed implementation of the lookup mechanism you currently have implemented.

The boundary "src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValue" has such an example:

1. The lookup table object is defined like this in the header file ".H":
   

   autoPtr<Function1<Type>> uniformValue_;

   • In your case:
     

     autoPtr<Function1<scalar>> visc;

2. Constructed like this:
   

   uniformValue_(Function1<Type>::New("uniformValue", dict))

   • In your case:
     

     visc(Function1<scalar>::New("viscTable", dict))

     Then the entry "viscTable" needs to be in "transportProperties", within the respective "*Coefs" section.
3. Then used like this:
   

       const scalar t = this->db().time().timeOutputValue();
       fvPatchField<Type>::operator==(uniformValue_->value(t));

   • Although in your case, it would be something like:
     

     visc->value(strainRate[celli])

4. Don't forget to include the header file:
   

   #include "Function1.H"

5. In the "transportProperties" file, you will need something like what's described here: http://openfoam.org/release/2-1-0/bo...ime-dependent/ - possibly this:
   

   viscTable     table
   (
       (0.0     0.014)
       (0.125   0.0123715)
       (0.25    0.011782)
   ...
       (????  0.00409976)
   );

If you still have trouble implementing it, please do the following steps, so it's easier to help you:

1. Provide the source code by following these steps:
   1. Run inside the source code folder of your custom library:
      

      wclean all

   2. Then go to the parent folder and package the source code folder, e.g. if the folder is named "splineTransport":
      

      tar -czf splineLib.tar.gz splineTransport

   3. Then attach the file "splineLib.tar.gz" to your next post.
2. Then we either need a test case, or instructions on how to change a tutorial in OpenFOAM to use your transport model.

Best regards,
Bruno