Custom viscosity model
 

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:


It compiles with no error, however after start of pimpleFoam it exits with the following log:
What can be wrong with my code ?

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:
Furthermore you commented that:
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:
This means that there was a https://en.wikipedia.org/wiki/SIGFPE - which quoting from there:
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:
  
  Code:

  ---

          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.:
Beyond this, the way that OpenFOAM goes through lists is by using forAll, e.g.:
Best regards,
Bruno

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.

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:
Error log is the following:
It seems to me that something wrong happens when I return the "strain" object from calcNu() method. Maybe some pointer violation or type mismatch...

->

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:
error log states that:
I've attached modified sources as well as full log below

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

Dear Zeppo, thank you very much for help !

You suggestion is right and I was able to implement just a simple Newtonian model:
This code works fine when nu is not depend on strainRate. However, if I need a viscosity as a function of strain rate:
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

Both strain_Rate and nu have the same number of elements, so the code like this
should be unconditionally valid. Commenting out the code you should isolate the problem to the single line which is the cause of the error

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 !

Which line triggers the error? This one?

Yes, you are right. Something wrong here.

run this

In the first step for every nu it prints i: 1 1 like and so on.

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

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?

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

Thanks a lot for help !

Quick answer:
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?

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.

Can this problem be treated ?

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":
   
   Code:

   ---

   autoPtr<Function1<Type>> uniformValue_;

   ---

   • In your case:
     
     Code:

     ---

     autoPtr<Function1<scalar>> visc;

     ---

2. Constructed like this:
   
   Code:

   ---

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

   ---

   • In your case:
     
     Code:

     ---

     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:
   
   Code:

   ---

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

   ---

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

     ---

     visc->value(strainRate[celli])

     ---

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

   ---

   #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:
   
   Code:

   ---

   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:
      
      Code:

      ---

      wclean all

      ---

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

      ---

      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