[juhuettn]

Hi all,

I am trying to implement a custom wall function to use them with the pisoFoam solver. I compiled the new modified wall functions with wmake libso without problems, but I cannot select them when I am running a case. I get an error telling me I am using an unknown patchField type.
I have included the link to my user library in the controlDict file of the case and I also get no errors during compiling. Does anybody have an idea what mistake I am doing?

The Make/files file for the RAS model:

/* Wall functions */
wallFunctions = derivedFvPatchFields/wallFunctions

nutWallFunctions = $(wallFunctions)/nutWallFunctions
$(nutWallFunctions)/nutDurbinWallFunction/nutDurbinWallFunctionFvPatchScalarField.C

epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
$(epsilonWallFunctions)/epsilonDurbinWallFunction/epsilonDurbinWallFunctionFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/mylibincompressibleRASModels

The Make/options file:

EXE_INC = \
    -I$(LIB_SRC)/turbulenceModels \
    -I$(LIB_SRC)/transportModels \
    -I$(LIB_SRC)/finiteVolume/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude

LIB_LIBS = \
    -lincompressibleTurbulenceModel \
    -lfiniteVolume \
    -lmeshTools

And the controlDict file of the case:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

libs ("mylibincompressibleRASModels.so");

application     pisoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.005;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;

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

finally, this is the error during running the case:

--> FOAM FATAL IO ERROR: 
Unknown patchField type epsilonDurbinWallFunction for patch type wall

Valid patchField types are :

82
(
advective
alphatJayatillekeWallFunction
atmBoundaryLayerInletEpsilon
calculated
codedFixedValue
codedMixed
cyclic
cyclicACMI
cyclicAMI
cyclicSlip
directionMixed
empty
epsilonLowReWallFunction
epsilonWallFunction
externalCoupled
fWallFunction
fan
fanPressure
fixedFluxPressure
fixedGradient
fixedInternalValue
fixedJump
fixedJumpAMI
fixedMean
fixedPressureCompressibleDensity
fixedValue
freestream
freestreamPressure
inletOutlet
inletOutletTotalTemperature
kLowReWallFunction
kqRWallFunction
mapped
mappedField
mappedFixedInternalValue
mappedFixedPushedInternalValue
mixed
nonuniformTransformCyclic
nutLowReWallFunction
nutTabulatedWallFunction
nutURoughWallFunction
nutUSpaldingWallFunction
nutUWallFunction
nutkAtmRoughWallFunction
nutkRoughWallFunction
nutkWallFunction
omegaWallFunction
oscillatingFixedValue
outletInlet
outletMappedUniformInlet
partialSlip
phaseHydrostaticPressure
prghPressure
processor
processorCyclic
rotatingTotalPressure
sliced
slip
symmetry
symmetryPlane
syringePressure
timeVaryingMappedFixedValue
totalPressure
totalTemperature
turbulentHeatFluxTemperature
turbulentInlet
turbulentIntensityKineticEnergyInlet
turbulentMixingLengthDissipationRateInlet
turbulentMixingLengthFrequencyInlet
uniformDensityHydrostaticPressure
uniformFixedGradient
uniformFixedValue
uniformInletOutlet
uniformJump
uniformJumpAMI
uniformTotalPressure
v2WallFunction
variableHeightFlowRate
waveSurfacePressure
waveTransmissive
wedge
zeroGradient
)


file: /home/julian/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough/0/epsilon.boundaryField.movingWall from line 26 to line 28.

    From function fvPatchField<Type>::New(const fvPatch&, const DimensionedField<Type, volMesh>&, const dictionary&)
    in file /home/opencfd/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/fvPatchFieldNew.C at line 143.

FOAM exiting

Thanks very much,
Julian

[alexeym]

Hi,

source code of the wall function will be quite helpful. Right now I can try to guess that you've forgotten to change wallFunction class name when registering it with addToRunTimeSelectionTable macro.

[juhuettn]

Hi Alexey, thanks for your reply!

The .C file:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     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 3 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, see <http://www.gnu.org/licenses/>.

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

#include "epsilonDurbinWallFunctionFvPatchScalarField.H"
#include "incompressible/turbulenceModel/turbulenceModel.H"
#include "fvPatchFieldMapper.H"
#include "fvMatrix.H"
#include "volFields.H"
#include "wallFvPatch.H"
#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
namespace incompressible
{

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

void epsilonDurbinWallFunctionFvPatchScalarField::checkType()
{
    if (!isA<wallFvPatch>(patch()))
    {
        FatalErrorIn("epsilonDurbinWallFunctionFvPatchScalarField::checkType()")
            << "Invalid wall function specification" << nl
            << "    Patch type for patch " << patch().name()
            << " must be wall" << nl
            << "    Current patch type is " << patch().type() << nl << endl
            << abort(FatalError);
    }
}


void epsilonDurbinWallFunctionFvPatchScalarField::writeLocalEntries(Ostream& os) const
{
    os.writeKeyword("Cmu") << Cmu_ << token::END_STATEMENT << nl;
    os.writeKeyword("kappa") << kappa_ << token::END_STATEMENT << nl;
    os.writeKeyword("E") << E_ << token::END_STATEMENT << nl;
}


void epsilonDurbinWallFunctionFvPatchScalarField::setMaster()
{
    if (master_ != -1)
    {
        return;
    }

    const volScalarField& epsilon =
        static_cast<const volScalarField&>(this->dimensionedInternalField());

    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();

    label master = -1;
    forAll(bf, patchI)
    {
        if (isA<epsilonDurbinWallFunctionFvPatchScalarField>(bf[patchI]))
        {
            epsilonDurbinWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);

            if (master == -1)
            {
                master = patchI;
            }

            epf.master() = master;
        }
    }
}


