[zordiack]

Hi, I'm currently programming a RAS version of kOmegaSSTSAS turbulence model and I ran into a problem when the cell volumes are needed in calculation. Cell volumes are needed in the calculation of high frequency eddy viscosity limiter Lvk2.

The LES model from current OpenFOAM-2.3.x has this:

tmp<volScalarField> kOmegaSSTSAS::Lvk2
(
    const volScalarField& S2
) const
{
    return max
    (
        kappa_*sqrt(S2)
       /(
            mag(fvc::laplacian(U()))
          + dimensionedScalar
            (
                "ROOTVSMALL",
                dimensionSet(0, -1 , -1, 0, 0, 0, 0),
                ROOTVSMALL
            )
        ),
        Cs_*delta()
    );
}

I'm quessing the delta() is implemented in LES models. The formulation for the limiter can be found for example in paper Y. Egorov and F. Menter - Development and Application of SST-SAS Turbulence Model in the DESIDER project (2008). This is what I've tried so far:

tmp<volScalarField> mykOmegaSSTSAS::Lvk
(
    const volScalarField& S2,
    const volScalarField& F1
) const

{
    return max
    (
        kappa_*sqrt(S2)
       /(
            mag(fvc::laplacian(U_))
          + dimensionedScalar
            (
                "ROOTVSMALL",
                dimensionSet(0, -1 , -1, 0, 0, 0, 0),
                ROOTVSMALL
            )
        ),
        Cs_*sqrt((kappa_*zeta2_)/((beta(F1)/Cmu_)-gamma(F1)))*pow(mesh_.V(),(1.0/3.0))
    );
}

I think the problem is that mesh_.V() is not a volScalarField. How should cell volumes be implemented in that equation? I'm not an experienced coder and I'm kind of clueless here

This is the error message:

mykOmegaSSTSAS/mykOmegaSSTSAS.C:89:12: error: no matching function for call to 'max'
    return max
           ^~~
/home/zordiack/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/label.H:295:1: note: candidate function not viable: no known conversion from 'tmp<GeometricField<scalar, fvPatchField, Foam::volMesh> >' to 'const char' for 1st argument
MAXMIN(char, char, char)
^
/home/zordiack/OpenFOAM/OpenFOAM-2.3.x/src/OpenFOAM/lnInclude/label.H:284:16: note: expanded from macro 'MAXMIN'
inline retType max(const type1 s1, const type2 s2) \
               ^

[tomislav_maric]

The relationship between the GeometricField and the DimensionedField is a clear cut example of why inheriting for implementation does not work.

On one hand side, there is a public inheritance from DimensionedField:

template<class Type, template<class> class PatchField, class GeoMesh>
class GeometricField
:
    public DimensionedField<Type, GeoMesh>
{

and the OO principles state that the GeometricField IS-A DimensionedField, and should behave like that. As a result, following the Liskov substitution principle, every time a GeometricField is passed as a parameter to a function that takes a (const) reference of a DimensionedField, the call would work. It doesn't for assignment: the assignment operators have all been overriden (overloaded for some types, but still they hide the assignment operators of the DimensionedField), in the GeometricField:

    // Member operators

        void operator=(const GeometricField<Type, PatchField, GeoMesh>&);
        void operator=(const tmp<GeometricField<Type, PatchField, GeoMesh> >&);
        void operator=(const dimensioned<Type>&);

        void operator==(const tmp<GeometricField<Type, PatchField, GeoMesh> >&);
        void operator==(const dimensioned<Type>&);

As a result, GeometricField is ... well... not working with DimensionedField. Such inheritance for implementation was maybe done on purpose, since the boundary fields are also dimensioned fields.

So if someone would code:

dimensionedScalarField p_1;
volScalarField p2 = p_1;

the semantical question would be "Should the internal field be assigned to, or one of the boundary fields?". To me, that question is retorical, the assignment operators should use those in DimensionedField, and the internal field should be assigned to. Besides that, 'operator==' was used precisely to introduce assignment to boundaries, so 'operator=' assignes to internal fields. The access to internal field is actually supported by the member function in GeometricField that gives a non-const access to the internal field:

        //- Return dimensioned internal field
        DimensionedInternalField& dimensionedInternalField();

which

1. breaks encapsulation
2. forces the user to write the assignment by force

Also, you have placed the entire function implementation in a return statement. There is a reason functions have bodies. Put some code there. Coding train wrecks like this one:

Cs_*sqrt((kappa_*zeta2_)/((beta(F1)/Cmu_)-gamma(F1)))*pow(mesh_.V(),(1.0/3.0))





will make it very difficult for anyone trying to maintain your code to debug it and extend it.

Also, if you are working with temporary objects, take a look at how they are used in the code, you can see how the interpolation/convection/diffusion schemes work, how they first initialize the temporary, operate on it, then return it. Make the code humanly readable, you will be grateful for it after enough time passes and you forget what you did - and so will all your users.

Here is a working example with that what is available in the GeometricField, off the top of my head. I am not sure if you can copy-paste-use that directly, but it is a start:

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:


tmp<volScalarField> Lvk(
    const volScalarField& S2,
    const volScalarField& F1
) 
{
    const Time& runTime = S2.time(); 
    const fvMesh& mesh = S2.mesh(); 

    // Initialize the temporary object that is to be returned.
    tmp<volScalarField> resultTmp ( 
        new volScalarField (
            IOobject (
                "Lvk",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh, 
            dimensionedScalar (
                "zero", 
                S2.dimensions(),  // Judging by your max() call, you are returning S2 dimensions
                0
            )
        )
    );

    // Get the non-const reference to the temporary object.
    volScalarField& result = resultTmp(); 

    // Keep dimensions of the result but re-name the field. 
    volScalarField h ("h", result);

    // Get the non-const ref to the internal field which is of dimensioned field type.
    DimensionedField<scalar, volMesh>& hInternal = h.dimensionedInternalField(); 

    // Set the internal field values. 
    hInternal = pow(mesh.V(), 1.0 / 3.0); 

    // Calculate the result. 
    // FIXME: separate your huge return call into separate clearly readable calculations and put 
    // them here.  
    
    result = Foam::max(S2, h); 

    return resultTmp;
}

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"


    volScalarField S2 
    (
        IOobject
        (
            "S2",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    volScalarField F1 
    (
        IOobject
        (
            "F1",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

    tmp<volScalarField> result = Lvk(S2, F1); 

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

    Info<< "\nEnd\n" << endl;
    return 0;
}

Edit: also, try extracting your problem into a minimal working example, like the application I posted above, you'll understand your own problem better and get answers sooner.

[zordiack]

Wow, thank you so much for the detailed explanation and example. As my coding skill level is currently mostly "copy-paste", this kind of information is really helpful and I might actually learn something by doing this.

I will try to implement this as soon as possible and I'll split my calculation into smaller pieces, I promise

[tomislav_maric]

@zordiack: No problem, I'm glad it helped!

[zordiack]

The solver compiled and it runs! Thanks again for the advice Now I have to test that it gives reasonable and physically sound results.