Rate this Entry

Excellent tutorial for debugging

Posted March 10, 2017 at 11:50 by kindle

banana and old-school printing out.

Cool !

Quote:

Originally Posted by wyldckat

OK, so as I mentioned in one of the posts above, having an example case as a basis and using the "banana" trick, here's what I've done (used OpenFOAM 2.3.x for this example):

Made a copy of the tutorial "lagrangian/DPMFoam/Goldschmidt" and ran the utility blockMesh.
Then I had to look at which file needed to modified. All pointed to the file "constant/kinematicCloudProperties".
Now I had to figure out where the Lift comes in. After going around looking at the other tutorials and a few threads, the conclusion was that this block of dictionary is where we can stack up the various forces at work:
Code:
```
    particleForces
    {
        ErgunWenYuDrag
        {
            alphac alpha.air;
        }
        gravity;
    }
```

And here's where the "banana" comes in:

Code:

    particleForces
    {
        ErgunWenYuDrag
        {
            alphac alpha.air;
        }
        gravity;
        banana;
    }

As soon as I run DPMFoam, it complains about "banana" not being valid:

Code:

--> FOAM FATAL ERROR: 
Unknown particle force type banana, constructor not in hash table

    Valid particle force types are:

13
(
ErgunWenYuDrag
PlessisMasliyahDrag
SRF
SaffmanMeiLiftForce
TomiyamaLift
WenYuDrag
gravity
nonInertialFrame
nonSphereDrag
paramagnetic
pressureGradient
sphereDrag
virtualMass
)

There's our pretty lift model "SaffmanMeiLiftForce".

OK, rename "banana" to "SaffmanMeiLiftForce" and added brackets with a "banana" inside

Code:

    particleForces
    {
        ErgunWenYuDrag
        {
            alphac alpha.air;
        }
        gravity;
        SaffmanMeiLiftForce
        {
            banana;
        }
    }

mmm... the solver crashed. So much for the "banana".

OK, if we look at the source code for this model:

Code:

template<class CloudType>
Foam::SaffmanMeiLiftForce<CloudType>::SaffmanMeiLiftForce
(
    CloudType& owner,
    const fvMesh& mesh,
    const dictionary& dict,
    const word& forceType
)
:
    LiftForce<CloudType>(owner, mesh, dict, forceType)
{}

there is no coefficient necessary.

Therefore, this should be enough:

Code:

    particleForces
    {
        ErgunWenYuDrag
        {
            alphac alpha.air;
        }
        gravity;
        SaffmanMeiLiftForce;
    }

Er, nope... it complains it should be a dictionary itself. Then back to this:

Code:

    particleForces
    {
        ErgunWenYuDrag
        {
            alphac alpha.air;
        }
        gravity;
        SaffmanMeiLiftForce
        {
        }
    }

So we go back to the crash. The error given is this:
Code:
```
--> FOAM FATAL ERROR: 

    request for volVectorField U from objectRegistry region0 failed
    available objects of type volVectorField are
1(U.air)
```
OK, it's looking for the "U" field and it only found "U.air". So, either the model is just wrong, or something else is missing.

Looking at the class from which this model derives, the constructor has a vital clue to what's missing:

Code:

template<class CloudType>
Foam::LiftForce<CloudType>::LiftForce
(
    CloudType& owner,
    const fvMesh& mesh,
    const dictionary& dict,
    const word& forceType
)
:
    ParticleForce<CloudType>(owner, mesh, dict, forceType, true),
    UName_(this->coeffs().template lookupOrDefault<word>("U", "U")),
    curlUcInterpPtr_(NULL)
{}

It's looking for a new name for "U", something like:

Code:

U U.air;

So let's try this:

Code:

    particleForces
    {
        ErgunWenYuDrag
        {
            alphac alpha.air;
        }
        gravity;
        SaffmanMeiLiftForce
        {
            U U.air;
        }
    }

Side note: curiously enough, Maysam got a similar answer a few minutes ago in post #4 at http://www.cfd-online.com/Forums/ope...tml#post512965

Running DPMFoam once again gives this error message:

Code:

--> FOAM FATAL IO ERROR: 
keyword curlUcDt is undefined in dictionary "/home/ofuser/OpenFOAM/ofuser-2.3.x/run/tutorials/lagrangian/DPMFoam/Goldschmidt/constant/kinematicCloudProperties.solution.interpolationSchemes"

file: /home/ofuser/OpenFOAM/ofuser-2.3.x/run/tutorials/lagrangian/DPMFoam/Goldschmidt/constant/kinematicCloudProperties.solution.interpolationSchemes from line 27 to line 29.

OK, we need to look at this this block of dictionary:

Code:

solution
{
    active          true;
    coupled         true;
    transient       yes;
    cellValueSourceCorrection off;

    interpolationSchemes
    {
        rho.air         cell;
        U.air           cellPoint;
        mu.air          cell;
    }

    integrationSchemes
    {
        U               Euler;
    }

    sourceTerms
    {
        schemes
        {
            U semiImplicit 1;
        }
    }
}

