CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

help diagnosing alphaContactAngle error message

Register Blogs Community New Posts Updated Threads Search

Like Tree1Likes
  • 1 Post By frantov

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 28, 2020, 13:43
Default help diagnosing alphaContactAngle error message
  #1
Senior Member
 
Josh McCraney
Join Date: Jun 2018
Posts: 220
Rep Power: 9
joshmccraney is on a distinguished road
Hi FOAMers!

I am on OpenFOAM 6. I found an alphaContactAngle model that I'm trying to implement (Kistler contact angle model), the file called dynamicKistlerAlphaContactAngleFvPatchScalarField. H (and .C). It compiles, but issues the warning:

Code:
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:160:33: warning: unused variable ‘sigmap’ [-Wunused-variable]
     const fvPatchField<scalar>& sigmap =
                                 ^~~~~~
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:163:33: warning: unused variable ‘mup’ [-Wunused-variable]
     const fvPatchField<scalar>& mup =
                                 ^~~
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:202:12: warning: unused variable ‘InvHoffFuncThetaAroot’ [-Wunused-variable]
     scalar InvHoffFuncThetaAroot = RRInvHoffFuncThetaA.root(0,65);
            ^~~~~~~~~~~~~~~~~~~~~
alphaContactAngle/dynamicDavisHockingAlphaContactAngle/dynamicDavisHockingAlphaContactAngleFvPatchScalarField.C:205:12: warning: unused variable ‘InvHoffFuncThetaRroot’ [-Wunused-variable]
     scalar InvHoffFuncThetaRroot = RRInvHoffFuncThetaR.root(0,65);
In order to use this model I have to use a different solver called kistlerInterFoam (basically interFoam with some small tweaks). This compiles fine with no warnings. When I run the case, I get the error

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  ? in "/lib/x86_64-linux-gnu/libm.so.6"
#4  pow in "/lib/x86_64-linux-gnu/libm.so.6"
#5  Foam::dynamicKistlerAlphaContactAngleFvPatchScalarField::HoffmanFunction(double const&) const at ??:?
#6  Foam::dynamicKistlerAlphaContactAngleFvPatchScalarField::theta(Foam::fvPatchField<Foam::Vector<double> > const&, Foam::fvsPatchField<Foam::Vector<double> > const&) const at ??:?
#7  Foam::interfaceProperties::correctContactAngle(Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>::Boundary&, Foam::GeometricField<Foam::Vector<double>, Foam::fvsPatchField, Foam::surfaceMesh>::Boundary const&) const at ??:?
#8  Foam::interfaceProperties::calculateK() at ??:?
#9  Foam::interfaceProperties::interfaceProperties(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::IOdictionary const&) at ??:?
#10  Foam::immiscibleIncompressibleTwoPhaseMixture::immiscibleIncompressibleTwoPhaseMixture(Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&, Foam::GeometricField<double, Foam::fvsPatchField, Foam::surfaceMesh> const&) at ??:?
#11  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/kistlerInterFoam"
#12  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#13  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/kistlerInterFoam"
./run.sh: line 9:  8374 Floating point exception(core dumped) kistlerInterFoam
Here's what I know

Line #1 contains an Fpe error, which google implies is division by zero.
Line #4 includes pow, which appears only once in the dynamicKistlerAlphaContactAngleFvPatchScalarField. C and .H files
Line #5 references a HoffmanFunction, which is also contained in dynamicKistlerAlphaContactAngleFvPatchScalarField. H and .C files.

Below is the associated .C file, where I highlight in green the Hoffman function and pow components and immediately before that a division step by sigmap, which again threw a warning in compiling.

Code:
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

\*---------------------------------------------------------------------------*/

#include "dynamicKistlerAlphaContactAngleFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchFieldMapper.H"
#include "fvPatchFields.H"
#include "volMesh.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //

const scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::convertToDeg =
    180.0/constant::mathematical::pi;

const scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::convertToRad =
    constant::mathematical::pi/180.0;

const scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::theta0 = 90.0;

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    alphaContactAngleFvPatchScalarField(p, iF),
    thetaA_(0.0),
    thetaR_(0.0),
    muName_("undefined"),
    sigmaName_("undefined")
{}

dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const dynamicKistlerAlphaContactAngleFvPatchScalarField& acpsf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    alphaContactAngleFvPatchScalarField(acpsf, p, iF, mapper),
    thetaA_(acpsf.thetaA_),
    thetaR_(acpsf.thetaR_),
    muName_(acpsf.muName_),
    sigmaName_(acpsf.sigmaName_)
{}


dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    alphaContactAngleFvPatchScalarField(p, iF),
    thetaA_(readScalar(dict.lookup("thetaA"))),
    thetaR_(readScalar(dict.lookup("thetaR"))),
    muName_(dict.lookup("muEffKistler")),
    sigmaName_(dict.lookup("sigmaKistler"))
{
    evaluate();
}


dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const dynamicKistlerAlphaContactAngleFvPatchScalarField& acpsf
)
:
    alphaContactAngleFvPatchScalarField(acpsf),
    thetaA_(acpsf.thetaA_),
    thetaR_(acpsf.thetaR_),
    muName_(acpsf.muName_),
    sigmaName_(acpsf.sigmaName_)
{}


dynamicKistlerAlphaContactAngleFvPatchScalarField::
dynamicKistlerAlphaContactAngleFvPatchScalarField
(
    const dynamicKistlerAlphaContactAngleFvPatchScalarField& acpsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    alphaContactAngleFvPatchScalarField(acpsf, iF),
    thetaA_(acpsf.thetaA_),
    thetaR_(acpsf.thetaR_),
    muName_(acpsf.muName_),
    sigmaName_(acpsf.sigmaName_)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

tmp<scalarField> dynamicKistlerAlphaContactAngleFvPatchScalarField::theta
(
    const fvPatchVectorField& Up,
    const fvsPatchVectorField& nHat
) const
{
    //eb - Lookup and return the patchField of dynamic viscosity of mixture
    //     and surface tension
/*    
    if((muName_ != "muEffKistler") || (sigmaName_ != "sigmaKistler"))
    {
        FatalErrorIn
        (
            "dynamicKistlerAlphaContactAngleFvPatchScalarField"
        )   << " muEffKistler or sigma set inconsitently, muEffKistler = "
            << muName_ << ", sigmaKistler = " << sigmaName_ << '.' << nl
            << "    Set both muEffKistler and sigmaKistler according to the "
            << "definition of dynamicKistlerAlphaContactAngle"
            << "\n    on patch " << this->patch().name()
            << " of field " << this->dimensionedInternalField().name()
            << " in file " << this->dimensionedInternalField().objectPath()
            << exit(FatalError);
    }
*/

    const fvPatchField<scalar>& sigmap =
        patch().lookupPatchField<volScalarField, scalar>(sigmaName_);

    const fvPatchField<scalar>& mup =
        patch().lookupPatchField<volScalarField, scalar>(muName_);

    vectorField nf = patch().nf();

    // Calculate the component of the velocity parallel to the wall
    vectorField Uwall = Up.patchInternalField() - Up;
    Uwall -= (nf & Uwall)*nf;

    // Find the direction of the interface parallel to the wall
    vectorField nWall = nHat - (nf & nHat)*nf;

    // Normalise nWall
    nWall /= (mag(nWall) + SMALL);

    // Calculate Uwall resolved normal to the interface parallel to
    // the wall
    scalarField uwall = nWall & Uwall;

    //eb - Calculate local Capillary number
    scalarField Ca = mup*mag(uwall)/sigmap;
 
    //eb - Instantiate function object InverseHoffmanFunction for thetaA and thetaR
    dynamicKistlerAlphaContactAngleFvPatchScalarField::InverseHoffmanFunction
    InvHoffFuncThetaA
    (
        convertToRad*thetaA_
    );

    dynamicKistlerAlphaContactAngleFvPatchScalarField::InverseHoffmanFunction
    InvHoffFuncThetaR
    (
        convertToRad*thetaR_
    );

   //eb - Calculate InverseHoffmanFunction for thetaA and thetaR using
    // RiddersRoot
    RiddersRoot RRInvHoffFuncThetaA(InvHoffFuncThetaA, 1.e-10);
    scalar InvHoffFuncThetaAroot = RRInvHoffFuncThetaA.root(0,65);

    RiddersRoot RRInvHoffFuncThetaR(InvHoffFuncThetaR, 1.e-10);
    scalar InvHoffFuncThetaRroot = RRInvHoffFuncThetaR.root(0,65);

    //eb - Calculate and return the value of contact angle on patch faces,
    //     a general approach: the product of Uwall and nWall is negative
    //     for advancing and positive for receding motion.
    //     thetaDp is initialized to 90 degrees corresponding to no wall
    //     adhesion
    scalarField thetaDp(patch().size(), convertToRad*theta0);
    forAll(uwall, pfacei)
    {
        if(uwall[pfacei] < 0.0)
        {
            //thetaDp[pfacei] = convertToRad*thetaA_;
            thetaDp[pfacei] = HoffmanFunction(   Ca[pfacei]
                                               + InvHoffFuncThetaAroot);
        }
        else if (uwall[pfacei] > 0.0)
        {
            //thetaDp[pfacei] = convertToRad*thetaR_;
            thetaDp[pfacei] = HoffmanFunction(   Ca[pfacei]
                                               - InvHoffFuncThetaRroot);
        }
    }
    
    return convertToDeg*thetaDp;
}


scalar dynamicKistlerAlphaContactAngleFvPatchScalarField::HoffmanFunction
(
    const scalar& x
) const
{
    return acos(1 - 2*tanh(5.16*pow(x/(1+1.31*pow(x,0.99)),0.706)));
}


void dynamicKistlerAlphaContactAngleFvPatchScalarField::write(Ostream& os) const
{
    fvPatchScalarField::write(os);
    os.writeKeyword("thetaA") << thetaA_ << token::END_STATEMENT << nl;
    os.writeKeyword("thetaR") << thetaR_ << token::END_STATEMENT << nl;
    os.writeKeyword("muEffKistler") << muName_ << token::END_STATEMENT << nl;
    os.writeKeyword("sigmaKistler") << sigmaName_ << token::END_STATEMENT << nl;
    writeEntry("value", os);
}


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

makePatchTypeField
(
    fvPatchScalarField,
    dynamicKistlerAlphaContactAngleFvPatchScalarField
);

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// ************************************************************************* //
In case this is relevant, the dynamicKistlerAlphaContactAngleFvPatchScalarField. H file is below, I've highlighted the HoffmanFunction part of the code in green.

Code:
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Class
    dynamicKistlerAlphaContactAngleFvPatchScalarField

Description
    Evaluation of the dynamic contact angle using the Hoffman function,
    based on Kistler (1993)

SourceFiles
    dynamicKistlerAlphaContactAngleFvPatchScalarField.C

Author
    Edin Berberovic 2008 - Updated by Alexander Rattner, 2012

\*---------------------------------------------------------------------------*/

#ifndef dynamicKistlerAlphaContactAngleFvPatchScalarField_H
#define dynamicKistlerAlphaContactAngleFvPatchScalarField_H

#include "alphaContactAngleFvPatchScalarField.H"
#include "mathematicalConstants.H"
#include "surfaceFields.H"
#include "RiddersRoot.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                  Class dynamicAlphaContactAngleFvPatch Declaration
\*---------------------------------------------------------------------------*/

class dynamicKistlerAlphaContactAngleFvPatchScalarField
:
    public alphaContactAngleFvPatchScalarField
{
    // Private data

        //- Limiting advancing contact angle
        scalar thetaA_;

        //- Limiting receeding contact angle
        scalar thetaR_;

        //eb Name of the dynamic viscosity field
        word muName_;

        //eb Name of the surface tension
        word sigmaName_;

public:

    //eb - Function object class, inverse of the Hoffman function
    class InverseHoffmanFunction : public FuncBase
    {

        // Private data

            //- Equilibrium contact angle (advancing or receding), as parameter
            scalar thetaE_;

    public:

        // Constructors

            //- Construct from data
            InverseHoffmanFunction(const scalar& thetaE)
            :
                thetaE_(thetaE)
            {}

        // operator()
            scalar operator()(scalar fHI) const
            {
                return
                  pow(fHI,0.706)
                -   pow((1+1.31*pow(fHI,0.99)),0.706)
                   *(1.0/5.16)*atanh((1-cos(thetaE_))/2);
            }
    };

public:

    //- Runtime type information
    TypeName("dynamicKistlerAlphaContactAngle");

    // Static data members

        //eb - Conversion factor for radians into degrees
        static const scalar convertToDeg;

        //eb - Conversion factor for degrees into radians
        static const scalar convertToRad;

        //eb - Contact angle corresponding to zero contact line velocity
        static const scalar theta0;

    // Constructors

        //- Construct from patch and internal field
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct from patch, internal field and dictionary
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const dictionary&
        );

        //- Construct by mapping given
        //  dynamicKistlerAlphaContactAngleFvPatchScalarField
        //  onto a new patch
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const dynamicKistlerAlphaContactAngleFvPatchScalarField&,
            const fvPatch&,
            const DimensionedField<scalar, volMesh>&,
            const fvPatchFieldMapper&
        );

        //- Construct as copy
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const dynamicKistlerAlphaContactAngleFvPatchScalarField&
        );

        //- Construct and return a clone
        virtual tmp<fvPatchScalarField> clone() const
        {
            return tmp<fvPatchScalarField>
            (
                new dynamicKistlerAlphaContactAngleFvPatchScalarField(*this)
            );
        }

        //- Construct as copy setting internal field reference
        dynamicKistlerAlphaContactAngleFvPatchScalarField
        (
            const dynamicKistlerAlphaContactAngleFvPatchScalarField&,
            const DimensionedField<scalar, volMesh>&
        );

        //- Construct and return a clone setting internal field reference
        virtual tmp<fvPatchScalarField> clone
        (
            const DimensionedField<scalar, volMesh>& iF
        ) const
        {
            return tmp<fvPatchScalarField>
            (
                new dynamicKistlerAlphaContactAngleFvPatchScalarField(*this, iF)
            );
        }


    // Member functions

        //- Evaluate and return dynamic contact-angle
        virtual tmp<scalarField> theta
        (
            const fvPatchVectorField& Up,
            const fvsPatchVectorField& nHat
        ) const;

        //eb - Evaluate the value of the Hoffman function
        scalar HoffmanFunction
        (
            const scalar& x
        ) const;

        //- Write
        virtual void write(Ostream&) const;
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
Any help is greatly appreciated!