void epsilonDurbinWallFunctionFvPatchScalarField::createAveragingWeights()
{
    if (initialised_)
    {
        return;
    }

    const volScalarField& epsilon =
        static_cast<const volScalarField&>(this->dimensionedInternalField());

    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();

    const fvMesh& mesh = epsilon.mesh();

    volScalarField weights
    (
        IOobject
        (
            "weights",
            mesh.time().timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false // do not register
        ),
        mesh,
        dimensionedScalar("zero", dimless, 0.0)
    );

    DynamicList<label> epsilonPatches(bf.size());
    forAll(bf, patchI)
    {
        if (isA<epsilonDurbinWallFunctionFvPatchScalarField>(bf[patchI]))
        {
            epsilonPatches.append(patchI);

            const labelUList& faceCells = bf[patchI].patch().faceCells();
            forAll(faceCells, i)
            {
                label cellI = faceCells[i];
                weights[cellI]++;
            }
        }
    }

    cornerWeights_.setSize(bf.size());
    forAll(epsilonPatches, i)
    {
        label patchI = epsilonPatches[i];
        const fvPatchScalarField& wf = weights.boundaryField()[patchI];
        cornerWeights_[patchI] = 1.0/wf.patchInternalField();
    }

    G_.setSize(dimensionedInternalField().size(), 0.0);
    epsilon_.setSize(dimensionedInternalField().size(), 0.0);

    initialised_ = true;
}


epsilonDurbinWallFunctionFvPatchScalarField&
epsilonDurbinWallFunctionFvPatchScalarField::epsilonPatch(const label patchI)
{
    const volScalarField& epsilon =
        static_cast<const volScalarField&>(this->dimensionedInternalField());

    const volScalarField::GeometricBoundaryField& bf = epsilon.boundaryField();

    const epsilonDurbinWallFunctionFvPatchScalarField& epf =
        refCast<const epsilonDurbinWallFunctionFvPatchScalarField>(bf[patchI]);

    return const_cast<epsilonDurbinWallFunctionFvPatchScalarField&>(epf);
}


void epsilonDurbinWallFunctionFvPatchScalarField::calculateTurbulenceFields
(
    const turbulenceModel& turbulence,
    scalarField& G0,
    scalarField& epsilon0
)
{
    // accumulate all of the G and epsilon contributions
    forAll(cornerWeights_, patchI)
    {
        if (!cornerWeights_[patchI].empty())
        {
            epsilonDurbinWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);

            const List<scalar>& w = cornerWeights_[patchI];

            epf.calculate(turbulence, w, epf.patch(), G0, epsilon0);
        }
    }

    // apply zero-gradient condition for epsilon
    forAll(cornerWeights_, patchI)
    {
        if (!cornerWeights_[patchI].empty())
        {
            epsilonDurbinWallFunctionFvPatchScalarField& epf = epsilonPatch(patchI);

            epf == scalarField(epsilon0, epf.patch().faceCells());
        }
    }
}


void epsilonDurbinWallFunctionFvPatchScalarField::calculate
(
    const turbulenceModel& turbulence,
    const List<scalar>& cornerWeights,
    const fvPatch& patch,
    scalarField& G,
    scalarField& epsilon
)
{
    const label patchI = patch.index();

    const scalarField& y = turbulence.y()[patchI];

    const scalar Cmu25 = pow025(Cmu_);

    const tmp<volScalarField> tk = turbulence.k();
    const volScalarField& k = tk();

    const tmp<volScalarField> tnu = turbulence.nu();
    const scalarField& nuw = tnu().boundaryField()[patchI];

    const tmp<volScalarField> tnut = turbulence.nut();
    const volScalarField& nut = tnut();
    const scalarField& nutw = nut.boundaryField()[patchI];

    const fvPatchVectorField& Uw = turbulence.U().boundaryField()[patchI];

    const scalarField magGradUw(mag(Uw.snGrad()));
    const scalarField magUp(mag(Uw.patchInternalField() - Uw));
    
    tmp<scalarField> tfmu = calcfmu(magUp);    //Connection to calcfmu
    scalarField& fmu = tfmu();

    // Set epsilon and G
    
    //decide what epsilon to calculate
    forAll(nutw, faceI) //calculate epsilon, wallfunction is takeing epsilon from last time step, not this one!
    {
        if (fmu[faceI] <= 0.95) //chose epsilon innen
        {     
            
            tmp<scalarField> tyeff = calcyeff(magUp);    //Connection to calcYeff
            scalarField& yeff = tyeff();
            
            label cellI = patch.faceCells()[faceI];

            scalar w = cornerWeights[faceI];
                    
            scalar Rey = sqrt(k[cellI])*y[faceI]/nuw[faceI];
            epsilon[cellI] = pow(k[cellI],1.5)/(2.5*yeff[cellI]*(1-exp(Rey/5)));

            G[cellI] +=
                w
                *(nutw[faceI] + nuw[faceI])
                *magGradUw[faceI]
                *Cmu25*sqrt(k[cellI])
                /(kappa_*y[faceI]);
        }
    
        else //chose epsilon außen (keep epsilon as calculated before)
        {
            label cellI = patch.faceCells()[faceI];

            scalar w = cornerWeights[faceI];
            
            epsilon[cellI] += 0.0;
            
            G[cellI] +=
                w
                *(nutw[faceI] + nuw[faceI])
                *magGradUw[faceI]
                *Cmu25*sqrt(k[cellI])
                /(kappa_*y[faceI]);
        }
    }
}

tmp<scalarField> epsilonDurbinWallFunctionFvPatchScalarField::calcfmu

(
    const scalarField& magUp
) const
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();
    const scalar Cmu25 = pow025(Cmu_);
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];

    // The flow velocity at the cell centre
    const scalarField magUpfmu(mag(Uw.patchInternalField() - Uw));

    tmp<scalarField> tfmu(new scalarField(patch().size(), 0.0));
    scalarField& fmu = tfmu();

    forAll(fmu, facei)
       {
            label faceCelli = patch().faceCells()[facei];
            scalar uPlus = Cmu25*sqrt(k[facei]);
               scalar KsPlus = uPlus*Ks_[facei]/nuw[facei];
            scalar Amu = max(1,62.5*(1-KsPlus/90));        // (5.20)
            scalar Rey = sqrt(k[faceCelli])*y[facei]/nuw[facei];    //(5.19)

            fmu[facei] = 1-exp(-(Rey/Amu));    //(5.18)
        }
    
    return tfmu;
    
}

