|
[Sponsors] |
November 2, 2016, 07:14 |
Custom viscosity model
|
#1 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
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: Code:
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; } Code:
#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 ??:? |
|
November 2, 2016, 19:01 |
|
#2 | ||
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,978
Blog Entries: 45
Rep Power: 128 |
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: 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; } Quote:
If we take a look at the printed stack trace, the first tell-tale sign is this: Code:
#1 Foam::sigFpe::sigHandler(int) at ??:? Quote:
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:
Code:
Info << "size visc: " << sizeof(visc) << endl; Info << "dn: " << dn << endl; Code:
forAll(strain, itemi) { Info << strain[itemi] << endl; } Bruno
__________________
|
|||
November 3, 2016, 07:12 |
|
#3 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
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. |
|
November 3, 2016, 09:04 |
|
#4 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
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: Code:
return dimensionedScalar("one", pow(dimLength, 2), 1.0)*strain; Code:
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). |
|
November 4, 2016, 17:35 |
|
#5 |
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 |
Code:
Foam::tmp<Foam::volScalarField> Foam::viscosityModels::Spline::calcNu() const { ... volScalarField strain(strainRate()); ... return dimensionedScalar("one", dimLength*dimLength, 1.0)*strain; } Code:
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; } |
|
November 4, 2016, 18:53 |
|
#6 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
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: Code:
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() ) ); Code:
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 |
|
November 5, 2016, 05:11 |
|
#7 |
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 |
Ok. Then replace your calcNu() with the code following below:
Code:
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; } |
|
November 5, 2016, 05:48 |
|
#8 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
Dear Zeppo, thank you very much for help !
You suggestion is right and I was able to implement just a simple Newtonian model: Code:
forAll(nu, i){ nu[i] = visc[0]; } Code:
forAll(nu, i){ nu[i] = function(strain_Rate[i]); } 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 |
|
November 5, 2016, 06:17 |
|
#9 |
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 |
Both strain_Rate and nu have the same number of elements, so the code like this
Code:
forAll(nu, i) { Info << nu[i] << " " << strain_Rate[i] <<endl; } |
|
November 5, 2016, 06:55 |
|
#10 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
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 ! |
|
November 5, 2016, 09:40 |
|
#11 |
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 |
Which line triggers the error? This one?
Code:
double n = strain_Rate[i] / dn; |
|
November 5, 2016, 10:43 |
|
#12 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
Yes, you are right. Something wrong here.
|
|
November 5, 2016, 12:11 |
|
#13 |
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 |
run this
Code:
forAll(nu, i) { Info << i << ": " << nu[i] << " " << strain_Rate[i] <<endl; } |
|
November 5, 2016, 12:55 |
|
#14 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
In the first step for every nu it prints i: 1 1 like
Code:
10380: 1 1 At the beginning of second step it crashes before any output starts. The full log is attached |
|
November 5, 2016, 14:21 |
|
#15 |
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 |
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? |
|
November 5, 2016, 14:26 |
|
#16 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
Yes, serial execution works fine.
I will post this issue on issue tracker. Thanks a lot for help ! |
|
November 5, 2016, 15:17 |
|
#18 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
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. |
|
November 6, 2016, 03:42 |
|
#19 |
Senior Member
Join Date: Jan 2015
Posts: 150
Rep Power: 11 |
Can this problem be treated ?
|
|
November 6, 2016, 04:21 |
|
#20 |
Retired Super Moderator
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,978
Blog Entries: 45
Rep Power: 128 |
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:
If you still have trouble implementing it, please do the following steps, so it's easier to help you:
Bruno
__________________
|
|
Tags |
openfoam-dev, viscosity model |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
Problem with divergence | TDK | FLUENT | 13 | December 14, 2018 06:00 |
New Viscosity Model - Best Practice for Temporary Values | MaSch | OpenFOAM Programming & Development | 6 | March 16, 2018 03:51 |
Time constant in Herschel-Bulkley viscosity model | Mikel6 | Main CFD Forum | 0 | October 17, 2016 04:52 |
Viscosity ratio in gamma-theta transition model based on k-w sst turb model | Qiaol618 | Main CFD Forum | 8 | June 9, 2012 06:43 |
How to modify the viscosity model | mpml | OpenFOAM Running, Solving & CFD | 4 | October 13, 2010 07:44 |