[Chikiton]

Hi foamers,
I'm doing a multhiphase task with 3 phases using multiphaseInterFoam and i want to get volume integral of the some variables and of their multiplication, for example, mass of one phase.
in the controlDict file, I want to use a function to calculate the integral for the quantity ( rho*alpha.water) during the whole calculation.
I has found this topic - integrate several fields , but in last version OpenFOAM some functions of swak4foam was included in the OpenFOAM.
I want use smth like this (code from link above)

functions
{
    makeMagU {
       type expressionField;
         functionObjectLibs
        (
            "libswakFunctionObjects.so"
        );      autowrite false;
       fieldName magU;
       expression "mag(U)";
    }
    temp_average
    {
        type volumeAverage;
        functionObjectLibs
        (
            "libsimpleFunctionObjects.so"
        );
        verbose true;
        fields (magU);
    }
}

but I don't know how to implement this in version 1912. Maybe anybody knows. Please help. Thanks.

[crubio.abujas]

Possibly there are other options, but you can try with coded functionObjects. Just define the following on your controlDict:

waterDensity
{
    functionObjectLibs ("libutilityFunctionObjects.so");
    type         coded;
    name       waterDensity;
    codeWrite
    #{
        volScalarField rhoWater
        (
            IOobject
            (   
                "rho.water",
                mesh().time().timeName(),
                mesh(),
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh().lookupObject<volScalarField>("alpha.water") * 
            mesh().lookupObject<volScalarField>("rho")
        );
        rhoWater.write(); 
    #};
}

I've tried it in OF7 and worked fine to me. I guess that you can make more complex field combinations with minor modifications on the code.

[HPE]

Hi,

Just a quick heads up: crubio.abujas' good example will create a new field every execution step, which may hamper the speed of the simulation, or would not allow you to use this function object output for the subsequent function objects, for example fieldAverage of the multiplied field.

I had had provided an example of coded function objects in which a field was created once in one of my threads - for someone from Loughborough Uni. - I have searched for it, but couldn't find in a short period of time. If fancy, please do search for that thread.

Good luck.

[Chikiton]

Originally Posted by crubio.abujas

Possibly there are other options, but you can try with coded functionObjects. Just define the following on your controlDict:

waterDensity
{
    functionObjectLibs ("libutilityFunctionObjects.so");
    type         coded;
    name       waterDensity;
    codeWrite
    #{
        volScalarField rhoWater
        (
            IOobject
            (   
                "rho.water",
                mesh().time().timeName(),
                mesh(),
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            mesh().lookupObject<volScalarField>("alpha.water") * 
            mesh().lookupObject<volScalarField>("rho")
        );
        rhoWater.write(); 
    #};
}

I've tried it in OF7 and worked fine to me. I guess that you can make more complex field combinations with minor modifications on the code.

Thank you crubio.abujas, but as HPE specified, it would not allow me to use this function object output for volIntegrate, code

volumeIntegrate
        type            volFieldValue;
        libs            ("libfieldFunctionObjects.so");

        log             true;
        writeControl    timeStep;
        writeFields     false;

        regionType      all;.
        operation       volIntegrate;

        fields
        (
         rho.water
............
        );
    }.

and then in folder "postProcessing" file "0" does not contain the necessary data. crubio.abujas, do you know solution to this problem? Thank you again.

HPE, thank you, I will definitely look for this topic because I have run into the problem you described.

[Chikiton]

Originally Posted by HPE

Hi,

Just a quick heads up: crubio.abujas' good example will create a new field every execution step, which may hamper the speed of the simulation, or would not allow you to use this function object output for the subsequent function objects, for example fieldAverage of the multiplied field.

I had had provided an example of coded function objects in which a field was created once in one of my threads - for someone from Loughborough Uni. - I have searched for it, but couldn't find in a short period of time. If fancy, please do search for that thread.

Good luck.

HPE, are you talking about this post?Uprime2Mean and R (LES)

Dou you think that using this code I can use volIntegrate for my field?
Thank you!

[crubio.abujas]

Originally Posted by HPE

Hi,

Just a quick heads up: crubio.abujas' good example will create a new field every execution step, which may hamper the speed of the simulation, or would not allow you to use this function object output for the subsequent function objects, for example fieldAverage of the multiplied field.

Good luck.

Concerning the first issue, it is true that the code as stated in the previous post will do it on every execution step. But I think that can be controller adding writeControl/writeInverval properties inside the definition of the function, so it behaves as any other functionObject.


Concerning the second issue: I think that,as the declaration of the rhoAir object is made on the fly, other functionObjects are not aware of their existence. I've to match the writeIntervals to check if the second functionObject will try to read the file, but I guess that the dbregister drops the field as soon as it leave the scope of the first one, so it cannot be located for other ones.

[crubio.abujas]

I've been trying to figure out a way, but right now I have nothing... There has to be a better way for sure, but maybe you can integrate the field using paraview, or you can try to calculate the volIntegral inside the coded functionObject. Not specially satisfying options, but It's all I have in mind right now.

[Chikiton]

Originally Posted by crubio.abujas

I've been trying to figure out a way, but right now I have nothing... There has to be a better way for sure, but maybe you can integrate the field using paraview, or you can try to calculate the volIntegral inside the coded functionObject. Not specially satisfying options, but It's all I have in mind right now.

Thanks anyway! I will look for solutions further.

[crubio.abujas]

Originally Posted by crubio.abujas

I've been trying to figure out a way, but right now I have nothing... There has to be a better way for sure, but maybe you can integrate the field using paraview, or you can try to calculate the volIntegral inside the coded functionObject. Not specially satisfying options, but It's all I have in mind right now.

Sorry for the quick succesions of replies. I was taking an eye onto the HPE post you linked and I've just seen that. If the problem is that when it leave the functionObject scope it drops the IOobject referente, just set it as static as HPE had done in his post. I've tried in a dummy-test and worked fine. Let's hope it works for you as well!

waterDensity
{
    functionObjectLibs ("libutilityFunctionObjects.so");
    type            coded;
    name            waterDensity;
    writeControl    adjustableRunTime;
    writeInterval   0.125;
    codeWrite
    #{

        Info << "Writing rhoWater field" << endl;

        if (!mesh().foundObject<volScalarField>("rho.water")) 
        {
            static autoPtr<volScalarField> rhoWater;
            rhoWater.set
            (
                new volScalarField
                (
                    IOobject
                    (   
                        "rho.water",
                        mesh().time().timeName(),
                        mesh(),
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh(), 
                    mesh().lookupObject<volScalarField>("rho").dimensions()
                )
            );
        }

        if(mesh().time().timeIndex() != 1)
        {
           volScalarField& rhoWater = 
                mesh().lookupObjectRef<volScalarField>("rho.water");
            rhoWater = 
                mesh().lookupObject<volScalarField>("alpha.water") * 
                mesh().lookupObject<volScalarField>("rho");

            if (mesh().time().outputTime())
            {
                rhoWater.write();
            }
        }
    #};

}

EDIT: The code I've posted ran with errors, but never updates the field. I've studied with more calm the HPE code and adapted for the current case. Now it just create a single volScalarField in the register the first time is called, and update it of the subsequent calls.

[Chikiton]

crubio.abujas, thank you so much!!! it works.
Another little question, I need integral : rho*Uz * dV, where Uz - vertical velocity component. How can I access this vertical component of velocity in your code?

[crubio.abujas]

Originally Posted by Chikiton

crubio.abujas, thank you so much!!! it works.
Another little question, I need integral : rho*Uz * dV, where Uz - vertical velocity component. How can I access this vertical component of velocity in your code?

Hi again Chikiton,
For accessing the component of a vector you can try this code:

zMomentum
{
    functionObjectLibs ("libutilityFunctionObjects.so");
    type            coded;
    name            zMomentum;
    writeControl    adjustableRunTime;
    writeInterval   0.125;
    codeWrite
    #{

        if (!mesh().foundObject<volScalarField>("zMom")) 
        {
            static autoPtr<volScalarField> zMom;
            zMom.set
            (
                new volScalarField
                (
                    IOobject
                    (   
                        "zMom",
                        mesh().time().timeName(),
                        mesh(),
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh(), 
                    dimMass/dimVolume*dimLength/dimTime
                )
            );
        }

       volScalarField& zMom = 
            mesh().lookupObjectRef<volScalarField>("zMom");
        zMom = 
            mesh().lookupObject<volScalarField>("rho") * 
            mesh().lookupObject<volVectorField>("U").component(2);

        zMom.write();
    #};
}

Just ensure that the new objects has the proper units. When refering the field in the assignment (zMom = ...) ensure as well to use the proper type, in this case Vector. As in the end you want a scalar magnitude you have to use just one of the components, I think component method starts with 0(x), and go up until 2(z). But check it.


Referring the volume multiplication. I tried to multiply straightforward, but it fails. The reason is that volume is not defined in the boundaries, so to have a magnitude of the mass could be a little bit more complicated inside the coded function objects. Maybe you can do that afterwards with another functionObject?


Hopes it works!

[Chikiton]

Originally Posted by crubio.abujas

Hi again Chikiton,
For accessing the component of a vector you can try this code:

zMomentum
{
    functionObjectLibs ("libutilityFunctionObjects.so");
    type            coded;
    name            zMomentum;
    writeControl    adjustableRunTime;
    writeInterval   0.125;
    codeWrite
    #{

        if (!mesh().foundObject<volScalarField>("zMom")) 
        {
            static autoPtr<volScalarField> zMom;
            zMom.set
            (
                new volScalarField
                (
                    IOobject
                    (   
                        "zMom",
                        mesh().time().timeName(),
                        mesh(),
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh(), 
                    dimMass/dimVolume*dimLength/dimTime
                )
            );
        }

       volScalarField& zMom = 
            mesh().lookupObjectRef<volScalarField>("zMom");
        zMom = 
            mesh().lookupObject<volScalarField>("rho") * 
            mesh().lookupObject<volVectorField>("U").component(2);

        zMom.write();
    #};
}

Just ensure that the new objects has the proper units. When refering the field in the assignment (zMom = ...) ensure as well to use the proper type, in this case Vector. As in the end you want a scalar magnitude you have to use just one of the components, I think component method starts with 0(x), and go up until 2(z). But check it.


Referring the volume multiplication. I tried to multiply straightforward, but it fails. The reason is that volume is not defined in the boundaries, so to have a magnitude of the mass could be a little bit more complicated inside the coded function objects. Maybe you can do that afterwards with another functionObject?


Hopes it works!

Thank you! I can do the integral with volIntegrate for new field zMom. Thank you so much again!!!

[olesen]

Originally Posted by crubio.abujas

Hi again Chikiton,
For accessing the component of a vector you can try this code:

zMomentum
{
    functionObjectLibs ("libutilityFunctionObjects.so");
    type            coded;
    name            zMomentum;
    writeControl    adjustableRunTime;
    writeInterval   0.125;
    codeWrite
    #{

        if (!mesh().foundObject<volScalarField>("zMom")) 
        {
            static autoPtr<volScalarField> zMom;
            zMom.set
            (
                new volScalarField
                (
                    IOobject
                    (   
                        "zMom",
                        mesh().time().timeName(),
                        mesh(),
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh(), 
                    dimMass/dimVolume*dimLength/dimTime
                )
            );
        }

       volScalarField& zMom = 
            mesh().lookupObjectRef<volScalarField>("zMom");
        zMom = 
            mesh().lookupObject<volScalarField>("rho") * 
            mesh().lookupObject<volVectorField>("U").component(2);

        zMom.write();
    #};
}

Just a note for anyone using this solution. Instead of using the static autoPtr to hold the pointer, it probably makes more sense to store it on the objectRegistry, since it is registered there anyhow.

        if (!mesh().foundObject<volScalarField>("zMom")) 
         {
            volScalarField*  pzmon = new volScalarField
                (
                    IOobject
                    (   
                        "zMom",
                        mesh().time().timeName(),
                        mesh(),
                        IOobject::NO_READ,
                        IOobject::AUTO_WRITE
                    ),
                    mesh(), 
                    dimMass/dimVolume*dimLength/dimTime
                );

            pzmon->store();
         }

       volScalarField& zMom = 
            mesh().lookupObjectRef<volScalarField>("zMom");
        zMom = 
            mesh().lookupObject<volScalarField>("rho") * 
            mesh().lookupObject<volVectorField>("U").component(2);

        zMom.write();
    #};
}

With a more recent version like 1906, 1912 etc I would generally write this more compactly and avoid the double lookups. Can also take advantage of auto deductions:

auto* pzmon = mesh().getObjectPtr("zMom");

if (!pzMon)
{
    pzmon = new volScalarField
    (
        IOobject
        (
            "zMom",
            mesh().time().timeName(),
            mesh(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh(), 
         dimMass/dimVolume*dimLength/dimTime
    );

      pzmon->store();
}

auto& zMom = *pzmon;

zMom = 
     mesh().lookupObject<volScalarField>("rho") * 
     mesh().lookupObject<volVectorField>("U").component(2);

zMom.write();

Not sure if you'd want AUTO_WRITE since you are taking over writing yourself.

[olesen]

Originally Posted by olesen

...
With a more recent version like 1906, 1912 etc I would generally write this more compactly and avoid the double lookups. Can also take advantage of auto deductions:
...

Maybe a bit funny adding more answers to my own post, but figured it would be worthwhile providing a link to the following code snippet:

https://develop.openfoam.com/Develop...erivedFields.C

It illustrates much the same approach as previously outlined, but also includes handling of incompressible/compressible without much additional coding and remains IMO fairly readable in the process.
For clarification, the getObjectPtr method is very much like the findObject method. They both return a pointer that can be tested (as a bool) or de-referenced. The getObjectPtr version just includes an additional const_cast since this is a frequently used idiom and it makes sense not to have all these operations cluttering things up.

[hhu_lulu]

Hi, crubio.abujas, thanks for your contribution!
I tried to adapt your code to my case to calculate the time-averaged flux：<Us>, in which the U and s denotes velocity and scalar field repectively.
The simulation successfully generated the field Us. However, this variable could not be used by fieldAverage. That is no say, there is no UsMean output.
Here's the code adapted from your code(the flux_1 is the variable in place of Us):

fluctuation_flux1
{
libs ("libutilityFunctionObjects.so");
type coded;
name fluctuation_flux1;
timeStart 0.2;
timeEnd 2;
writeControl adjustableRunTime;
writeInterval 2;
codeWrite
#{

Info << "Writing flux1 field" << endl;

if (!mesh().foundObject<volVectorField>("flux_1"))
{
static autoPtr<volVectorField> flux1;
flux1.set
(
new volVectorField
(
IOobject
(
"flux_1",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh(),
mesh().lookupObject<volVectorField>("U").dimension s()
)
);
}

if(mesh().time().timeIndex() != 1)
{
volVectorField& flux1 =
mesh().lookupObjectRef<volVectorField>("flux_1");
flux1 =
mesh().lookupObject<volVectorField>("U") *
mesh().lookupObject<volScalarField>("s");

if (mesh().time().outputTime())
{
flux1.write();
}
}
#};
}

Also the fieldAverage function:

fieldAverage1
{
type fieldAverage;
libs ( "libfieldFunctionObjects.so" );
resetOnOutput false;
startTime 0.5;
endTime 2;
writeControl adjustableRunTime;
writeInterval 2;

fields
(
flux_1
{
mean on;
prime2Mean on;
base time;
}
);
}

And the log from HPC platfoam shows:

request for volVectorField flux_1Mean from objectRegistry region0 failed

[crubio.abujas]

Hi Shengjie Lu,

I’ve made some tests in OF7 and fieldAverage works fine with custom fields. The code should try to reconstruct the fields defined on each time folder, so if flux exist it should process it properly.

In my trials I have been unable to replicate the error you mentioned. Could you detail what commands did you execute and the detailed output? Did you execute “postprocess” or “pimpleFoam -postProcess”? Which version of OpenFOAM did you use?

[hhu_lulu]

Hi crubio.abujas,
I use OF7.0 in clusters, here's the commands in script:

source /data0/software/intel/2017u5/bin/compilervars.sh intel64
source /data0/software/openfoam/7.0/OpenFOAM-7/etc/bashrc
decomposePar
mpirun -np 96 pimpleFoam -parallel
scontrol show job $SLURM_JOBID
reconstructPar

BTW, I only need to output the result in the last timestep. In your original code, the rho.water and rhoWater were specified,but in another code ,you used the same zMom. What's could be the difference? Should I use rho.water as a field name for fieldAverage?

Could you share the controlDict file of your test case if possible.

Thanks again for your quick reply!