tmp<scalarField> epsilonDurbinWallFunctionFvPatchScalarField::calcyeff
(
    const scalarField& magUp
) const
{
    const label patchi = patch().index();

    const turbulenceModel& turbModel =
        db().lookupObject<turbulenceModel>("turbulenceModel");
    const scalarField& y = turbModel.y()[patchi];
    const tmp<volScalarField> tnu = turbModel.nu();
    const volScalarField& nu = tnu();
    const scalarField& nuw = nu.boundaryField()[patchi];
    const tmp<volScalarField> tk = turbModel.k();
    const volScalarField& k = tk();
    const scalar Cmu25 = pow025(Cmu_);
    const fvPatchVectorField& Uw = turbModel.U().boundaryField()[patchi];

    const scalarField magUpyeff(mag(Uw.patchInternalField() - Uw));

    tmp<scalarField> tyeff(new scalarField(patch().size(), 0.0));
    scalarField& yeff = tyeff();
    
forAll(yeff, facei)
    {
    
            const scalar magUpara = magUp[facei];
            //const scalarField& y0Plus;
            label faceCelli = patch().faceCells()[facei];
                            
            tmp<scalarField> tyZeroPlus(new scalarField(patch().size(), 0.0));
            scalarField& yZeroPlus = tyZeroPlus();

            scalar uPlus = Cmu25*sqrt(k[faceCelli]);
               scalar KsPlus = uPlus*Ks_[facei]/nuw[facei];

                    if (KsPlus > 90)
                    {
                        yZeroPlus = 0.035*KsPlus-1.548;    //(5.24, part 2)
                    }
                    else //if (KsPlus <= 90)
                    {
                        yZeroPlus = (-1.16+0.85*sqrt(KsPlus)-0.071*KsPlus)/(1-0.16*sqr(KsPlus)+0.01*KsPlus); //(5.24, part 1)
                    }

            yeff[facei] = y[facei] + yZeroPlus[facei]*nuw[facei]*uPlus/magUpara;    //Formula developed from formula on p.124 bottom,
                                                                            //utau=u/uplus and yeff=y+y0
    }
    
    return tyeff;
}


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

epsilonDurbinWallFunctionFvPatchScalarField::epsilonDurbinWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedValueFvPatchField<scalar>(p, iF),
    Cmu_(0.09),
    kappa_(0.41),
    E_(9.8),
    G_(),
    epsilon_(),
    initialised_(false),
    master_(-1),
    cornerWeights_()
{
    checkType();
}


epsilonDurbinWallFunctionFvPatchScalarField::epsilonDurbinWallFunctionFvPatchScalarField
(
    const epsilonDurbinWallFunctionFvPatchScalarField& ptf,
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const fvPatchFieldMapper& mapper
)
:
    fixedValueFvPatchField<scalar>(ptf, p, iF, mapper),
    Cmu_(ptf.Cmu_),
    kappa_(ptf.kappa_),
    E_(ptf.E_),
    G_(),
    epsilon_(),
    initialised_(false),
    master_(-1),
    cornerWeights_()
{
    checkType();
}


epsilonDurbinWallFunctionFvPatchScalarField::epsilonDurbinWallFunctionFvPatchScalarField
(
    const fvPatch& p,
    const DimensionedField<scalar, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchField<scalar>(p, iF, dict),
    Cmu_(dict.lookupOrDefault<scalar>("Cmu", 0.09)),
    kappa_(dict.lookupOrDefault<scalar>("kappa", 0.41)),
    E_(dict.lookupOrDefault<scalar>("E", 9.8)),
    G_(),
    epsilon_(),
    initialised_(false),
    master_(-1),
    cornerWeights_()
{
    checkType();

    // apply zero-gradient condition on start-up
    this->operator==(patchInternalField());
}


epsilonDurbinWallFunctionFvPatchScalarField::epsilonDurbinWallFunctionFvPatchScalarField
(
    const epsilonDurbinWallFunctionFvPatchScalarField& ewfpsf
)
:
    fixedValueFvPatchField<scalar>(ewfpsf),
    Cmu_(ewfpsf.Cmu_),
    kappa_(ewfpsf.kappa_),
    E_(ewfpsf.E_),
    G_(),
    epsilon_(),
    initialised_(false),
    master_(-1),
    cornerWeights_()
{
    checkType();
}


epsilonDurbinWallFunctionFvPatchScalarField::epsilonDurbinWallFunctionFvPatchScalarField
(
    const epsilonDurbinWallFunctionFvPatchScalarField& ewfpsf,
    const DimensionedField<scalar, volMesh>& iF
)
:
    fixedValueFvPatchField<scalar>(ewfpsf, iF),
    Cmu_(ewfpsf.Cmu_),
    kappa_(ewfpsf.kappa_),
    E_(ewfpsf.E_),
    G_(),
    epsilon_(),
    initialised_(false),
    master_(-1),
    cornerWeights_()
{
    checkType();
}


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

scalarField& epsilonDurbinWallFunctionFvPatchScalarField::G(bool init)
{
    if (patch().index() == master_)
    {
        if (init)
        {
            G_ = 0.0;
        }

        return G_;
    }

    return epsilonPatch(master_).G();
}


scalarField& epsilonDurbinWallFunctionFvPatchScalarField::epsilon(bool init)
{
    if (patch().index() == master_)
    {
        if (init)
        {
            epsilon_ = 0.0;
        }

        return epsilon_;
    }

    return epsilonPatch(master_).epsilon(init);
}


void epsilonDurbinWallFunctionFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const turbulenceModel& turbulence =
        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);

    setMaster();

    if (patch().index() == master_)
    {
        createAveragingWeights();
        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
    }

    const scalarField& G0 = this->G();
    const scalarField& epsilon0 = this->epsilon();

    typedef DimensionedField<scalar, volMesh> FieldType;

    FieldType& G =
        const_cast<FieldType&>
        (
            db().lookupObject<FieldType>(turbulence.GName())
        );

    FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());

    forAll(*this, faceI)
    {
        label cellI = patch().faceCells()[faceI];

        G[cellI] = G0[cellI];
        epsilon[cellI] = epsilon0[cellI];
    }

    fvPatchField<scalar>::updateCoeffs();
}