Let's try:

Code:

    interpolationSchemes
    {
        rho.air         cell;
        U.air           cellPoint;
        mu.air          cell;
        curlUcDt        banana;
    }

The error message stated by DPMFoam:

Code:

--> FOAM FATAL ERROR: 
Unknown interpolation type banana for field curlUcDt

Valid interpolation types : 

6
(
cell
cellPatchConstrained
cellPoint
cellPointFace
cellPointWallModified
pointMVC
)

I have no idea. Let's try the most conservative/complex one "cellPointFace"... nope, it crashes. Let's try all others... nope, all of them crash.

The error message of the crash is something like this:

Code:

Solving 3-D cloud kinematicCloud
#0  Foam::error::printStack(Foam::Ostream&) in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigFpe::sigHandler(int) in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::SaffmanMeiLiftForce<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > >::Cl(Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > const&, Foam::Vector<double> const&, double, double) const in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/liblagrangianIntermediate.so"
#4  Foam::LiftForce<Foam::KinematicCloud<Foam::Cloud<Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > > > >::calcCoupled(Foam::CollidingParcel<Foam::KinematicParcel<Foam::particle> > const&, double, double, double, double) const in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/lib/liblagrangianIntermediate.so"
#5  
 in "/home/ofuser/OpenFOAM/OpenFOAM-2.3.x/platforms/linux64GccDPOpt/bin/DPMFoam"
#6

and it continues with more information, but doesn't help much more...
So there was a "sigFpe", hence division by zero or infinite or something like that, in the method "::Cl" (it's in line #3).

OK, looking at the source code for this method:
there are a few possible locations where it might have broken.

OK, let's go old-school debugging:

Code:

template<class CloudType>
Foam::scalar Foam::SaffmanMeiLiftForce<CloudType>::SaffmanMeiLiftForce::Cl
(
    const typename CloudType::parcelType& p,
    const vector& curlUc,
    const scalar Re,
    const scalar muc
) const
{
    scalar Rew = p.rhoc()*mag(curlUc)*sqr(p.d())/(muc + ROOTVSMALL);
    Info << "Rew: " << Rew << endl;
    
    scalar beta = 0.5*(Rew/(Re + ROOTVSMALL));
    Info << "Re: " << Re << endl;
    Info << "beta: " << beta << endl;

    scalar alpha = 0.3314*sqrt(beta);
    Info << "alpha: " << alpha << endl;

    scalar f = (1.0 - alpha)*exp(-0.1*Re) + alpha;
    Info << "f: " << f << endl;

    scalar Cld = 0.0;
    if (Re < 40)
    {
        Cld = 6.46*f;
    }
    else
    {
        Cld = 6.46*0.0524*sqrt(beta*Re);
    }

    Info << "Cld: " << Cld << endl;

    return 3.0/(mathematical::twoPi*sqrt(Rew))*Cld;
}

we print out everything so that we can have a look at each value. It's good to know what's happening under-the-hood, instead of having to guess: http://openfoamwiki.net/index.php/HowTo_debugging

And we have to build the modified library "src/lagrangian/intermediate":
Code:
```
wmake libso $FOAM_SRC/lagrangian/intermediate
```
Let's have another go with DPMFoam... the output we now get:
Code:
```
Solving 3-D cloud kinematicCloud
Rew: 0
Re: 71.5702
beta: 0
alpha: 0
f: 0.00077937
Cld: 0
```
that's waaaay too many zeros... this leads to this expression to go haywire:
Code:
```
return 3.0/(mathematical::twoPi*sqrt(Rew))*Cld;
```
Because it becomes:
Code:
```
3.0/(2*Pi*sqrt(0.0))*0.0
```
where it divides by zero.
OK, we can try and insert a division protection like it's done for the others:
Code:
```
return 3.0/(mathematical::twoPi*sqrt(Rew+SMALL))*Cld;
```
Build the library again and gave DPMFoam another run... didn't crash, but gave a lot of junk output. So I Ctrl+C to cancel the run and went back to the code and removed the "Info" lines I had added. And built the library again.
Another run with DPMFoam and... and it's finally working! As to whether it's giving good or bad values, I have no idea.

In the meantime I'll report this bug, although I'm not sure if the fix I presented is correct or not

.

-------------------
edit: Bug reported here: http://www.openfoam.org/mantisbt/view.php?id=1408

Additional note: Letting the case run first without any particles should initialize the "U.air" field, should no longer require you to modify the source code, since the flow field would all not be zero.
-------------------

Best regards,
Bruno

Posted in programming

Views 1519 Comments 0

« Prev Main Next »

Total Comments 0

Archive
< April 2024
Su	Mo	Tu	We	Th	Fr	Sa
24	25	26	27	28	29	30
31	1	2	3	4	5	6
7	8	9	10	11	12	13
14	15	16	17	18	19	20
21	22	23	24	25	26	27
28	29	30	1	2	3	4

CFD Online Discussion Forum Message

Cancel Changes