Last edited by joshmccraney; May 28, 2020 at 17:56.
joshmccraney is offline   Reply With Quote

Old   August 18, 2020, 03:39
Default
  #2
Member
 
Francisco T
Join Date: Nov 2011
Location: Melbourne, Australia
Posts: 64
Blog Entries: 1
Rep Power: 15
frantov is on a distinguished road
Hello, have you solved your problem?


Im also trying to implemen Kistler contact angle model in OF7 but I get this error.error:



‘Foam::RiddersRoot’ is not a template
194 | RiddersRoot<InverseHoffmanFunction> RRInvHoffFuncThetaA(InvHoffFuncThetaA, 1.e-10);



Im thinking the RiddlersRoot.H and C. files I have need to be updated. Did you updated this files?
Cheers
Fran
frantov is offline   Reply With Quote

Old   August 18, 2020, 04:50
Default
  #3
Member
 
Francisco T
Join Date: Nov 2011
Location: Melbourne, Australia
Posts: 64
Blog Entries: 1
Rep Power: 15
frantov is on a distinguished road
I found some updated .H .C files and did some conversion to OF7.
The library now compiled well..no errors
Now Im trying to link that library with the solver....
Cheers
joshmccraney likes this.
frantov is offline   Reply With Quote

Old   September 16, 2020, 10:33
Default
  #4