void epsilonDurbinWallFunctionFvPatchScalarField::updateCoeffs
(
    const scalarField& weights
)
{
    if (updated())
    {
        return;
    }

    const turbulenceModel& turbulence =
        db().lookupObject<turbulenceModel>(turbulenceModel::typeName);

    setMaster();

    if (patch().index() == master_)
    {
        createAveragingWeights();
        calculateTurbulenceFields(turbulence, G(true), epsilon(true));
    }

    const scalarField& G0 = this->G();
    const scalarField& epsilon0 = this->epsilon();

    typedef DimensionedField<scalar, volMesh> FieldType;

    FieldType& G =
        const_cast<FieldType&>
        (
            db().lookupObject<FieldType>(turbulence.GName())
        );

    FieldType& epsilon = const_cast<FieldType&>(dimensionedInternalField());

    scalarField& epsilonf = *this;

    // only set the values if the weights are < 1 - tolerance
    forAll(weights, faceI)
    {
        scalar w = weights[faceI];

        if (w < 1.0 - 1e-6)
        {
            label cellI = patch().faceCells()[faceI];

            G[cellI] = w*G[cellI] + (1.0 - w)*G0[cellI];
            epsilon[cellI] = w*epsilon[cellI] + (1.0 - w)*epsilon0[cellI];
            epsilonf[faceI] = epsilon[cellI];
        }
    }

    fvPatchField<scalar>::updateCoeffs();
}


void epsilonDurbinWallFunctionFvPatchScalarField::manipulateMatrix
(
    fvMatrix<scalar>& matrix
)
{
    if (manipulatedMatrix())
    {
        return;
    }

    matrix.setValues(patch().faceCells(), patchInternalField());

    fvPatchField<scalar>::manipulateMatrix(matrix);
}


void epsilonDurbinWallFunctionFvPatchScalarField::manipulateMatrix
(
    fvMatrix<scalar>& matrix,
    const Field<scalar>& weights
)
{
    if (manipulatedMatrix())
    {
        return;
    }

    // filter weights so that we only apply the constraint where the
    // weight > SMALL
    DynamicList<label> constraintCells(weights.size());
    DynamicList<scalar> constraintEpsilon(weights.size());
    const labelUList& faceCells = patch().faceCells();

    const DimensionedField<scalar, volMesh>& epsilon
        = dimensionedInternalField();

    label nConstrainedCells = 0;


    forAll(weights, faceI)
    {
        // only set the values if the weights are < 1 - tolerance
        if (weights[faceI] < (1.0 - 1e-6))
        {
            nConstrainedCells++;

            label cellI = faceCells[faceI];

            constraintCells.append(cellI);
            constraintEpsilon.append(epsilon[cellI]);
        }
    }

    if (debug)
    {
        Pout<< "Patch: " << patch().name()
            << ": number of constrained cells = " << nConstrainedCells
            << " out of " << patch().size()
            << endl;
    }

    matrix.setValues
    (
        constraintCells,
        scalarField(constraintEpsilon.xfer())
    );

    fvPatchField<scalar>::manipulateMatrix(matrix);
}


void epsilonDurbinWallFunctionFvPatchScalarField::write(Ostream& os) const
{
    fixedValueFvPatchField<scalar>::write(os);
    writeLocalEntries(os);
    writeEntry("value", os);
}


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

makePatchTypeField
(
    fvPatchScalarField,
    epsilonDurbinWallFunctionFvPatchScalarField
);

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

} // End namespace incompressible
} // End namespace Foam

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

[juhuettn]

And the .H file:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     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 3 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, see <http://www.gnu.org/licenses/>.

Class
    Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField

Group
    grpIcoWallFunctions

Description
    This boundary condition provides a turbulence dissipation wall function
    condition for high Reynolds number, turbulent flow cases.

    The condition can be applied to wall boundaries, whereby it
    - calculates \c epsilon and \c G
    - inserts near wall epsilon values directly into the epsilon equation
        to act as a constraint

    where

    \vartable
        epsilon | turblence dissipation field
        G       | turblence generation field
    \endvartable

    \heading Patch usage

    \table
        Property     | Description             | Required    | Default value
        Cmu          | model coefficient       | no          | 0.09
        kappa        | Von Karman constant     | no          | 0.41
        E            | model coefficient       | no          | 9.8
    \endtable

    Example of the boundary condition specification:
    \verbatim
    myPatch
    {
        type            epsilonDurbinWallFunction;
    }
    \endverbatim

SeeAlso
    Foam::fixedInternalValueFvPatchField

SourceFiles
    epsilonDurbinWallFunctionFvPatchScalarField.C

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

#ifndef epsilonDurbinWallFunctionFvPatchScalarField_H
#define epsilonDurbinWallFunctionFvPatchScalarField_H

#include "fixedValueFvPatchField.H"

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

