CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

How to include cell volumes in calculation of volScalarField?

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 2 Post By tomislav_maric

Reply
 
LinkBack Thread Tools Display Modes
Old   April 9, 2014, 08:55
Default How to include cell volumes in calculation of volScalarField?
  #1
Member
 
Pekka Pasanen
Join Date: Feb 2012
Location: Finland
Posts: 38
Rep Power: 5
zordiack is on a distinguished road
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:

Code:
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:

Code:
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:

Code:
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) \
               ^
zordiack is offline   Reply With Quote

Old   April 9, 2014, 14:28
Default
  #2
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 260
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
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:
Code:
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:

Code:
    // 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:

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:
Code:
        //- 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:

Code:
#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 and aylalisa like this.
tomislav_maric is offline   Reply With Quote

Old   April 10, 2014, 00:36
Default
  #3
Member
 
Pekka Pasanen
Join Date: Feb 2012
Location: Finland
Posts: 38
Rep Power: 5
zordiack is on a distinguished road
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
zordiack is offline   Reply With Quote

Old   April 10, 2014, 12:11
Default
  #4
Senior Member
 
Tomislav Maric
Join Date: Mar 2009
Location: Darmstadt, Germany
Posts: 260
Blog Entries: 5
Rep Power: 11
tomislav_maric is on a distinguished road
@zordiack: No problem, I'm glad it helped!
__________________
Blog: sourceflux.de/blog

"The OpenFOAM Technology Primer": sourceflux.de/book

Twitter: @sourceflux_de

When asking a question, prepare a SSCCE.
tomislav_maric is offline   Reply With Quote

Old   April 11, 2014, 05:57
Default
  #5
Member
 
Pekka Pasanen
Join Date: Feb 2012
Location: Finland
Posts: 38
Rep Power: 5
zordiack is on a distinguished road
The solver compiled and it runs! Thanks again for the advice Now I have to test that it gives reasonable and physically sound results.
zordiack is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Questions on dynamicTopoFvMesh danvica OpenFOAM Running, Solving & CFD 70 June 12, 2015 04:33
interFoam running blowing up sandy13 OpenFOAM Running, Solving & CFD 2 May 5, 2015 07:16
volScalarField for cell volumes and face surfaces AlmostSurelyRob OpenFOAM 2 December 13, 2010 06:24
OpenFOAM15 paraFoam bug koen OpenFOAM Bugs 19 June 30, 2009 10:46
Error #285UNSUPPORTED CELL TYPE IN CALCULATION OF Nahidh Hamid Sharif CD-adapco 0 September 19, 2007 02:00


All times are GMT -4. The time now is 07:11.