Senior Member
 
Josh McCraney
Join Date: Jun 2018
Posts: 220
Rep Power: 9
joshmccraney is on a distinguished road
Quote:
Originally Posted by frantov View Post
I found some updated .H .C files and did some conversion to OF7.
The library now compiled well..no errors
Now Im trying to link that library with the solver....
Cheers
Sorry, I'm just now seeing these two posts. I ended up moving to OpenFOAM 4.X. If you got it working in OF7, could you post your scripts somewhere? Would be very helpful to see!
joshmccraney is offline   Reply With Quote

Reply


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
alphaContactAngle for reactingMultiphaseEulerFoam amin.z OpenFOAM Running, Solving & CFD 1 September 28, 2021 06:47
Duplicate entry alphaContactAngle in runtime selection table fvPatchField Roesch OpenFOAM Programming & Development 2 January 12, 2021 04:35
alphaContactAngle parameters: in which file to define? makaveli_lcf OpenFOAM Running, Solving & CFD 0 April 8, 2013 09:25
alphaContactAngle not recognized robosys541 OpenFOAM Running, Solving & CFD 0 August 21, 2012 17:41
alphaContactAngle patch type does not work in OF-1.6 sundaero OpenFOAM Bugs 10 April 13, 2012 19:08


All times are GMT -4. The time now is 10:38.