namespace Foam
{
namespace incompressible
{

class turbulenceModel;

/*---------------------------------------------------------------------------*\
           Class epsilonDurbinWallFunctionFvPatchScalarField Declaration
\*---------------------------------------------------------------------------*/

class epsilonDurbinWallFunctionFvPatchScalarField
:
    public fixedValueFvPatchField<scalar>
{
protected:

    // Protected data

        //- Cmu coefficient
        scalar Cmu_;

        //- Von Karman constant
        scalar kappa_;

        //- E coefficient
        scalar E_;

        //- Local copy of turbulence G field
        scalarField G_;

        //- Local copy of turbulence epsilon field
        scalarField epsilon_;

        //- Initialised flag
        bool initialised_;

        //- Master patch ID
        label master_;

        //- List of averaging corner weights
        List<List<scalar> > cornerWeights_;
        
        //- Roughness height
        scalarField Ks_;


    // Protected Member Functions

        //- Check the type of the patch
        virtual void checkType();

        //- Write local wall function variables
        virtual void writeLocalEntries(Ostream&) const;

        //- Set the master patch - master is responsible for updating all
        //  wall function patches
        virtual void setMaster();

        //- Create the averaging weights for cells which are bounded by
        //  multiple wall function faces
        virtual void createAveragingWeights();

        //- Helper function to return non-const access to an epsilon patch
        virtual epsilonDurbinWallFunctionFvPatchScalarField& epsilonPatch
        (
            const label patchI
        );

        //- Main driver to calculate the turbulence fields
        virtual void calculateTurbulenceFields
        (
            const turbulenceModel& turbulence,
            scalarField& G0,
            scalarField& epsilon0
        );

        //- Calculate the epsilon and G
        virtual void calculate
        (
            const turbulenceModel& turbulence,
            const List<scalar>& cornerWeights,
            const fvPatch& patch,
            scalarField& G,
            scalarField& epsilon
        );

        //- Return non-const access to the master patch ID
        virtual label& master()
        {
            return master_;
        }
        
        //- Calculate fmu
        virtual tmp<scalarField> calcfmu(const scalarField& magUp) const;
        
        //- Calculate yeff
        virtual tmp<scalarField> calcyeff(const scalarField& magUp) const;


public:

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


    // Constructors

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

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

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

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

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

        //- Construct as copy setting internal field reference
        epsilonDurbinWallFunctionFvPatchScalarField
        (
            const epsilonDurbinWallFunctionFvPatchScalarField&,
            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 epsilonDurbinWallFunctionFvPatchScalarField(*this, iF)
            );
        }

    //- Destructor
    virtual ~epsilonDurbinWallFunctionFvPatchScalarField()
    {}


    // Member functions

        // Access

            //- Return non-const access to the master's G field
            scalarField& G(bool init = false);

            //- Return non-const access to the master's epsilon field
            scalarField& epsilon(bool init = false);


        // Evaluation functions

            //- Update the coefficients associated with the patch field
            virtual void updateCoeffs();

            //- Update the coefficients associated with the patch field
            virtual void updateCoeffs(const scalarField& weights);

            //- Manipulate matrix
            virtual void manipulateMatrix(fvMatrix<scalar>& matrix);

            //- Manipulate matrix with given weights
            virtual void manipulateMatrix
            (
                fvMatrix<scalar>& matrix,
                const scalarField& weights
            );


        // I-O

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


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

} // End namespace incompressible
} // End namespace Foam

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

#endif

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

[alexeym]

Hi,

I wasn't able to reproduce your behavior. I've copied the files you've provided, compiled library, added lib("...") in controlDict, changed epsilon BCs to epsilonDurbinWallFunction and got this output (the case is pitzDaily from incompressible/pimpleFoam tutorial):

smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 8.07931e-06, No Iterations 4
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 4.8644e-06, No Iterations 5
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.00925834, No Iterations 12
time step continuity errors : sum local = 1.80001e-05, global = -1.4582e-06, cumulative = -1.4582e-06
GAMG:  Solving for p, Initial residual = 0.0217289, Final residual = 9.10856e-08, No Iterations 20
time step continuity errors : sum local = 1.80266e-08, global = 1.31557e-09, cumulative = -1.45689e-06
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigSegv::sigHandler(int) at ??:?
#2   in "/lib64/libc.so.6"
#3  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::calcfmu(Foam::Field<double> const&) const at ??:?
#4  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::calculate(Foam::incompressible::turbulenceModel const&, Foam::List<double> const&, Foam::fvPatch const&, Foam::Field<double>&, Foam::Field<double>&) at ??:?
#5  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::calculateTurbulenceFields(Foam::incompressible::turbulenceModel const&, Foam::Field<double>&, Foam::Field<double>&) at ??:?
#6  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::updateCoeffs() at ??:?
#7  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::updateCoeffs() at ??:?
#8  Foam::incompressible::RASModels::kEpsilon::correct() at ??:?
#9  
 at ??:?
#10  __libc_start_main in "/lib64/libc.so.6"
#11  
 at ??:?
Erreur de segmentation (core dumped)

I.e. run time selection mechanism was able to create wall function class.

[juhuettn]

Hi Alexey.

thank you for trying it out. I also did all the steps again for the pitzDaily case and still getting the same error. Is there any other issue during the compiling process which can be the reason for my compiled but not recognized wall function? E.g. does it matter in which folder I copy my custom wall function or something like that?

[alexeym]

Hi,

Originally Posted by juhuettn

E.g. does it matter in which folder I copy my custom wall function or something like that?

Not sure, my placement of the sources was rather arbitrary. Finally compiled library matters, not location of the source files.

There were two differences in my setup:

1. I've used pimpleFoam instead of pisoFoam (though it seems, you were able to reproduce error with pitzDaily and pimpleFoam)
2. I've placed lib ("...") in the end of controlDict.

Can you post whole log-file of the run? Maybe there were errors during load of mylibincompressibleRASModels.so.

[juhuettn]

Hi Alexey,

I also tried the pitzDaily case with putting the lib ("...") in the end of the file, but no difference. The error message I posted first is already the complete log file of the case.

Might there be an issue with my user library connection to the main program? The library is built and I can also find it in the /user-2.3.0/platforms/linux64GccDPOpt/lib folder, but somehow it is not used during running a case. Is there something general I need to do before using user libraries? Because, it is the first library I am trying to link.

Thanks

[alexeym]

Hi,

let me learn a little bit more about your installation (maybe it's the problem with the build or something similar), i.e. what's your Linux distributive? did you install OpenFOAM using package manager or did you build it from sources?

The error message I posted first is already the complete log file of the case.

I'm not quite sure, are the message lacks standard OpenFOAM solver output header and further messages about selection of the models. Also there could be a message about missing library (due to typo in the name or something else).

P.S. I hope it was clear that by libs("...") I meant libs ("mylibincompressibleRASModels.so")

[juhuettn]

Hi Alexey,

sorry I ran the case with the additional command ">>log &". so I did not see the complete log directly. But there is not much more in the log than the error message, have a look:

julian@julian-ubuntu:~/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough$ pisoFoam
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.3.0-f5222ca19ce6
Exec   : pisoFoam
Date   : Oct 13 2014
Time   : 16:50:23
Host   : "julian-ubuntu"
PID    : 2250
Case   : /home/julian/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RASModel
Selecting RAS turbulence model kEpsilon


--> FOAM FATAL IO ERROR: 
Unknown patchField type epsilonDurbinWallFunction for patch type wall

Valid patchField types are :

82
(
advective
alphatJayatillekeWallFunction
atmBoundaryLayerInletEpsilon
calculated
codedFixedValue
codedMixed
cyclic
cyclicACMI
cyclicAMI
cyclicSlip
directionMixed
empty
epsilonLowReWallFunction
epsilonWallFunction
externalCoupled
fWallFunction
fan
fanPressure
fixedFluxPressure
fixedGradient
fixedInternalValue
fixedJump
fixedJumpAMI
fixedMean
fixedPressureCompressibleDensity
fixedValue
freestream
freestreamPressure
inletOutlet
inletOutletTotalTemperature
kLowReWallFunction
kqRWallFunction
mapped
mappedField
mappedFixedInternalValue
mappedFixedPushedInternalValue
mixed
nonuniformTransformCyclic
nutLowReWallFunction
nutTabulatedWallFunction
nutURoughWallFunction
nutUSpaldingWallFunction
nutUWallFunction
nutkAtmRoughWallFunction
nutkRoughWallFunction
nutkWallFunction
omegaWallFunction
oscillatingFixedValue
outletInlet
outletMappedUniformInlet
partialSlip
phaseHydrostaticPressure
prghPressure
processor
processorCyclic
rotatingTotalPressure
sliced
slip
symmetry
symmetryPlane
syringePressure
timeVaryingMappedFixedValue
totalPressure
totalTemperature
turbulentHeatFluxTemperature
turbulentInlet
turbulentIntensityKineticEnergyInlet
turbulentMixingLengthDissipationRateInlet
turbulentMixingLengthFrequencyInlet
uniformDensityHydrostaticPressure
uniformFixedGradient
uniformFixedValue
uniformInletOutlet
uniformJump
uniformJumpAMI
uniformTotalPressure
v2WallFunction
variableHeightFlowRate
waveSurfacePressure
waveTransmissive
wedge
zeroGradient
)


file: /home/julian/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough/0/epsilon.boundaryField.movingWall from line 26 to line 28.

    From function fvPatchField<Type>::New(const fvPatch&, const DimensionedField<Type, volMesh>&, const dictionary&)
    in file /home/opencfd/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/fvPatchFieldNew.C at line 143.

FOAM exiting

Concerning the libs("...") command, I was getting that right.

[alexeym]

Well, let's try to look at your case files could you attach archive?

[juhuettn]

Hi Alexey,

I attached the archive of my case. I also found something else in another thread. People having problems with their libraries, sometimes needed to add the line

export LD_LIBRARY_PATH=$FOAM_USER_LIBBIN

to their /.bashrc file. I already tried this with no success. Maybe you can see if there is something wrong with my library paths from the output of "echo $LD_LIBRARY_PATH":

julian@julian-ubuntu:~$ echo $LD_LIBRARY_PATH
/opt/ThirdParty-2.3.0/platforms/linux64Gcc/gperftools-svn/lib:/opt/paraviewopenfoam410/lib/paraview-4.1:/opt/openfoam230/platforms/linux64GccDPOpt/lib/openmpi-system:/opt/ThirdParty-2.3.0/platforms/linux64GccDPOpt/lib/openmpi-system:/usr/lib/openmpi/lib:/home/julian/OpenFOAM/julian-2.3.0/platforms/linux64GccDPOpt/lib:/opt/site/2.3.0/platforms/linux64GccDPOpt/lib:/opt/openfoam230/platforms/linux64GccDPOpt/lib:/opt/ThirdParty-2.3.0/platforms/linux64GccDPOpt/lib:/opt/openfoam230/platforms/linux64GccDPOpt/lib/dummy

Regards,
Julian

[alexeym]

Hi,

it's rather funny but ... controlDict in the attached archive has a line

lib ("mylibincompressibleRASModels.so");

while it should be

libs ("mylibincompressibleRASModels.so");

[juhuettn]

Hi Alexey,

sorry I played around with that line to try "lib" and "libs" and I forgot to change it back. But also when using libs ("..."), I still get the same error.

Maybe I do something wrong during compiling? I copy the complete RAS model from the source folder and then compile the wall functions from the Make folder inside the RAS folder is that correct? Do I maybe need to recompile the whole RAS model?

[alexeym]

Hi,

from your PATHs I've figured out that it is Ubuntu, what version? Did you install binary version of OpenFOAM or built it from sources (though as it's in /opt folder I guess it's binary installation)?

I've tested your source files on OS X and Linux (unfortunately neither was binary installation). In both cases I've just created a folder (let's call it cfdonline-src) with 2 folders inside: Make and epsilonDurbinWallFunction, put source files into epsilonDurbinWallFunction folder (two files you've posted), in Make folder I've created files and options files, then I just ran wmake in cfdonline-src folder. And that is all. On Linux I've tested library with pimpleFoam and pitzDaily tutorial case (and it has run without any complains), on OS X I have tested the archive you have posted and it has run after I had changed lib to libs in controlDict.

I can check if the problem persists on Ubuntu with binary OpenFOAM installation as soon as you post version number you're using (maybe not so soon as I'll need to create VM with Ubuntu ).

[juhuettn]

Hi Alexey,

Yes I am running binary OpenFOAM installation on Ubuntu. I just figured out follwoing issue:

I also compile another wallFunction "nutDurbinWallFunction" together with the "epsilonDurbinWallfunction". Although I do not use the "nutDurbinWallFunction" in the moment, it somehow seems to make the library not working.

Maybe you can try it as well. Compile the attached files with this Make/files file from my first post:

/* Wall functions */
wallFunctions = derivedFvPatchFields/wallFunctions

nutWallFunctions = $(wallFunctions)/nutWallFunctions
$(nutWallFunctions)/nutDurbinWallFunction/nutDurbinWallFunctionFvPatchScalarField.C

epsilonWallFunctions = $(wallFunctions)/epsilonWallFunctions
$(epsilonWallFunctions)/epsilonDurbinWallFunction/epsilonDurbinWallFunctionFvPatchScalarField.C

LIB = $(FOAM_USER_LIBBIN)/mylibincompressibleRASModels

]I just recognised, that I get an additional Warning in the beginning of the log when running the case with "libs" instead of "lib" in the controlDict (first few lines):

Create time

--> FOAM Warning : 
    From function dlOpen(const fileName&, const bool)
    in file POSIX.C at line 1179
    dlopen error : /home/julian/OpenFOAM/julian-2.3.0/platforms/linux64GccDPOpt/lib/mylibincompressibleRASModels.so: undefined symbol: _ZTIN4Foam14incompressible39nutDurbinWallFunctionFvPatchScalarFieldE
--> FOAM Warning : 
    From function dlLibraryTable::open(const fileName&, const bool)
    in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 99
    could not load "mylibincompressibleRASModels.so"
Create mesh for time = 0

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RASModel
Selecting RAS turbulence model kEpsilon


--> FOAM FATAL IO ERROR: 
Unknown patchField type epsilonDurbinWallFunction for patch type wall

Valid patchField types are :

82
(
advective
alphatJayatillekeWallFunction
atmBoundaryLayerInletEpsilon
calculated
codedFixedValue
codedMixed
cyclic
cyclicACMI
cyclicAMI
cyclicSlip
directionMixed
empty
epsilonLowReWallFunction
epsilonWallFunction
externalCoupled
fWallFunction
fan
fanPressure
fixedFluxPressure
fixedGradient
fixedInternalValue
fixedJump
fixedJumpAMI
fixedMean
fixedPressureCompressibleDensity
fixedValue
freestream
freestreamPressure
inletOutlet
inletOutletTotalTemperature
kLowReWallFunction
kqRWallFunction
mapped
mappedField
mappedFixedInternalValue
mappedFixedPushedInternalValue
mixed
nonuniformTransformCyclic
nutLowReWallFunction
nutTabulatedWallFunction
nutURoughWallFunction
nutUSpaldingWallFunction
nutUWallFunction
nutkAtmRoughWallFunction
nutkRoughWallFunction
nutkWallFunction
omegaWallFunction
oscillatingFixedValue
outletInlet
outletMappedUniformInlet
partialSlip
phaseHydrostaticPressure
prghPressure
processor
processorCyclic
rotatingTotalPressure
sliced
slip
symmetry
symmetryPlane
syringePressure
timeVaryingMappedFixedValue
totalPressure
totalTemperature
turbulentHeatFluxTemperature
turbulentInlet
turbulentIntensityKineticEnergyInlet
turbulentMixingLengthDissipationRateInlet
turbulentMixingLengthFrequencyInlet
uniformDensityHydrostaticPressure
uniformFixedGradient
uniformFixedValue
uniformInletOutlet
uniformJump
uniformJumpAMI
uniformTotalPressure
v2WallFunction
variableHeightFlowRate
waveSurfacePressure
waveTransmissive
wedge
zeroGradient
)


file: /home/julian/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough/0/epsilon.boundaryField.movingWall from line 26 to line 28.

    From function fvPatchField<Type>::New(const fvPatch&, const DimensionedField<Type, volMesh>&, const dictionary&)
    in file /home/opencfd/OpenFOAM/OpenFOAM-2.3.0/src/finiteVolume/lnInclude/fvPatchFieldNew.C at line 143.

FOAM exiting

So I guess something is wrong with this other wallFunction "nutDurbinWallFunction". I get no error during compiling and I can also not find anything in the file itself. But without compiling "nutDurbinWallFuntcion", the library works and I can use the "epsilonDurbinWallfunction".

Thanks for your help!

[alexeym]

Hi,

you've got two unimplemented methods in nutDurbinWallFunctionFvPatchScalarField.C:

calcYPlus and writeLocalEntries so nutDurbinWallFunctionFvPatchScalarField can't be instantiated. After addition of empty implementations of these methods dlopen error has gone.

Also I not quite get the idea behind:

virtual tmp<scalarField> calcfmu(const scalarField& magUp) const;

if you'd like to recalculate member field fmu_ why method has const specifier?

the same thing with (i.e. I do not understand the logic behind)

	//- fmu
		scalarField fmu_;
...            
            //- Calculate and return fmu
            virtual tmp<scalarField> fmu() const;

why fmu_ is public?

[juhuettn]

The idea is to make fmu_ also available for other objects, such as epsilonDurbinWallFunction. If you look at the if loop in "epsilonDurbinWallFunction", it uses fmu as value to decide which calculation method to choose:

    forAll(nutw, faceI)
    {
        if (fmu[faceI] <= 0.95) //chose epsilon innen
        {     
            
            tmp<scalarField> tyeff = calcyeff(magUp);
            scalarField& yeff = tyeff();
            
            label cellI = patch.faceCells()[faceI];

            scalar w = cornerWeights[faceI];
                    
            scalar Rey = sqrt(k[cellI])*y[faceI]/nuw[faceI];
            epsilon[cellI] = pow(k[cellI],1.5)/(2.5*yeff[cellI]*(1-exp(Rey/5)));

            G[cellI] +=
                w
                *(nutw[faceI] + nuw[faceI])
                *magGradUw[faceI]
                *Cmu25*sqrt(k[cellI])
                /(kappa_*y[faceI]);
        }
    
        else //chose epsilon außen (keep epsilon as calculated before)
        {
            label cellI = patch.faceCells()[faceI];

            scalar w = cornerWeights[faceI];
            
            epsilon[cellI] += 0.0;
            
            G[cellI] +=
                w
                *(nutw[faceI] + nuw[faceI])
                *magGradUw[faceI]
                *Cmu25*sqrt(k[cellI])
                /(kappa_*y[faceI]);
        }
    }

Now, the fmu is still calculated in epsilonDurbinWallFunction as well, as I did not manage to "connect" the two epsilonDurbin and nutDurbin wall functions. Aim of course is to calculate it only once in "nutDurbinWallFunction". So I thought if I make that fmu public, it should be usable for other objects as well. I think this is a complete wrong method, but I just started working with object orientated C++ and I am actual mechanical engineer, so please excuse me :-).

Also the two unused methods were just a fault of mine. I tried to build up the .H file together with the .C step by step until it was compiling without error but did not see other things which could cause problems later on.

Now I can also run the case using the two wall functions without the dlopen error when following your hint, but I get multiple other errors when running a case. I now want to get the epsilonDurbinWallfunction running first, what do you think is the best method to find the errors and debug it? The error message is as follows:

julian@julian-ubuntu:~/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough$ pisoFoam
/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 2.3.0-f5222ca19ce6
Exec   : pisoFoam
Date   : Oct 15 2014
Time   : 13:39:20
Host   : "julian-ubuntu"
PID    : 6424
Case   : /home/julian/OpenFOAM/julian-2.3.0/run/tutorials/incompressible/pisoFoam/ras/cavityRough
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Disallowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type RASModel
Selecting RAS turbulence model kEpsilon
kEpsilonCoeffs
{
    Cmu             0.09;
    C1              1.44;
    C2              1.92;
    sigmaEps        1.3;
}


Starting time loop

Time = 5e-05

Courant Number mean: 0 max: 0
smoothSolver:  Solving for Ux, Initial residual = 1, Final residual = 1.14667e-07, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 1, Final residual = 6.98716e-08, No Iterations 3
GAMG:  Solving for p, Initial residual = 1, Final residual = 0.0405156, No Iterations 2
time step continuity errors : sum local = 4.99544e-09, global = 1.34417e-24, cumulative = 1.34417e-24
GAMG:  Solving for p, Initial residual = 0.0201159, Final residual = 7.52232e-07, No Iterations 7
time step continuity errors : sum local = 2.47471e-13, global = 3.46291e-25, cumulative = 1.69046e-24
#0  Foam::error::printStack(Foam::Ostream&) in "/opt/openfoam230/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#1  Foam::sigFpe::sigHandler(int) in "/opt/openfoam230/platforms/linux64GccDPOpt/lib/libOpenFOAM.so"
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::calculate(Foam::incompressible::turbulenceModel const&, Foam::List<double> const&, Foam::fvPatch const&, Foam::Field<double>&, Foam::Field<double>&) in "/home/julian/OpenFOAM/julian-2.3.0/platforms/linux64GccDPOpt/lib/mylibincompressibleRASModels.so"
#4  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::calculateTurbulenceFields(Foam::incompressible::turbulenceModel const&, Foam::Field<double>&, Foam::Field<double>&) in "/home/julian/OpenFOAM/julian-2.3.0/platforms/linux64GccDPOpt/lib/mylibincompressibleRASModels.so"
#5  Foam::incompressible::epsilonDurbinWallFunctionFvPatchScalarField::updateCoeffs() in "/home/julian/OpenFOAM/julian-2.3.0/platforms/linux64GccDPOpt/lib/mylibincompressibleRASModels.so"
#6  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::GeometricBoundaryField::updateCoeffs() in "/opt/openfoam230/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
#7  Foam::incompressible::RASModels::kEpsilon::correct() in "/opt/openfoam230/platforms/linux64GccDPOpt/lib/libincompressibleRASModels.so"
#8  
 in "/opt/openfoam230/platforms/linux64GccDPOpt/bin/pisoFoam"
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  
 in "/opt/openfoam230/platforms/linux64GccDPOpt/bin/pisoFoam"
Floating point exception (core dumped)

Maybe you can give me a tip what is a good way how to debug at this stage. Is there some library for looking up what the error messages mean or something similar? Because I cannot really get a lot of information out of it which is useful for me.

Thanks for your help so far and for your patience!

[alexeym]

Originally Posted by juhuettn

The idea is to make fmu_ also available for other objects, such as epsilonDurbinWallFunction. If you look at the if loop in "epsilonDurbinWallFunction", it uses fmu as value to decide which calculation method to choose:

But for some reason for Ks_ you've implemented accessor method (which returns constant reference for the field) but not for fmu_.

I think this is a complete wrong method, but I just started working with object orientated C++ and I am actual mechanical engineer, so please excuse me :-).

It depends if you plan to release the code, yeah, making the field public is not so bright idea, while if you need this code just for internal usage - whatever works; if you'd like the field to be public make it public, but still the method which recalculates the field can't have const specifier cause it changes fmu_ field of the object.

Also the two unused methods were just a fault of mine. I tried to build up the .H file together with the .C step by step until it was compiling without error but did not see other things which could cause problems later on.

I found the problem while converting nutUSpaldingWallFunction code into nutDurbinWallFunction. Library that I've built using new code complained explicitly about missing methods during dlopen, while library with old nutDurbinWallFunction just complained about missing constructor.

Now I can also run the case using the two wall functions without the dlopen error when following your hint, but I get multiple other errors when running a case. I now want to get the epsilonDurbinWallfunction running first, what do you think is the best method to find the errors and debug it?

Ideally, you should install (or build) debug version of OpenFOAM and execute solver under gdb or another debugger, set breakpoint in the method with error, examine variables and so on till bug is found and fixed. Debug version will print more information about the error (for example, line number where error happened), so maybe running solver under gdb won't be necessary.

Also you can try to use common sense. As the error is FPE usually it's division by zero or fractional power of zero or negative value. For example in calculate method of epsilonDurbinWallFunction you divide by nuw[faceI] which is basically nu taken from turbulence model, are you sure it's non-zero?

[juhuettn]

Thanks a lot for that hints, they are very useful for me! I will work through it and let you know how I am getting on. I think this is a good exercise for me to gain experience :-)...