Merge branch 'master' of github.com-OpenFOAM:OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/multiphaseEulerFoam/functionObjects/Make/files

@@ -1,3 +1,4 @@
+moments/moments.C
 sizeDistribution/sizeDistribution.C
 phaseForces/phaseForces.C
 phaseMap/phaseMap.C

applications/solvers/multiphase/multiphaseEulerFoam/functionObjects/moments/moments.C

@@ -0,0 +1,972 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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 "moments.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(moments, 0);
+    addToRunTimeSelectionTable(functionObject, moments, dictionary);
+}
+}
+
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum
+    <
+        Foam::functionObjects::moments::momentType,
+        4
+    >::names[] = {"integerMoment", "mean", "variance", "stdDev"};
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::moments::momentType,
+    4
+>
+Foam::functionObjects::moments::momentTypeNames_;
+
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum
+    <
+        Foam::functionObjects::moments::coordinateType,
+        3
+    >::names[] = {"volume", "area", "diameter"};
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::moments::coordinateType,
+    3
+>
+Foam::functionObjects::moments::coordinateTypeNames_;
+
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum
+    <
+        Foam::functionObjects::moments::weightType,
+        3
+    >::names[] =
+    {
+        "numberConcentration",
+        "volumeConcentration",
+        "areaConcentration"
+    };
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::moments::weightType,
+    3
+>
+Foam::functionObjects::moments::weightTypeNames_;
+
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum
+    <
+        Foam::functionObjects::moments::meanType,
+        3
+    >::names[] = {"arithmetic", "geometric", "notApplicable"};
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::moments::meanType,
+    3
+>
+Foam::functionObjects::moments::meanTypeNames_;
+
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+Foam::word Foam::functionObjects::moments::coordinateTypeSymbolicName()
+{
+    word coordinateTypeSymbolicName(word::null);
+
+    switch (coordinateType_)
+    {
+        case coordinateType::volume:
+        {
+            coordinateTypeSymbolicName = "v";
+
+            break;
+        }
+        case coordinateType::area:
+        {
+            coordinateTypeSymbolicName = "a";
+
+            break;
+        }
+        case coordinateType::diameter:
+        {
+            coordinateTypeSymbolicName = "d";
+
+            break;
+        }
+    }
+
+    return coordinateTypeSymbolicName;
+}
+
+
+Foam::word Foam::functionObjects::moments::weightTypeSymbolicName()
+{
+    word weightTypeSymbolicName(word::null);
+
+    switch (weightType_)
+    {
+        case weightType::numberConcentration:
+        {
+            weightTypeSymbolicName = "N";
+
+            break;
+        }
+        case weightType::volumeConcentration:
+        {
+            weightTypeSymbolicName = "V";
+
+            break;
+        }
+        case weightType::areaConcentration:
+        {
+            weightTypeSymbolicName = "A";
+
+            break;
+        }
+    }
+
+    return weightTypeSymbolicName;
+}
+
+
+Foam::word Foam::functionObjects::moments::defaultFldName()
+{
+    word meanName
+    (
+        meanType_ == meanType::geometric
+      ? word(meanTypeNames_[meanType_]).capitalise()
+      : word("")
+    );
+
+    return
+        word
+        (
+            IOobject::groupName
+            (
+                "weighted"
+              + meanName
+              + word(momentTypeNames_[momentType_]).capitalise()
+              + "("
+              + weightTypeSymbolicName()
+              + ","
+              + coordinateTypeSymbolicName()
+              + ")",
+                popBal_.name()
+            )
+        );
+}
+
+
+Foam::word Foam::functionObjects::moments::integerMomentFldName()
+{
+    return
+        word
+        (
+            IOobject::groupName
+            (
+                word(momentTypeNames_[momentType_])
+              + Foam::name(order_)
+              + "("
+              + weightTypeSymbolicName()
+              + ","
+              + coordinateTypeSymbolicName()
+              + ")",
+                popBal_.name()
+            )
+        );
+}
+
+
+void Foam::functionObjects::moments::setDimensions
+(
+    volScalarField& fld,
+    momentType momType
+)
+{
+    switch (momType)
+    {
+        case momentType::integerMoment:
+        {
+            switch (coordinateType_)
+            {
+                case coordinateType::volume:
+                {
+                    fld.dimensions().reset
+                    (
+                        pow(dimVolume, order_)/dimVolume
+                    );
+
+                    break;
+                }
+                case coordinateType::area:
+                {
+                    fld.dimensions().reset
+                    (
+                        pow(dimArea, order_)/dimVolume
+                    );
+
+                    break;
+                }
+                case coordinateType::diameter:
+                {
+                    fld.dimensions().reset
+                    (
+                        pow(dimLength, order_)/dimVolume
+                    );
+
+                    break;
+                }
+            }
+
+            switch (weightType_)
+            {
+                case weightType::volumeConcentration:
+                {
+                    fld.dimensions().reset(fld.dimensions()*dimVolume);
+
+                    break;
+                }
+                case weightType::areaConcentration:
+                {
+                    fld.dimensions().reset(fld.dimensions()*dimArea);
+
+                    break;
+                }
+                default:
+                {
+                    break;
+                }
+            }
+
+            break;
+        }
+        case momentType::mean:
+        {
+            switch (coordinateType_)
+            {
+                case coordinateType::volume:
+                {
+                    fld.dimensions().reset(dimVolume);
+
+                    break;
+                }
+                case coordinateType::area:
+                {
+                    fld.dimensions().reset(dimArea);
+
+                    break;
+                }
+                case coordinateType::diameter:
+                {
+                    fld.dimensions().reset(dimLength);
+
+                    break;
+                }
+            }
+
+            break;
+        }
+        case momentType::variance:
+        {
+            switch (coordinateType_)
+            {
+                case coordinateType::volume:
+                {
+                    fld.dimensions().reset(sqr(dimVolume));
+
+                    break;
+                }
+                case coordinateType::area:
+                {
+                    fld.dimensions().reset(sqr(dimArea));
+
+                    break;
+                }
+                case coordinateType::diameter:
+                {
+                    fld.dimensions().reset(sqr(dimLength));
+
+                    break;
+                }
+            }
+
+            if (meanType_ == meanType::geometric)
+            {
+                fld.dimensions().reset(dimless);
+            }
+
+            break;
+        }
+        case momentType::stdDev:
+        {
+            switch (coordinateType_)
+            {
+                case coordinateType::volume:
+                {
+                    fld.dimensions().reset(dimVolume);
+
+                    break;
+                }
+                case coordinateType::area:
+                {
+                    fld.dimensions().reset(dimArea);
+
+                    break;
+                }
+                case coordinateType::diameter:
+                {
+                    fld.dimensions().reset(dimLength);
+
+                    break;
+                }
+            }
+
+            if (meanType_ == meanType::geometric)
+            {
+                fld.dimensions().reset(dimless);
+            }
+
+            break;
+        }
+    }
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::moments::totalConcentration()
+{
+    tmp<volScalarField> tTotalConcentration
+    (
+        volScalarField::New
+        (
+            "totalConcentration",
+            mesh_,
+            dimensionedScalar(inv(dimVolume), Zero)
+        )
+    );
+
+    volScalarField& totalConcentration = tTotalConcentration.ref();
+
+    switch (weightType_)
+    {
+        case weightType::volumeConcentration:
+        {
+            totalConcentration.dimensions().reset
+            (
+                totalConcentration.dimensions()*dimVolume
+            );
+
+            break;
+        }
+        case weightType::areaConcentration:
+        {
+            totalConcentration.dimensions().reset
+            (
+                totalConcentration.dimensions()*dimArea
+            );
+
+            break;
+        }
+        default:
+        {
+            break;
+        }
+    }
+
+    forAll(popBal_.sizeGroups(), i)
+    {
+        const Foam::diameterModels::sizeGroup& fi = popBal_.sizeGroups()[i];
+
+        switch (weightType_)
+        {
+            case weightType::numberConcentration:
+            {
+                totalConcentration += fi*fi.phase()/fi.x();
+
+                break;
+            }
+            case weightType::volumeConcentration:
+            {
+                totalConcentration += fi*fi.phase();
+
+                break;
+            }
+            case weightType::areaConcentration:
+            {
+                totalConcentration += fi.a()*fi*fi.phase()/fi.x();
+
+                break;
+            }
+        }
+    }
+
+    return tTotalConcentration;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::functionObjects::moments::mean()
+{
+    tmp<volScalarField> tMean
+    (
+        volScalarField::New
+        (
+            "mean",
+            mesh_,
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+
+    volScalarField& mean = tMean.ref();
+
+    setDimensions(mean, momentType::mean);
+
+    volScalarField totalConcentration(this->totalConcentration());
+
+    forAll(popBal_.sizeGroups(), i)
+    {
+        const Foam::diameterModels::sizeGroup& fi = popBal_.sizeGroups()[i];
+
+        volScalarField concentration(fi*fi.phase()/fi.x());
+
+        switch (weightType_)
+        {
+            case weightType::volumeConcentration:
+            {
+                concentration *= fi.x();
+
+                break;
+            }
+            case weightType::areaConcentration:
+            {
+                concentration *= fi.a();
+
+                break;
+            }
+            default:
+            {
+                break;
+            }
+        }
+
+        switch (meanType_)
+        {
+            case meanType::geometric:
+            {
+                mean.dimensions().reset(dimless);
+
+                switch (coordinateType_)
+                {
+                    case coordinateType::volume:
+                    {
+                        dimensionedScalar unitVolume(dimVolume, 1);
+
+                        mean +=
+                            Foam::log(fi.x()/unitVolume)
+                           *concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::area:
+                    {
+                        dimensionedScalar unitArea(dimArea, 1);
+
+                        mean +=
+                            Foam::log(fi.a()/unitArea)
+                           *concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::diameter:
+                    {
+                        dimensionedScalar unitLength(dimLength, 1);
+
+                        mean +=
+                            Foam::log(fi.d()/unitLength)
+                           *concentration/totalConcentration;
+
+                        break;
+                    }
+                }
+
+                break;
+            }
+            default:
+            {
+                switch (coordinateType_)
+                {
+                    case coordinateType::volume:
+                    {
+                        mean += fi.x()*concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::area:
+                    {
+                        mean += fi.a()*concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::diameter:
+                    {
+                        mean += fi.d()*concentration/totalConcentration;
+
+                        break;
+                    }
+                }
+
+                break;
+            }
+        }
+    }
+
+    if (meanType_ == meanType::geometric)
+    {
+        mean = exp(mean);
+
+        setDimensions(mean, momentType::mean);
+    }
+
+    return tMean;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::functionObjects::moments::variance()
+{
+    tmp<volScalarField> tVariance
+    (
+        volScalarField::New
+        (
+            "variance",
+            mesh_,
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+
+    volScalarField& variance = tVariance.ref();
+
+    setDimensions(variance, momentType::variance);
+
+    volScalarField totalConcentration(this->totalConcentration());
+    volScalarField mean(this->mean());
+
+    forAll(popBal_.sizeGroups(), i)
+    {
+        const Foam::diameterModels::sizeGroup& fi = popBal_.sizeGroups()[i];
+
+        volScalarField concentration(fi*fi.phase()/fi.x());
+
+        switch (weightType_)
+        {
+            case weightType::volumeConcentration:
+            {
+                concentration *= fi.x();
+
+                break;
+            }
+            case weightType::areaConcentration:
+            {
+                concentration *= fi.a();
+
+                break;
+            }
+            default:
+            {
+                break;
+            }
+        }
+
+        switch (meanType_)
+        {
+            case meanType::geometric:
+            {
+                switch (coordinateType_)
+                {
+                    case coordinateType::volume:
+                    {
+                        variance +=
+                            sqr(Foam::log(fi.x()/mean))
+                           *concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::area:
+                    {
+                        variance +=
+                            sqr(Foam::log(fi.a()/mean))
+                           *concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::diameter:
+                    {
+                        variance +=
+                            sqr(Foam::log(fi.d()/mean))
+                           *concentration/totalConcentration;
+
+                        break;
+                    }
+                }
+
+                break;
+            }
+            default:
+            {
+                switch (coordinateType_)
+                {
+                    case coordinateType::volume:
+                    {
+                        variance +=
+                            sqr(fi.x() - mean)*concentration/totalConcentration;
+
+
+                        break;
+                    }
+                    case coordinateType::area:
+                    {
+                        variance +=
+                            sqr(fi.a() - mean)*concentration/totalConcentration;
+
+                        break;
+                    }
+                    case coordinateType::diameter:
+                    {
+                        variance +=
+                            sqr(fi.d() - mean)*concentration/totalConcentration;
+
+                        break;
+                    }
+                }
+
+                break;
+            }
+        }
+    }
+
+    return tVariance;
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::moments::stdDev()
+{
+    switch (meanType_)
+    {
+        case meanType::geometric:
+        {
+            return exp(sqrt(this->variance()));
+        }
+        default:
+        {
+            return sqrt(this->variance());
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::moments::moments
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    popBal_
+    (
+        obr_.lookupObject<Foam::diameterModels::populationBalanceModel>
+        (
+            dict.lookup("populationBalance")
+        )
+    ),
+    momentType_(momentTypeNames_.read(dict.lookup("momentType"))),
+    coordinateType_(coordinateTypeNames_.read(dict.lookup("coordinateType"))),
+    weightType_
+    (
+        dict.found("weight")
+      ? weightTypeNames_.read(dict.lookup("weightType"))
+      : weightType::numberConcentration
+    ),
+    meanType_(meanType::notApplicable),
+    order_(-1),
+    fldPtr_(nullptr)
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::functionObjects::moments::~moments()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::moments::read(const dictionary& dict)
+{
+    fvMeshFunctionObject::read(dict);
+
+    switch (momentType_)
+    {
+        case momentType::integerMoment:
+        {
+            order_ = dict.lookup<int>("order");
+
+            break;
+        }
+        default:
+        {
+            meanType_ =
+                dict.found("meanType")
+              ? meanTypeNames_.read(dict.lookup("meanType"))
+              : meanType::arithmetic;
+
+            break;
+        }
+    }
+
+    switch (momentType_)
+    {
+        case momentType::integerMoment:
+        {
+            fldPtr_.set
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        this->integerMomentFldName(),
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    mesh_,
+                    dimensionedScalar(dimless, Zero)
+                )
+            );
+
+            volScalarField& integerMoment = fldPtr_();
+
+            setDimensions(integerMoment, momentType::integerMoment);
+
+            break;
+        }
+        case momentType::mean:
+        {
+            fldPtr_.set
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        this->defaultFldName(),
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    this->mean()
+                )
+            );
+
+            break;
+        }
+        case momentType::variance:
+        {
+            fldPtr_.set
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        this->defaultFldName(),
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    this->variance()
+                )
+            );
+
+            break;
+        }
+        case momentType::stdDev:
+        {
+            fldPtr_.set
+            (
+                new volScalarField
+                (
+                    IOobject
+                    (
+                        this->defaultFldName(),
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    this->stdDev()
+                )
+            );
+
+            break;
+        }
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::moments::execute()
+{
+    switch (momentType_)
+    {
+        case momentType::integerMoment:
+        {
+            volScalarField& integerMoment = fldPtr_();
+
+            integerMoment = Zero;
+
+            forAll(popBal_.sizeGroups(), i)
+            {
+                const Foam::diameterModels::sizeGroup& fi =
+                    popBal_.sizeGroups()[i];
+
+                volScalarField concentration(fi*fi.phase()/fi.x());
+
+                switch (weightType_)
+                {
+                    case weightType::volumeConcentration:
+                    {
+                        concentration *= fi.x();
+
+                        break;
+                    }
+                    case weightType::areaConcentration:
+                    {
+                        concentration *= fi.a();
+
+                        break;
+                    }
+                    default:
+                    {
+                        break;
+                    }
+                }
+
+                switch (coordinateType_)
+                {
+                    case coordinateType::volume:
+                    {
+                        integerMoment +=
+                            pow(fi.x(), order_)*concentration;
+
+                        break;
+                    }
+                    case coordinateType::area:
+                    {
+                        integerMoment +=
+                            pow(fi.a(), order_)*concentration;
+
+                        break;
+                    }
+                    case coordinateType::diameter:
+                    {
+                        integerMoment +=
+                            pow(fi.d(), order_)*concentration;
+
+                        break;
+                    }
+                }
+            }
+
+            break;
+        }
+        case momentType::mean:
+        {
+            fldPtr_() = this->mean();
+
+            break;
+        }
+        case momentType::variance:
+        {
+            fldPtr_() = this->variance();
+
+            break;
+        }
+        case momentType::stdDev:
+        {
+            fldPtr_() = sqrt(this->variance());
+
+            break;
+        }
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::moments::write()
+{
+    writeObject(fldPtr_->name());
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/functionObjects/moments/moments.H

@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2022 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::functionObjects::moments
+
+Description
+    Calculates and writes out integral (integer moments) or mean properties
+    (mean, variance, standard deviation) of a size distribution computed with
+    multiphaseEulerFoam. Requires solver post-processing.
+
+    The following function object specification for example returns the first
+    moment of the volume-based number density function which is equivalent to
+    the phase fraction of the particulate phase:
+
+    \verbatim
+    moments
+    {
+        type                moments;
+        libs                ("libmultiphaseEulerFoamFunctionObjects.so");
+        executeControl      timeStep;
+        writeControl        writeTime;
+        populationBalance   bubbles;
+        momentType          integerMoment;
+        coordinateType      volume;
+        order               1;
+    }
+    \endverbatim
+
+Usage
+    \table
+        Property          | Description              | Required  | Default
+        populationBalance | population balance name  | yes       |
+        momentType        | desired moment of the distribution\\
+                                                     | yes       |
+        coordinateType    | particle property        | yes       |
+        weightType        | number/volume/area concentration\\
+                                                     | no\\
+                                                     | numberConcentration
+        order             | order of integer moment  | for integer moments |
+        meanType          | arithmetic or geometric  | for non-integer moments\\
+                                                                 | arithmetic
+    \endtable
+
+See also
+    Foam::diameterModels::populationBalanceModel
+    Foam::functionObjects::fvMeshFunctionObject
+    Foam::functionObject
+
+SourceFiles
+    moments.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_moments_H
+#define functionObjects_moments_H
+
+#include "fvMeshFunctionObject.H"
+#include "populationBalanceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class moments Declaration
+\*---------------------------------------------------------------------------*/
+
+class moments
+:
+    public fvMeshFunctionObject
+{
+public:
+
+    //- Enumeration for the moment types
+    enum class momentType
+    {
+        integerMoment,
+        mean,
+        variance,
+        stdDev
+    };
+
+    //- Names of the moment types
+    static const NamedEnum<momentType, 4> momentTypeNames_;
+
+    //- Enumeration for the coordinate types
+    enum class coordinateType
+    {
+        volume,
+        area,
+        diameter
+    };
+
+    //- Names of the coordinate types
+    static const NamedEnum<coordinateType, 3> coordinateTypeNames_;
+
+    //- Enumeration for the weight types
+    enum class weightType
+    {
+        numberConcentration,
+        volumeConcentration,
+        areaConcentration
+    };
+
+    //- Names of the weight types
+    static const NamedEnum<weightType, 3> weightTypeNames_;
+
+    //- Enumeration for the mean types
+    enum class meanType
+    {
+        arithmetic,
+        geometric,
+        notApplicable
+    };
+
+    //- Names of the mean types
+    static const NamedEnum<meanType, 3> meanTypeNames_;
+
+
+private:
+
+    // Private Data
+
+        //- Reference to population balance
+        const Foam::diameterModels::populationBalanceModel& popBal_;
+
+        //- Moment type
+        momentType momentType_;
+
+        //- Coordinate type
+        coordinateType coordinateType_;
+
+        //- Weight type
+        weightType weightType_;
+
+        //- Mean type
+        meanType meanType_;
+
+        //- Integer moment order
+        int order_;
+
+        //- Result field
+        autoPtr<volScalarField> fldPtr_;
+
+
+    // Private Member Functions
+
+        //- Coordinate type symbolic name for shorter field names
+        word coordinateTypeSymbolicName();
+
+        //- Weight type symbolic name for shorter field names
+        word weightTypeSymbolicName();
+
+        //- Default field name
+        word defaultFldName();
+
+        //- Integer moment field name
+        word integerMomentFldName();
+
+        //- Set dimensions
+        void setDimensions(volScalarField& fld, momentType momType);
+
+        //- Total concentration
+        tmp<volScalarField> totalConcentration();
+
+        //- Mean value
+        tmp<volScalarField> mean();
+
+        //- Variance
+        tmp<volScalarField> variance();
+
+        //- Standard deviation
+        tmp<volScalarField> stdDev();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("moments");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        moments
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary&
+        );
+
+        //- Disallow default bitwise copy construction
+        moments(const moments&) = delete;
+
+
+    //- Destructor
+    virtual ~moments();
+
+
+    // Member Functions
+
+        //- Read the data
+        virtual bool read(const dictionary&);
+
+        //- Return the list of fields required
+        virtual wordList fields() const
+        {
+            return wordList::null();
+        }
+
+        //- Calculate the force fields
+        virtual bool execute();
+
+        //- Write the force fields
+        virtual bool write();
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const moments&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/functionObjects/sizeDistribution/sizeDistribution.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,10 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "sizeDistribution.H"
-#include "Time.H"
-#include "fvMesh.H"
 #include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -45,19 +42,21 @@ const char*
 Foam::NamedEnum
 <
     Foam::functionObjects::sizeDistribution::functionType,
-    4
+    6
 >::names[] =
 {
-    "moments",
+    "numberConcentration",
-    "standardDeviation",
+    "numberDensity",
-    "number",
+    "volumeConcentration",
-    "volume"
+    "volumeDensity",
+    "areaConcentration",
+    "areaDensity"
 };
 
 const Foam::NamedEnum
 <
     Foam::functionObjects::sizeDistribution::functionType,
-    4
+    6
 > Foam::functionObjects::sizeDistribution::functionTypeNames_;
 
 template<>
@@ -80,10 +79,122 @@ const Foam::NamedEnum
     4
 > Foam::functionObjects::sizeDistribution::coordinateTypeNames_;
 
+
+namespace Foam
+{
+    template<>
+    const char* NamedEnum
+    <
+        Foam::functionObjects::sizeDistribution::weightType,
+        4
+    >::names[] =
+    {
+        "numberConcentration",
+        "volumeConcentration",
+        "areaConcentration",
+        "cellVolume"
+    };
+}
+
+
+const Foam::NamedEnum
+<
+    Foam::functionObjects::sizeDistribution::weightType,
+    4
+>
+Foam::functionObjects::sizeDistribution::weightTypeNames_;
+
 using Foam::constant::mathematical::pi;
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+Foam::word Foam::functionObjects::sizeDistribution::functionTypeSymbolicName()
+{
+    word functionTypeSymbolicName(word::null);
+
+    switch (functionType_)
+    {
+        case functionType::numberConcentration:
+        {
+            functionTypeSymbolicName = "N";
+
+            break;
+        }
+        case functionType::numberDensity:
+        {
+            functionTypeSymbolicName = "n";
+
+            break;
+        }
+        case functionType::volumeConcentration:
+        {
+            functionTypeSymbolicName = "V";
+
+            break;
+        }
+        case functionType::volumeDensity:
+        {
+            functionTypeSymbolicName = "v";
+
+            break;
+        }
+        case functionType::areaConcentration:
+        {
+            functionTypeSymbolicName = "A";
+
+            break;
+        }
+        case functionType::areaDensity:
+        {
+            functionTypeSymbolicName = "a";
+
+            break;
+        }
+    }
+
+    return functionTypeSymbolicName;
+}
+
+
+Foam::word Foam::functionObjects::sizeDistribution::coordinateTypeSymbolicName
+(
+    const coordinateType& cType
+)
+{
+    word coordinateTypeSymbolicName(word::null);
+
+    switch (cType)
+    {
+        case coordinateType::volume:
+        {
+            coordinateTypeSymbolicName = "v";
+
+            break;
+        }
+        case coordinateType::area:
+        {
+            coordinateTypeSymbolicName = "a";
+
+            break;
+        }
+        case coordinateType::diameter:
+        {
+            coordinateTypeSymbolicName = "d";
+
+            break;
+        }
+        case coordinateType::projectedAreaDiameter:
+        {
+            coordinateTypeSymbolicName = "dPa";
+
+            break;
+        }
+    }
+
+    return coordinateTypeSymbolicName;
+}
+
 
 Foam::tmp<Foam::scalarField>
 Foam::functionObjects::sizeDistribution::filterField
@@ -102,252 +213,120 @@ Foam::functionObjects::sizeDistribution::filterField
 }
 
 
-void Foam::functionObjects::sizeDistribution::correctVolAverages()
+Foam::scalar Foam::functionObjects::sizeDistribution::averageCoordinateValue
-{
-    forAll(N_, i)
-    {
-        const Foam::diameterModels::sizeGroup& fi = popBal_.sizeGroups()[i];
-
-        scalarField Ni(filterField(fi*fi.phase()/fi.x()));
-        scalarField V(filterField(mesh_.V()));
-        scalarField ai(filterField(fi.a()));
-        scalarField di(Ni.size());
-
-        switch (coordinateType_)
-        {
-            case ctProjectedAreaDiameter:
-            {
-                di = sqrt(ai/pi);
-
-                break;
-            }
-            default:
-            {
-                di = fi.d();
-
-                break;
-            }
-        }
-
-        N_[i] = gSum(V*Ni)/this->V();
-        V_[i] = fi.x().value();
-        a_[i] = gSum(V*ai)/this->V();
-        d_[i] = gSum(V*di)/this->V();
-    }
-}
-
-
-void Foam::functionObjects::sizeDistribution::writeMoments()
-{
-    logFiles::write();
-
-    Log << "    writing moments of size distribution." << endl;
-
-    if (Pstream::master())
-    {
-        writeTime(file());
-    }
-
-    const scalarField& bin = this->bin();
-
-    for (label k = 0; k <= maxOrder_; k++)
-    {
-        scalar result = 0;
-
-        forAll(N_, i)
-        {
-            result += pow(bin[i], k)*N_[i];
-        }
-
-        if (Pstream::master())
-        {
-            file() << tab << result;
-        }
-    }
-
-    if (Pstream::master())
-    {
-        file() << endl;
-    }
-}
-
-
-void Foam::functionObjects::sizeDistribution::writeStdDev()
-{
-    logFiles::write();
-
-    Log << "    writing standard deviation of size distribution."
-        << endl;
-
-    if (Pstream::master())
-    {
-        writeTime(file());
-    }
-
-    scalar stdDev = 0;
-    scalar mean = 0;
-    scalar var = 0;
-
-    const scalarField& bin = this->bin();
-
-    if (sum(N_) != 0)
-    {
-        if (geometric_)
-        {
-            mean = exp(sum(Foam::log(bin)*N_/sum(N_)));
-
-            var =
-                sum(sqr(Foam::log(bin) - Foam::log(mean))
-               *N_/sum(N_));
-
-            stdDev = exp(sqrt(var));
-        }
-        else
-        {
-            mean = sum(bin*N_/sum(N_));
-
-            var = sum(sqr(bin - mean)*N_/sum(N_));
-
-            stdDev = sqrt(var);
-        }
-    }
-
-    if (Pstream::master())
-    {
-        file() << tab << stdDev << tab << mean << tab << var << endl;
-    }
-}
-
-
-void Foam::functionObjects::sizeDistribution::writeDistribution()
-{
-    scalarField result(N_);
-
-    const scalarField& bin = this->bin();
-
-    switch (functionType_)
-    {
-        case ftNumber:
-        {
-            Log << "    writing number distribution. "
-                << endl;
-
-            break;
-        }
-        case ftVolume:
-        {
-            Log << "    writing volume distribution. "
-                << endl;
-
-            result *= V_;
-
-            break;
-        }
-        default:
-        {
-            break;
-        }
-
-    }
-
-    if (normalise_)
-    {
-        if(sum(result) != 0)
-        {
-            result /= sum(result);
-        }
-
-    }
-
-    if (densityFunction_)
-    {
-        List<scalar> bndrs(N_.size() + 1);
-
-        bndrs.first() = bin.first();
-        bndrs.last() = bin.last();
-
-        for (label i = 1; i < N_.size(); i++)
-        {
-            bndrs[i] = (bin[i]+ bin[i-1])/2.0;
-        }
-
-        forAll(result, i)
-        {
-            if (geometric_)
-            {
-                result[i] /=
-                    (Foam::log(bndrs[i+1]) - Foam::log(bndrs[i]));
-            }
-            else
-            {
-                result[i] /= (bndrs[i+1] - bndrs[i]);
-            }
-        }
-    }
-
-    if (Pstream::master())
-    {
-        formatterPtr_->write
 (
-            file_.baseTimeDir(),
+    const Foam::diameterModels::sizeGroup& fi,
-            name(),
+    const coordinateType& cType
-            coordSet(true, "volume", V_),
-            "area",
-            a_,
-            "diameter",
-            d_,
-            word(functionTypeNames_[functionType_])
-          + (densityFunction_ ? "Density" : "Concentration"),
-            result
-        );
-    }
-}
-
-
-void Foam::functionObjects::sizeDistribution::writeFileHeader
-(
-    const label i
 )
 {
-    volRegion::writeFileHeader(*this, file());
+    scalar averageCoordinateValue(Zero);
 
-    writeHeaderValue
+    switch (cType)
+    {
+        case coordinateType::volume:
+        {
+            averageCoordinateValue = fi.x().value();
+
+            break;
+        }
+        case coordinateType::area:
+        {
+            averageCoordinateValue =
+                weightedAverage(fi.a(), fi);
+
+            break;
+        }
+        case coordinateType::diameter:
+        {
+            averageCoordinateValue =
+                weightedAverage(fi.d(), fi);
+
+            break;
+        }
+        case coordinateType::projectedAreaDiameter:
+        {
+            averageCoordinateValue =
+                weightedAverage(sqrt(fi.a()/pi), fi);
+
+            break;
+        }
+    }
+
+    return averageCoordinateValue;
+}
+
+
+Foam::scalar Foam::functionObjects::sizeDistribution::weightedAverage
 (
-        file(),
+    const Foam::scalarField& fld,
-        "Coordinate",
+    const Foam::diameterModels::sizeGroup& fi
-        word(coordinateTypeNames_[coordinateType_])
+)
-    );
+{
+    scalar weightedAverage(Zero);
 
-    word str("Time");
+    switch (weightType_)
+    {
+        case weightType::numberConcentration:
+        {
+            scalarField Ni(filterField(fi*fi.phase()/fi.x().value()));
 
-    switch (functionType_)
+            if (gSum(Ni) == 0)
             {
-        case ftMoments:
+                weightedAverage =
+                    gSum(filterField(mesh_.V()*fld))/this->V();
+            }
+            else
             {
-            for (label k = 0; k <= maxOrder_; k++)
+                weightedAverage =
-            {
+                    gSum(Ni*filterField(fld))/gSum(Ni);
-                str += (" k=" + std::to_string(k));
             }
 
             break;
         }
-
+        case weightType::volumeConcentration:
-        case ftStdDev:
         {
-            str += " standardDeviation mean variance";
+            scalarField Vi(filterField(fi*fi.phase()));
+
+            if (gSum(Vi) == 0)
+            {
+                weightedAverage =
+                    gSum(filterField(mesh_.V()*fld))/this->V();
+            }
+            else
+            {
+                weightedAverage =
+                    gSum(Vi*filterField(fld))/gSum(Vi);
+            }
 
             break;
         }
-
+        case weightType::areaConcentration:
-        default:
         {
+            scalarField Ai(filterField(fi.a().ref()*fi.phase()));
+
+            if (gSum(Ai) == 0)
+            {
+                weightedAverage =
+                    gSum(filterField(mesh_.V()*fld))/this->V();
+            }
+            else
+            {
+                weightedAverage =
+                    gSum(Ai*filterField(fld))/gSum(Ai);
+            }
+
+            break;
+        }
+        case weightType::cellVolume:
+        {
+            weightedAverage =
+                gSum(filterField(mesh_.V()*fld))/this->V();
+
             break;
         }
     }
 
-    writeCommented(file(), str);
+    return weightedAverage;
-
-    file() << endl;
 }
 
 
@@ -362,9 +341,8 @@ Foam::functionObjects::sizeDistribution::sizeDistribution
 :
     fvMeshFunctionObject(name, runTime, dict),
     volRegion(fvMeshFunctionObject::mesh_, dict),
-    logFiles(obr_, name),
-    mesh_(fvMeshFunctionObject::mesh_),
     file_(obr_, name),
+    mesh_(fvMeshFunctionObject::mesh_),
     popBal_
     (
         obr_.lookupObject<Foam::diameterModels::populationBalanceModel>
@@ -374,10 +352,26 @@ Foam::functionObjects::sizeDistribution::sizeDistribution
     ),
     functionType_(functionTypeNames_.read(dict.lookup("functionType"))),
     coordinateType_(coordinateTypeNames_.read(dict.lookup("coordinateType"))),
-    N_(popBal_.sizeGroups().size(), 0),
+    allCoordinates_
-    V_(popBal_.sizeGroups().size(), 0),
+    (
-    a_(popBal_.sizeGroups().size(), 0),
+        dict.lookupOrDefault<Switch>("allCoordinates", false)
-    d_(popBal_.sizeGroups().size(), 0)
+    ),
+    normalise_(dict.lookupOrDefault<Switch>("normalise", false)),
+    logTransform_
+    (
+        dict.lookupOrDefaultBackwardsCompatible<Switch>
+        (
+            {"logTransform", "geometric"},
+            false
+        )
+    ),
+    weightType_
+    (
+        dict.found("weightType")
+      ? weightTypeNames_.read(dict.lookup("weightType"))
+      : weightType::numberConcentration
+    ),
+    formatterPtr_(nullptr)
 {
     read(dict);
 }
@@ -393,17 +387,12 @@ Foam::functionObjects::sizeDistribution::~sizeDistribution()
 
 bool Foam::functionObjects::sizeDistribution::read(const dictionary& dict)
 {
+    Log << type() << " " << name() << ":" << nl;
+
     fvMeshFunctionObject::read(dict);
 
-    normalise_ = dict.lookupOrDefault<Switch>("normalise", false);
-    densityFunction_ = dict.lookupOrDefault<Switch>("densityFunction", false);
-    geometric_ = dict.lookupOrDefault<Switch>("geometric", false);
-    maxOrder_ = dict.lookupOrDefault("maxOrder", 3);
-
     formatterPtr_ = setWriter::New(dict.lookup("setFormat"), dict);
 
-    resetName(name());
-
     return false;
 }
 
@@ -414,43 +403,255 @@ bool Foam::functionObjects::sizeDistribution::execute()
 }
 
 
-bool Foam::functionObjects::sizeDistribution::end()
-{
-    return true;
-}
-
-
 bool Foam::functionObjects::sizeDistribution::write()
 {
     Log << type() << " " << name() << " write:" << nl;
 
-    correctVolAverages();
+    const UPtrList<diameterModels::sizeGroup>& sizeGroups =
+        popBal_.sizeGroups();
+
+    scalarField coordinateValues(sizeGroups.size());
+    scalarField boundaryValues(sizeGroups.size() + 1);
+    scalarField resultValues(sizeGroups.size());
+
+    forAll(sizeGroups, i)
+    {
+        const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+        coordinateValues[i] = averageCoordinateValue(fi, coordinateType_);
+    }
+
+    if
+    (
+        functionType_ == functionType::numberDensity
+     || functionType_ == functionType::volumeDensity
+     || functionType_ == functionType::areaDensity
+    )
+    {
+        if (logTransform_)
+        {
+            boundaryValues.first() = Foam::log(coordinateValues.first());
+            boundaryValues.last() = Foam::log(coordinateValues.last());
+
+            for (label i = 1; i < boundaryValues.size() - 1; i++)
+            {
+                boundaryValues[i] =
+                    0.5
+                   *(
+                        Foam::log(coordinateValues[i])
+                      + Foam::log(coordinateValues[i-1])
+                    );
+            }
+        }
+        else
+        {
+            boundaryValues.first() = coordinateValues.first();
+            boundaryValues.last() = coordinateValues.last();
+
+            for (label i = 1; i < boundaryValues.size() - 1; i++)
+            {
+                boundaryValues[i] =
+                    (coordinateValues[i] + coordinateValues[i-1])/2;
+            }
+        }
+    }
 
     switch (functionType_)
     {
-        case ftMoments:
+        case functionType::numberConcentration:
         {
-            writeMoments();
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                resultValues[i] =
+                    gSum(filterField(mesh_.V()*fi*fi.phase()/fi.x()))*this->V();
+            }
+
+            if (normalise_ && sum(resultValues) != 0)
+            {
+                resultValues /= sum(resultValues);
+            }
 
             break;
         }
-
+        case functionType::numberDensity:
-        case ftStdDev:
         {
-            writeStdDev();
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                resultValues[i] =
+                    gSum(filterField(mesh_.V()*fi*fi.phase()/fi.x()))*this->V();
+            }
+
+            if (normalise_ && sum(resultValues) != 0)
+            {
+                resultValues /= sum(resultValues);
+            }
+
+            forAll(resultValues, i)
+            {
+                resultValues[i] /= (boundaryValues[i+1] - boundaryValues[i]);
+            }
 
             break;
         }
-
+        case functionType::volumeConcentration:
-        default:
         {
-            writeDistribution();
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                resultValues[i] =
+                    gSum(filterField(mesh_.V()*fi*fi.phase()))*this->V();
+            }
+
+            if (normalise_ && sum(resultValues) != 0)
+            {
+                resultValues /= sum(resultValues);
+            }
+
+            break;
+        }
+        case functionType::volumeDensity:
+        {
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                resultValues[i] =
+                    gSum(filterField(mesh_.V()*fi*fi.phase()))*this->V();
+            }
+
+            if (normalise_ && sum(resultValues) != 0)
+            {
+                resultValues /= sum(resultValues);
+            }
+
+            forAll(resultValues, i)
+            {
+                resultValues[i] /= (boundaryValues[i+1] - boundaryValues[i]);
+            }
+
+            break;
+        }
+        case functionType::areaConcentration:
+        {
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                resultValues[i] =
+                    gSum
+                    (
+                        filterField(mesh_.V()*fi.a().ref()*fi*fi.phase()/fi.x())
+                    )
+                   *this->V();
+            }
+
+            if (normalise_ && sum(resultValues) != 0)
+            {
+                resultValues /= sum(resultValues);
+            }
+
+            break;
+        }
+        case functionType::areaDensity:
+        {
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                resultValues[i] =
+                    gSum
+                    (
+                        filterField(mesh_.V()*fi.a().ref()*fi*fi.phase()/fi.x())
+                    )
+                   *this->V();
+            }
+
+            if (normalise_ && sum(resultValues) != 0)
+            {
+                resultValues /= sum(resultValues);
+            }
+
+            forAll(resultValues, i)
+            {
+                resultValues[i] /= (boundaryValues[i+1] - boundaryValues[i]);
+            }
 
             break;
         }
     }
 
-    Log << endl;
+
+    if (allCoordinates_)
+    {
+        wordList otherCoordinateSymbolicNames(coordinateTypeNames_.size());
+        PtrList<scalarField> otherCoordinateValues(coordinateTypeNames_.size());
+        typedef NamedEnum<coordinateType, 4> namedEnumCoordinateType;
+
+        forAllConstIter(namedEnumCoordinateType, coordinateTypeNames_, iter)
+        {
+            const coordinateType cType = coordinateTypeNames_[iter.key()];
+
+            otherCoordinateSymbolicNames[cType] =
+                coordinateTypeSymbolicName(cType);
+
+            otherCoordinateValues.set
+            (
+                cType,
+                new scalarField(popBal_.sizeGroups().size())
+            );
+
+            forAll(sizeGroups, i)
+            {
+                const diameterModels::sizeGroup& fi = sizeGroups[i];
+
+                otherCoordinateValues[cType][i] =
+                    averageCoordinateValue(fi, cType);
+            }
+        }
+
+        if (Pstream::master())
+        {
+            formatterPtr_->write
+            (
+                file_.baseTimeDir(),
+                name(),
+                coordSet
+                (
+                    true,
+                    coordinateTypeSymbolicName(coordinateType_),
+                    coordinateValues
+                ),
+                functionTypeSymbolicName(),
+                resultValues,
+                otherCoordinateSymbolicNames,
+                otherCoordinateValues
+            );
+        }
+    }
+    else
+    {
+        if (Pstream::master())
+        {
+            formatterPtr_->write
+            (
+                file_.baseTimeDir(),
+                name(),
+                coordSet
+                (
+                    true,
+                    coordinateTypeSymbolicName(coordinateType_),
+                    coordinateValues
+                ),
+                functionTypeSymbolicName(),
+                resultValues
+            );
+        }
+    }
 
     return true;
 }

applications/solvers/multiphase/multiphaseEulerFoam/functionObjects/sizeDistribution/sizeDistribution.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,14 +25,11 @@ Class
     Foam::functionObjects::sizeDistribution
 
 Description
-    This function object calculates and outputs volume-averaged information
+    Writes out the size distribution computed with multiphaseEulerFoam for the
-    about the size distribution of the dispersed phase, such as the number
+    entire domain or a volume region. Requires solver post-processing.
-    density function or its moments. It is designed to be used exclusively with
+
-    the population balance modeling functionality of the multiphaseEulerFoam
+    The following function object specification for example returns the volume-
-    solver. It can be applied to a specific cellZone or the entire domain. The
+    based number density function:
-    function type determines whether the density function and its moments are
-    based on the number of dispersed phase elements in a size group or their
-    total volume.
 
     Example of function object specification:
     \verbatim
@@ -40,27 +37,30 @@ Description
     {
         type                sizeDistribution;
         libs                ("libmultiphaseEulerFoamFunctionObjects.so");
-        ...
+        writeControl        writeTime;
         populationBalance   bubbles;
-        regionType          cellZone;
+        functionType        numberDensity;
-        name                zone0;
-        functionType        number;
         coordinateType      volume;
-        densityFunction     yes;
+        setFormat           raw;
     }
     \endverbatim
 
 Usage
     \table
-        Property     | Description               | Required | Default value
+        Property          | Description                    | Required  | Default
-        type         | type name: sizeDistribution | yes    |
+        populationBalance | population balance name        | yes       |
-        populationBalance | corresponding populationBalance | yes |
+        functionType      | function type                  | yes       |
-        functionType | number/volume/moments/stdDev | yes   |
+        coordinateType    | particle property              | yes       |
-        coordinateType | used for density/moment calculation | yes |
+        allCoordinates    | write all coordinate values    | no        | false
-        normalise    | normalise concentrations  | no       | no
+        normalise         | divide by total concentration  | no        | false
-        densityFunction | compute densityFunction | no      | no
+        logTransform      | class width based on log of coordinate\\
-        logBased     | use log of coordinate for density | no | no
+                                                           | no        | false
-        maxOrder     | maxim order of moment output | no    | 3
+        weightType        | weighting in case of field-dependent particle\\
+                            properties                     | no\\
+                                                           | numberConcentration
+        regionType        | cellZone or all                | no        | all
+        name              | name of cellZone if required   | no        |
+        setFormat         | output format                  | yes       |
     \endtable
 
 SourceFiles
@@ -73,7 +73,6 @@ SourceFiles
 
 #include "fvMeshFunctionObject.H"
 #include "volRegion.H"
-#include "logFiles.H"
 #include "populationBalanceModel.H"
 #include "writeFile.H"
 #include "setWriter.H"
@@ -92,8 +91,7 @@ namespace functionObjects
 class sizeDistribution
 :
     public fvMeshFunctionObject,
-    public volRegion,
+    public volRegion
-    public logFiles
 {
 public:
 
@@ -102,111 +100,101 @@ public:
         //- Function type enumeration
         enum functionType
         {
-            ftMoments,
+            numberConcentration,
-            ftStdDev,
+            numberDensity,
-            ftNumber,
+            volumeConcentration,
-            ftVolume
+            volumeDensity,
+            areaConcentration,
+            areaDensity
         };
 
-        //- Ordinate type names
+        //- Function type names
-        static const NamedEnum<functionType, 4> functionTypeNames_;
+        static const NamedEnum<functionType, 6> functionTypeNames_;
 
         //- Coordinate type enumeration
         enum coordinateType
         {
-            ctVolume,
+            volume,
-            ctArea,
+            area,
-            ctDiameter,
+            diameter,
-            ctProjectedAreaDiameter
+            projectedAreaDiameter
         };
 
         //- Coordinate type names
         static const NamedEnum<coordinateType, 4> coordinateTypeNames_;
 
+        //- Enumeration for the weight types
+        enum class weightType
+        {
+            numberConcentration,
+            volumeConcentration,
+            areaConcentration,
+            cellVolume
+        };
 
-protected:
+        //- Names of the weight types
+        static const NamedEnum<weightType, 4> weightTypeNames_;
 
-    // Protected Data
 
-        //- Reference to fvMesh
+private:
-        const fvMesh& mesh_;
 
-        //- File containing data for all functionTypes except moments
+    // Private Data
+
+        //- Write file
         writeFile file_;
 
-        //- Output formatter
+        //- Reference to mesh
-        autoPtr<setWriter> formatterPtr_;
+        const fvMesh& mesh_;
 
-        //- Reference to populationBalanceModel
+        //- Reference to population balance
         const Foam::diameterModels::populationBalanceModel& popBal_;
 
-        //- Function to evaluate
+        //- Function type
         functionType functionType_;
 
-        //- Abscissa type
+        //- Coordinate type
         coordinateType coordinateType_;
 
-        //- List of volume-averaged number concentrations
+        //- Add values for all coordinate types to output
-        scalarField N_;
+        Switch allCoordinates_;
 
-        //- ???
+        //- Normalise result through division by sum
-        scalarField V_;
-
-        //- List of volume-averaged surface areas
-        scalarField a_;
-
-        //- List of volume-averaged diameters
-        scalarField d_;
-
-        //- Normalise number/volume concentrations
         Switch normalise_;
 
-        //- Determines whether density function is calculated
+        //- Log transform
-        Switch densityFunction_;
+        Switch logTransform_;
 
-        //- Geometric standard deviation/density function
+        //- Weight types, relevant if particle properties are field dependent
-        Switch geometric_;
+        weightType weightType_;
 
-        //- Highest moment order
+        //- Set formatter
-        label maxOrder_;
+        autoPtr<setWriter> formatterPtr_;
 
 
-    // Protected Member Functions
+    // Private Member Functions
 
-        //- Filter field according to cellIds
+        //- Function type symbolic name for shorter file header
+        word functionTypeSymbolicName();
+
+        //- Coordinate type symbolic name for shorter file header
+        word coordinateTypeSymbolicName(const coordinateType& cType);
+
+        //- Filter a field according to cellIds
         tmp<scalarField> filterField(const scalarField& field) const;
 
-        //- Bin component used according to chosen coordinate type
+        //- Field averaged coordinate value
-        inline const scalarField& bin() const
+        scalar averageCoordinateValue
-        {
+        (
-            switch (coordinateType_)
+            const diameterModels::sizeGroup& fi,
-            {
+            const coordinateType& cType
-                case ctVolume:
+        );
-                    return V_;
-                case ctArea:
-                    return a_;
-                case ctDiameter:
-                    return d_;
-                case ctProjectedAreaDiameter:
-                    return d_;
-            }
-            return scalarField::null();
-        }
 
-        //- Correct volume averages
+        //- Weighted average
-        void correctVolAverages();
+        scalar weightedAverage
-
+        (
-        //- Write moments
+            const scalarField& fld,
-        void writeMoments();
+            const diameterModels::sizeGroup& fi
-
+        );
-        //- Write standard deviation
-        void writeStdDev();
-
-        //- Write distribution
-        void writeDistribution();
-
-        //- Output file header information for functionType moments
-        virtual void writeFileHeader(const label i);
 
 
 public:
@@ -235,21 +223,18 @@ public:
 
     // Member Functions
 
-        //- Read the sizeDistribution data
-        virtual bool read(const dictionary&);
-
         //- Return the list of fields required
         virtual wordList fields() const
         {
             return wordList::null();
         }
 
+        //- Read the sizeDistribution data
+        virtual bool read(const dictionary&);
+
         //- Execute, currently does nothing
         virtual bool execute();
 
-        //- Execute at the final time-loop, currently does nothing
-        virtual bool end();
-
         //- Calculate and write the size distribution
         virtual bool write();
 

etc/caseDicts/postProcessing/multiphase/moments

@@ -0,0 +1,34 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Calculates and writes out integral (integer moments) or mean properties
+    (mean, variance, standard deviation) of a size distribution computed with
+    multiphaseEulerFoam. Requires solver post-processing.
+
+\*---------------------------------------------------------------------------*/
+
+type                moments;
+libs                ("libmultiphaseEulerFoamFunctionObjects.so");
+
+populationBalance   <populationBalance>;
+momentType          <moment>;             // integerMoment, mean, variance,
+                                          // stdDev
+coordinateType      <coordinate>;         // volume, area, diameter
+weightType          numberConcentration;  // volumeConcentration,
+                                          // areaConcentration
+                                          // defaults to numberConcentration
+order               0;                    // relevant for integer moments only
+mean                arithmetic;           // geometric
+                                          // relevant for non-integer moments,
+                                          // defaults to arithmetic
+
+executeControl      timeStep;
+writeControl        writeTime;
+
+
+// ************************************************************************* //

etc/caseDicts/postProcessing/multiphase/sizeDistribution

@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+Description
+    Writes out the size distribution computed with multiphaseEulerFoam for the
+    entire domain or a volume region. Requires solver post-processing.
+
+\*---------------------------------------------------------------------------*/
+
+type                sizeDistribution;
+libs                ("libmultiphaseEulerFoamFunctionObjects.so");
+
+populationBalance   <populationBalance>;
+functionType        <function>;           // numberConcentration, numberDensity
+                                          // volumeConcentration, volumeDensity
+                                          // areaConcentration, areaDensity
+coordinateType      <coordinate>;         // volume, area, diameter,
+                                          // projectedAreaDiameter
+allCoordinates      false;                // defaults to false
+normalise           false;                // defaults to false
+logTransform        false;                // defaults to false, only relevant
+                                          // for density functions
+weightType          numberConcentration;  // volumeConcentration,
+                                          // areaConcentration, cellVolume
+                                          // relevant for field-dependent
+                                          // particle properties, defaults to
+                                          // numberConcentration
+regionType          all;                  // cellZone
+                                          // defaults to all
+name                cellZoneName;         // relevant for regionType all
+
+
+setFormat           raw;
+writeControl        writeTime;
+
+
+// ************************************************************************* //

src/sampling/sampledSet/writers/setWriter.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,24 +69,26 @@ public:
 
         //- Helper for variadic write
         template<class Types, class Type>
-        static inline void setTypeValueSet
+        static inline void appendTypeValueSet
         (
             UPtrList<const Field<Types>>& valueSets,
-            const label valueSeti,
-            const Field<Type>& valueSet
-        )
-        {}
-
-        //- Helper for variadic write
-        template<class Type>
-        static inline void setTypeValueSet
-        (
-            UPtrList<const Field<Type>>& valueSets,
-            const label valueSeti,
             const Field<Type>& valueSet
         )
         {
-            valueSets.set(valueSeti, &valueSet);
+            valueSets.resize(valueSets.size() + 1);
+            valueSets.set(valueSets.size() - 1, nullptr);
+        }
+
+        //- Helper for variadic write
+        template<class Type>
+        static inline void appendTypeValueSet
+        (
+            UPtrList<const Field<Type>>& valueSets,
+            const Field<Type>& valueSet
+        )
+        {
+            valueSets.resize(valueSets.size() + 1);
+            valueSets.set(valueSets.size() - 1, &valueSet);
         }
 
         //- Helper for variadic write
@@ -114,14 +116,11 @@ public:
             Args& ... args
         )
         {
-            const label valueSeti =
+            valueSetNames.append(valueSetName);
-                valueSetNames.size() - 1 - sizeof...(Args)/2;
+            #define AppendTypeValueSet(Type, nullArg) \
-
+                appendTypeValueSet(Type##ValueSets, valueSet);
-            valueSetNames[valueSeti] = valueSetName;
+            FOR_ALL_FIELD_TYPES(AppendTypeValueSet);
-            #define SetTypeValueSet(Type, nullArg) \
+            #undef AppendTypeValueSet
-                setTypeValueSet(Type##ValueSets, valueSeti, valueSet);
-            FOR_ALL_FIELD_TYPES(SetTypeValueSet);
-            #undef SetTypeValueSet
 
             unpackTypeValueSets
             (
@@ -134,6 +133,68 @@ public:
             );
         }
 
+        //- Helper for variadic write
+        template<class Type, class ... Args>
+        static inline void unpackTypeValueSets
+        (
+            wordList& valueSetNames
+            #define TypeValueSetsNonConstArg(Type, nullArg) \
+                , UPtrList<const Field<Type>>& Type##ValueSets
+            FOR_ALL_FIELD_TYPES(TypeValueSetsNonConstArg),
+            #undef TypeValueSetsNonConstArg
+            const wordList& valueSetNamesPart,
+            const UPtrList<const Field<Type>>& valueSetsPart,
+            Args& ... args
+        )
+        {
+            forAll(valueSetNamesPart, i)
+            {
+                valueSetNames.append(valueSetNamesPart[i]);
+                #define AppendTypeValueSet(Type, nullArg) \
+                    appendTypeValueSet(Type##ValueSets, valueSetsPart[i]);
+                FOR_ALL_FIELD_TYPES(AppendTypeValueSet);
+                #undef AppendTypeValueSet
+            }
+
+            unpackTypeValueSets
+            (
+                valueSetNames
+                #define TypeValueSetsParameter(Type, nullArg) \
+                    , Type##ValueSets
+                FOR_ALL_FIELD_TYPES(TypeValueSetsParameter),
+                #undef TypeValueSetsParameter
+                args ...
+            );
+        }
+
+        //- Helper for variadic write
+        template<class Type, class ... Args>
+        static inline void unpackTypeValueSets
+        (
+            wordList& valueSetNames
+            #define TypeValueSetsNonConstArg(Type, nullArg) \
+                , UPtrList<const Field<Type>>& Type##ValueSets
+            FOR_ALL_FIELD_TYPES(TypeValueSetsNonConstArg),
+            #undef TypeValueSetsNonConstArg
+            const wordList& valueSetNamesPart,
+            const PtrList<Field<Type>>& valueSetsPart,
+            Args& ... args
+        )
+        {
+            unpackTypeValueSets
+            (
+                valueSetNames
+                #define TypeValueSetsParameter(Type, nullArg) \
+                    , Type##ValueSets
+                FOR_ALL_FIELD_TYPES(TypeValueSetsParameter),
+                #undef TypeValueSetsParameter
+                valueSetNamesPart,
+                reinterpret_cast<const UPtrList<const Field<Type>>&>
+                (valueSetsPart),
+                args ...
+            );
+        }
+
 
 protected:
 
@@ -345,11 +406,9 @@ public:
             const Args& ... args
         ) const
         {
-            const label nFields = sizeof...(Args)/2;
+            wordList valueSetNames;
-
-            wordList valueSetNames(nFields);
             #define DeclareTypeValueSets(Type, nullArg) \
-                UPtrList<const Field<Type>> Type##ValueSets(nFields);
+                UPtrList<const Field<Type>> Type##ValueSets;
             FOR_ALL_FIELD_TYPES(DeclareTypeValueSets);
             #undef DeclareTypeValueSets
 

test/multiphase/multiphaseEulerFoam/populationBalance/binaryBreakup/system/controlDict

@@ -52,21 +52,13 @@ maxDeltaT       1;
 
 functions
 {
-    numberDensity
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        functionType=numberDensity,
-
+        coordinateType=volume,
-        writeControl        outputTime;
+        funcName=numberDensity
-        writeInterval       1;
+    )
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        functionType        number;
-        coordinateType      volume;
-        densityFunction     yes;
-    }
 }
 
 // ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/binaryBreakup/validation/createGraphs

@@ -31,9 +31,9 @@ gnuplot<<EOF
 
     plot '<printf "1 0.05\n1 1e-5\n"' w l lc rgb 'black' t 'Init. Cond.',\
     n(x,10.0) dt 2 lc rgb 'black' t 'Ziff and McGrady (1985)',\
-    numberDensity every ::0::19 u 1:4 w p pt 1 lc rgb 'black' t 'air1',\
+    numberDensity every ::0::19 u 1:2 w p pt 1 lc rgb 'black' t 'air1',\
-    numberDensity every ::20::24 u 1:4 w p pt 5 lc rgb 'black' t 'air2',\
+    numberDensity every ::20::24 u 1:2 w p pt 5 lc rgb 'black' t 'air2',\
-    numberDensity every ::25::28 u 1:4 w p pt 9 lc rgb 'black' t 'air3'
+    numberDensity every ::25::28 u 1:2 w p pt 9 lc rgb 'black' t 'air3'
 EOF
 
 #------------------------------------------------------------------------------

test/multiphase/multiphaseEulerFoam/populationBalance/breakup/system/controlDict

@@ -52,21 +52,13 @@ maxDeltaT       1;
 
 functions
 {
-    numberDensity
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        functionType=numberDensity,
-
+        coordinateType=volume,
-        writeControl        outputTime;
+        funcName=numberDensity
-        writeInterval       1;
+    )
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        functionType        number;
-        coordinateType      volume;
-        densityFunction     yes;
-    }
 }
 
 // ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/breakup/validation/createGraphs

@@ -31,9 +31,9 @@ gnuplot<<EOF
 
     plot '<printf "1 0.05\n1 1e-5\n"' w l lc rgb 'black' t 'Init. Cond.',\
     n(x,10.0) dt 2 lc rgb 'black' t 'Ziff and McGrady (1985)',\
-    numberDensity every ::0::19 u 1:4 w p pt 1 lc rgb 'black' t 'air1',\
+    numberDensity every ::0::19 u 1:2 w p pt 1 lc rgb 'black' t 'air1',\
-    numberDensity every ::20::24 u 1:4 w p pt 5 lc rgb 'black' t 'air2',\
+    numberDensity every ::20::24 u 1:2 w p pt 5 lc rgb 'black' t 'air2',\
-    numberDensity every ::25::28 u 1:4 w p pt 9 lc rgb 'black' t 'air3'
+    numberDensity every ::25::28 u 1:2 w p pt 9 lc rgb 'black' t 'air3'
 EOF
 
 #------------------------------------------------------------------------------

test/multiphase/multiphaseEulerFoam/populationBalance/coalescence/system/controlDict

@@ -52,21 +52,13 @@ maxDeltaT       1;
 
 functions
 {
-    numberDensity
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        functionType=numberDensity,
-
+        coordinateType=volume,
-        writeControl        outputTime;
+        funcName=numberDensity
-        writeInterval       1;
+    )
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        functionType        number;
-        coordinateType      volume;
-        densityFunction     yes;
-    }
 }
 
 // ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/coalescence/validation/createGraphs

@@ -51,9 +51,9 @@ gnuplot<<EOF
 
     plot n0(x) ls -1 t 'Initial Condition',\
     n(x,10.0) dt 2 lc rgb 'black' t 'Scott (1968)',\
-    numberDensity every ::0::46 u 1:4 w p pt 1 lc rgb 'black' t 'air1',\
+    numberDensity every ::0::46 u 1:2 w p pt 1 lc rgb 'black' t 'air1',\
-    numberDensity every ::47::55 u 1:4 w p pt 5 lc rgb 'black' t 'air2',\
+    numberDensity every ::47::55 u 1:2 w p pt 5 lc rgb 'black' t 'air2',\
-    numberDensity every ::56::70 u 1:4 w p pt 9 lc rgb 'black' t 'air3'
+    numberDensity every ::56::70 u 1:2 w p pt 9 lc rgb 'black' t 'air3'
 EOF
 
 #------------------------------------------------------------------------------

test/multiphase/multiphaseEulerFoam/populationBalance/drift/Allclean

@@ -6,6 +6,6 @@ cd "${0%/*}" || exit 1
 # Source clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
-cleanCase && rm -f moments.eps numberDensity.eps
+cleanCase && rm -rf *.eps 0/d.air 0/uniform/ 0/integerMoment*
 
 #------------------------------------------------------------------------------

test/multiphase/multiphaseEulerFoam/populationBalance/drift/system/controlDict

@@ -28,7 +28,7 @@ deltaT          0.01;
 
 writeControl    runTime;
 
-writeInterval   6;
+writeInterval   0.2;
 
 purgeWrite      0;
 
@@ -52,37 +52,38 @@ maxDeltaT       1;
 
 functions
 {
-    numberDensity
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        functionType=numberDensity,
+        coordinateType=volume,
+        funcName=numberDensity
+    )
 
-        writeControl        outputTime;
+    #includeFunc moments
-        writeInterval       1;
+    (
+        populationBalance=bubbles,
+        momentType=integerMoment,
+        coordinateType=volume,
+        order=1
+    )
 
-        setFormat           raw;
+    #includeFunc moments
+    (
+        populationBalance=bubbles,
+        momentType=integerMoment,
+        coordinateType=volume,
+        order=0
+    )
 
-        populationBalance   bubbles;
+    #includeFunc probes
-        functionType        number;
+    (
-        coordinateType      volume;
+        funcName=probes,
-        densityFunction     yes;
+        points=((0.5 0.5 0.5)),
-    }
+        integerMoment0(N,v).bubbles,
-
+        integerMoment1(N,v).bubbles,
-    moments
+        writeControl=writeTime
-    {
+    )
-        type                sizeDistribution;
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
-
-        writeControl        runTime;
-        writeInterval       0.1;
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        functionType        moments;
-        coordinateType      volume;
-        maxOrder            1;
-    }
 }
 
 // ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/drift/validation/createGraphs

@@ -23,7 +23,7 @@ gnuplot<<EOF
 
     plot '<printf "1.876749927 0\n1.876749927 0.01\n4.144510443 0.01\n4.144510443 0\n"' w l dt 2 lc rgb 'black' t 'Initial condition',\
     '<printf "7.957828406 0\n7.957828406 0.01\n10.14451044 0.01\n10.14451044 0\n"' w l lc rgb 'black' t 'Exact',\
-    '../postProcessing/numberDensity/6/numberDensity.xy' u 1:4 w lp ls 4 lc rgb 'black' t 'Numerical'
+    '../postProcessing/numberDensity/6/numberDensity.xy' u 1:2 w lp ls 4 lc rgb 'black' t 'Numerical'
 
     set output '../moments.eps'
 
@@ -36,9 +36,9 @@ gnuplot<<EOF
     set yrange [0.5:3.5]
     set key top left
 
-    plot '../postProcessing/moments/0/moments.dat' every 2 u 1:(\$2/0.02186682) w p pt 5 lc rgb 'black' t 'j=0 (Numbers), numerical',\
+    plot '../postProcessing/probes/0/integerMoment0(N,v).bubbles' u 1:(\$2/0.02186682) w p pt 5 lc rgb 'black' t 'j=0 (Numbers), numerical',\
     '<printf "0 1\n6 1\n"' w l lc rgb 'black' t 'j=0 (Numbers), exact',\
-    '../postProcessing/moments/0/moments.dat' u 1:(\$3/0.066553351) every 2 w p pt 6 lc rgb 'black' t 'j=1 (Volume), numerical',\
+    '../postProcessing/probes/0/integerMoment1(N,v).bubbles' u 1:(\$2/0.066553351) w p pt 6 lc rgb 'black' t 'j=1 (Volume), numerical',\
     '<printf "0 1\n6 3\n"' w l dt 2 lc rgb 'black' t 'j=1 (Volume), exact'
 EOF
 

test/multiphase/multiphaseEulerFoam/populationBalance/isothermalGrowth/validation/d.bubbles

@@ -1,299 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  sprint
-     \\/     M anipulation  |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    format      ascii;
-    class       volScalarField;
-    location    "5";
-    object      d.bubbles;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 0 0 0 0 0];
-
-internalField   nonuniform List<scalar>
-50
-(
-0.00300040563
-0.00300110488
-0.00300180018
-0.00300249358
-0.00300318666
-0.00300387818
-0.00300457055
-0.00300526257
-0.00300595539
-0.00300664879
-0.00300734302
-0.00300803811
-0.00300873419
-0.00300943126
-0.00301012937
-0.00301082851
-0.00301152868
-0.00301222987
-0.00301293204
-0.00301363515
-0.00301433916
-0.00301504402
-0.00301574966
-0.00301645602
-0.00301716301
-0.00301787056
-0.00301857858
-0.00301928697
-0.00301999563
-0.00302070448
-0.00302141342
-0.00302212238
-0.00302283129
-0.00302354009
-0.00302424876
-0.00302495726
-0.0030256656
-0.00302637377
-0.00302708181
-0.00302778974
-0.0030284976
-0.00302920541
-0.00302991328
-0.00303062095
-0.00303132989
-0.00303203401
-0.00303275902
-0.00303342196
-0.00303447928
-0.00303514582
-)
-;
-
-boundaryField
-{
-    inlet
-    {
-        type            calculated;
-        value           uniform 0.00300000033;
-    }
-    outlet
-    {
-        type            calculated;
-        value           uniform 0.00303514582;
-    }
-    walls
-    {
-        type            calculated;
-        value           nonuniform List<scalar>
-200
-(
-0.00300040563
-0.00300110488
-0.00300180018
-0.00300249358
-0.00300318666
-0.00300387818
-0.00300457055
-0.00300526257
-0.00300595539
-0.00300664879
-0.00300734302
-0.00300803811
-0.00300873419
-0.00300943126
-0.00301012937
-0.00301082851
-0.00301152868
-0.00301222987
-0.00301293204
-0.00301363515
-0.00301433916
-0.00301504402
-0.00301574966
-0.00301645602
-0.00301716301
-0.00301787056
-0.00301857858
-0.00301928697
-0.00301999563
-0.00302070448
-0.00302141342
-0.00302212238
-0.00302283129
-0.00302354009
-0.00302424876
-0.00302495726
-0.0030256656
-0.00302637377
-0.00302708181
-0.00302778974
-0.0030284976
-0.00302920541
-0.00302991328
-0.00303062095
-0.00303132989
-0.00303203401
-0.00303275902
-0.00303342196
-0.00303447928
-0.00303514582
-0.00300040563
-0.00300110488
-0.00300180018
-0.00300249358
-0.00300318666
-0.00300387818
-0.00300457055
-0.00300526257
-0.00300595539
-0.00300664879
-0.00300734302
-0.00300803811
-0.00300873419
-0.00300943126
-0.00301012937
-0.00301082851
-0.00301152868
-0.00301222987
-0.00301293204
-0.00301363515
-0.00301433916
-0.00301504402
-0.00301574966
-0.00301645602
-0.00301716301
-0.00301787056
-0.00301857858
-0.00301928697
-0.00301999563
-0.00302070448
-0.00302141342
-0.00302212238
-0.00302283129
-0.00302354009
-0.00302424876
-0.00302495726
-0.0030256656
-0.00302637377
-0.00302708181
-0.00302778974
-0.0030284976
-0.00302920541
-0.00302991328
-0.00303062095
-0.00303132989
-0.00303203401
-0.00303275902
-0.00303342196
-0.00303447928
-0.00303514582
-0.00300040563
-0.00300110488
-0.00300180018
-0.00300249358
-0.00300318666
-0.00300387818
-0.00300457055
-0.00300526257
-0.00300595539
-0.00300664879
-0.00300734302
-0.00300803811
-0.00300873419
-0.00300943126
-0.00301012937
-0.00301082851
-0.00301152868
-0.00301222987
-0.00301293204
-0.00301363515
-0.00301433916
-0.00301504402
-0.00301574966
-0.00301645602
-0.00301716301
-0.00301787056
-0.00301857858
-0.00301928697
-0.00301999563
-0.00302070448
-0.00302141342
-0.00302212238
-0.00302283129
-0.00302354009
-0.00302424876
-0.00302495726
-0.0030256656
-0.00302637377
-0.00302708181
-0.00302778974
-0.0030284976
-0.00302920541
-0.00302991328
-0.00303062095
-0.00303132989
-0.00303203401
-0.00303275902
-0.00303342196
-0.00303447928
-0.00303514582
-0.00300040563
-0.00300110488
-0.00300180018
-0.00300249358
-0.00300318666
-0.00300387818
-0.00300457055
-0.00300526257
-0.00300595539
-0.00300664879
-0.00300734302
-0.00300803811
-0.00300873419
-0.00300943126
-0.00301012937
-0.00301082851
-0.00301152868
-0.00301222987
-0.00301293204
-0.00301363515
-0.00301433916
-0.00301504402
-0.00301574966
-0.00301645602
-0.00301716301
-0.00301787056
-0.00301857858
-0.00301928697
-0.00301999563
-0.00302070448
-0.00302141342
-0.00302212238
-0.00302283129
-0.00302354009
-0.00302424876
-0.00302495726
-0.0030256656
-0.00302637377
-0.00302708181
-0.00302778974
-0.0030284976
-0.00302920541
-0.00302991328
-0.00303062095
-0.00303132989
-0.00303203401
-0.00303275902
-0.00303342196
-0.00303447928
-0.00303514582
-)
-;
-    }
-}
-
-
-// ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/Allclean

@@ -6,6 +6,6 @@ cd "${0%/*}" || exit 1
 # Source clean functions
 . $WM_PROJECT_DIR/bin/tools/CleanFunctions
 
-cleanCase && rm -f moments.eps numberDensity.eps
+cleanCase && rm -rf *.eps 0/d.air 0/uniform/ 0/integerMoment*
 
 #------------------------------------------------------------------------------

test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/system/controlDict

@@ -28,7 +28,7 @@ deltaT          0.01;
 
 writeControl    runTime;
 
-writeInterval   4;
+writeInterval   0.2;
 
 purgeWrite      0;
 
@@ -52,37 +52,38 @@ maxDeltaT       1;
 
 functions
 {
-    numberDensity
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        functionType=numberDensity,
+        coordinateType=volume,
+        funcName=numberDensity
+    )
 
-        writeControl        outputTime;
+    #includeFunc moments
-        writeInterval       1;
+    (
+        populationBalance=bubbles,
+        momentType=integerMoment,
+        coordinateType=volume,
+        order=1
+    )
 
-        setFormat           raw;
+    #includeFunc moments
+    (
+        populationBalance=bubbles,
+        momentType=integerMoment,
+        coordinateType=volume,
+        order=0
+    )
 
-        populationBalance   bubbles;
+    #includeFunc probes
-        functionType        number;
+    (
-        coordinateType      volume;
+        funcName=probes,
-        densityFunction     yes;
+        points=((0.5 0.5 0.5)),
-    }
+        integerMoment0(N,v).bubbles,
-
+        integerMoment1(N,v).bubbles,
-    moments
+        writeControl=writeTime
-    {
+    )
-        type                sizeDistribution;
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
-
-        writeControl        runTime;
-        writeInterval       0.1;
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        functionType        moments;
-        coordinateType      volume;
-        maxOrder            1;
-    }
 }
 
 // ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/negativeDrift/validation/createGraphs

@@ -23,7 +23,7 @@ gnuplot<<EOF
 
     plot '<printf "7.548417545 0\n7.548417545 0.01\n10.62698298 0.01\n10.62698298 0\n"' w l dt 2 lc rgb 'black' t 'Initial condition',\
     '<printf "3.548417545 0\n3.548417545 0.01\n6.626982979 0.01\n6.626982979 0\n"' w l lc rgb 'black' t 'Exact',\
-    '../postProcessing/numberDensity/4/numberDensity.xy' u 1:4 w lp ls 4 lc rgb 'black' t 'Numerical'
+    '../postProcessing/numberDensity/4/numberDensity.xy' u 1:2 w lp ls 4 lc rgb 'black' t 'Numerical'
 
     set output '../moments.eps'
 
@@ -36,9 +36,9 @@ gnuplot<<EOF
     set yrange [0:1.5]
     set key top right
 
-    plot '../postProcessing/moments/0/moments.dat' every 2 u 1:(\$2/0.030785654) w p pt 5 lc rgb 'black' t 'j=0 (Numbers), numerical',\
+    plot '../postProcessing/probes/0/integerMoment0(N,v).bubbles' u 1:(\$2/0.030785654) w p pt 5 lc rgb 'black' t 'j=0 (Numbers), numerical',\
     '<printf "0 1\n6 1\n"' w l lc rgb 'black' t 'j=0 (Numbers), exact',\
-    '../postProcessing/moments/0/moments.dat' u 1:(\$3/0.279433081) every 2 w p pt 6 lc rgb 'black' t 'j=1 (Volume), numerical',\
+    '../postProcessing/probes/0/integerMoment1(N,v).bubbles' u 1:(\$2/0.279433081) w p pt 6 lc rgb 'black' t 'j=1 (Volume), numerical',\
     '<printf "0 1\n6 0.33333\n"' w l dt 2 lc rgb 'black' t 'j=1 (Volume), exact'
 EOF
 

test/multiphase/multiphaseEulerFoam/populationBalance/simultaneousCoalescenceAndBreakup/system/controlDict

@@ -52,39 +52,23 @@ maxDeltaT       1;
 
 functions
 {
-    numberDensity
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        functionType=numberDensity,
+        coordinateType=diameter,
+        normalise=yes,
+        funcName=numberDensity
+    )
 
-        writeControl        outputTime;
+    #includeFunc sizeDistribution
-        writeInterval       1;
+    (
-
+        populationBalance=bubbles,
-        setFormat           raw;
+        functionType=volumeDensity,
-
+        coordinateType=diameter,
-        populationBalance   bubbles;
+        normalise=yes,
-        functionType        number;
+        funcName=volumeDensity
-        coordinateType      diameter;
+    )
-        normalise           yes;
-        densityFunction     yes;
-    }
-
-    volumeDensity
-    {
-        type                sizeDistribution;
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
-
-        writeControl        outputTime;
-        writeInterval       1;
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        functionType        volume;
-        coordinateType      diameter;
-        normalise           yes;
-        densityFunction     yes;
-    }
 }
 
 // ************************************************************************* //

test/multiphase/multiphaseEulerFoam/populationBalance/simultaneousCoalescenceAndBreakup/validation/createGraphs

@@ -28,12 +28,12 @@ gnuplot<<EOF
 
     plot '-' w l lc rgb 'black' t 'NDF, Vanni (2000)',\
     '-' w l dt 2 lc rgb 'black' t 'VDF, Vanni (2000)',\
-    numberDensity every ::0::12 u (\$3/kVConversion):(\$4*kVConversion) w p pt 1 lc rgb 'black' t 'air1',\
+    numberDensity every ::0::12 u (\$1/kVConversion):(\$2*kVConversion) w p pt 1 lc rgb 'black' t 'air1',\
-    numberDensity every ::13::22 u (\$3/kVConversion):(\$4*kVConversion) w p pt 5 lc rgb 'black' t 'air2',\
+    numberDensity every ::13::22 u (\$1/kVConversion):(\$2*kVConversion) w p pt 5 lc rgb 'black' t 'air2',\
-    numberDensity every ::23::30 u (\$3/kVConversion):(\$4*kVConversion) w p pt 9 lc rgb 'black' t 'air3',\
+    numberDensity every ::23::30 u (\$1/kVConversion):(\$2*kVConversion) w p pt 9 lc rgb 'black' t 'air3',\
-    volumeDensity every ::0::12 u (\$3/kVConversion):(\$4*kVConversion) w p pt 2 lc rgb 'black' t 'air1',\
+    volumeDensity every ::0::12 u (\$1/kVConversion):(\$2*kVConversion) w p pt 2 lc rgb 'black' t 'air1',\
-    volumeDensity every ::13::22 u (\$3/kVConversion):(\$4*kVConversion) w p pt 6 lc rgb 'black' t 'air2',\
+    volumeDensity every ::13::22 u (\$1/kVConversion):(\$2*kVConversion) w p pt 6 lc rgb 'black' t 'air2',\
-    volumeDensity every ::23::30 u (\$3/kVConversion):(\$4*kVConversion) w p pt 10 lc rgb 'black' t 'air3'
+    volumeDensity every ::23::30 u (\$1/kVConversion):(\$2*kVConversion) w p pt 10 lc rgb 'black' t 'air3'
     1.000   0.042
     1.460   0.062
     1.950   0.083

tutorials/multiphase/multiphaseEulerFoam/RAS/bubblePipe/system/controlDict

@@ -73,43 +73,27 @@ functions
         )
     )
 
-    probabilityDensity.injection
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=bubbles,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        regionType=cellZone,
+        name=injection,
+        functionType=volumeDensity,
+        coordinateType=diameter,
+        normalise=yes,
+        funcName=probabilityDensity.injection
+    )
 
-        writeControl        outputTime;
+    #includeFunc sizeDistribution
-        writeInterval       1;
+    (
-
+        populationBalance=bubbles,
-        setFormat           raw;
+        regionType=cellZone,
-
+        name=outlet,
-        populationBalance   bubbles;
+        functionType=volumeDensity,
-        regionType          cellZone;
+        coordinateType=diameter,
-        name                injection;
+        normalise=yes,
-        functionType        volume;
+        funcName=probabilityDensity.outlet
-        coordinateType      diameter;
+    )
-        densityFunction     yes;
-        normalise           yes;
-    }
-
-    probabilityDensity.outlet
-    {
-        type                sizeDistribution;
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
-
-        writeControl        outputTime;
-        writeInterval       1;
-
-        setFormat           raw;
-
-        populationBalance   bubbles;
-        regionType          cellZone;
-        name                outlet;
-        functionType        volume;
-        coordinateType      diameter;
-        densityFunction     yes;
-        normalise           yes;
-    }
 }
 
 // ************************************************************************* //

tutorials/multiphase/multiphaseEulerFoam/RAS/bubblePipe/validation/createGraphs

@@ -43,8 +43,8 @@ gnuplot<<EOF
     set ylabel '({/Symbol a}/{/Symbol a}_{tot})/{/Symbol D}d (mm^{-1})'
     set xtics 1
 
-    plot "$injectionFile" u (\$3/1e-3):(\$4/1e3) w histeps lw 3 lc rgb 'green' t 'injection',\
+    plot "$injectionFile" u (\$1/1e-3):(\$2/1e3) w histeps lw 3 lc rgb 'green' t 'injection',\
-         "$outletFile" u (\$3/1e-3):(\$4/1e3) w histeps lw 3 lc rgb 'red' t 'outlet'
+         "$outletFile" u (\$1/1e-3):(\$2/1e3) w histeps lw 3 lc rgb 'red' t 'outlet'
 EOF
 
 #------------------------------------------------------------------------------

tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperse/system/controlDict.orig

@@ -60,7 +60,7 @@ functions
         writeInterval   0.5;
         setFormat       raw;
         log             yes;
-        functionType    volume;
+        functionType    volumeConcentration;
         coordinateType  diameter;
         regionType      cellZone;
         name            x349;
@@ -74,7 +74,7 @@ functions
         writeInterval   0.5;
         setFormat       raw;
         log             yes;
-        functionType    volume;
+        functionType    volumeConcentration;
         coordinateType  diameter;
         regionType      cellZone;
         name            x349_bulk;
@@ -88,7 +88,7 @@ functions
         writeInterval   0.5;
         setFormat       raw;
         log             yes;
-        functionType    volume;
+        functionType    volumeConcentration;
         coordinateType  diameter;
         regionType      cellZone;
         name            x349_wall;

tutorials/multiphase/multiphaseEulerFoam/RAS/wallBoilingPolydisperseTwoGroups/system/controlDict.orig

@@ -60,7 +60,7 @@ functions
         writeInterval   0.5;
         setFormat       raw;
         log             yes;
-        functionType    volume;
+        functionType    volumeConcentration;
         coordinateType  diameter;
         regionType      cellZone;
         name            x349;
@@ -74,7 +74,7 @@ functions
         writeInterval   0.5;
         setFormat       raw;
         log             yes;
-        functionType    volume;
+        functionType    volumeConcentration;
         coordinateType  diameter;
         regionType      cellZone;
         name            x349_bulk;
@@ -88,7 +88,7 @@ functions
         writeInterval   0.5;
         setFormat       raw;
         log             yes;
-        functionType    volume;
+        functionType    volumeConcentration;
         coordinateType  diameter;
         regionType      cellZone;
         name            x349_wall;

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/system/controlDict

@@ -60,26 +60,18 @@ functions
         fields=(TiCl4.vapor O2.vapor Cl2.vapor alpha.particles)
     )
 
-    numberConcentration
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=aggregates,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        regionType=cellZone,
-
+        name=outlet,
-        writeControl        outputTime;
+        functionType=numberDensity,
-        writeInterval       1;
+        coordinateType=projectedAreaDiameter,
-        setFormat           raw;
+        allCoordinates=yes,
-        log                 yes;
+        normalise=yes,
-
+        logTransform=yes,
-        regionType          cellZone;
+        funcName=numberDensity
-        name                outlet;
+    )
-
-        populationBalance   aggregates;
-        functionType        number;
-        coordinateType      projectedAreaDiameter;
-        densityFunction     yes;
-        normalise           yes;
-        geometric           yes;
-    }
 
     #includeFunc writeObjects
     (

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesis/validation/createGraphs

@@ -27,7 +27,7 @@ gnuplot<<EOF
          "$graphFile" u 1:5 axis x1y2 w l t '{/Symbol a}_{particles}'
 EOF
 
-resultFile=../postProcessing/numberConcentration/$time/numberConcentration.xy
+resultFile=../postProcessing/numberDensity/$time/numberDensity.xy
 
 gnuplot<<EOF
     set terminal postscript eps color enhanced font "Helvetica,20"
@@ -43,9 +43,9 @@ gnuplot<<EOF
     set xrange [1e-3:1]
     set yrange [1e-3:3]
 
-    plot "$resultFile" u (\$3/1e-6):4 w lp t 'Size distribution',\
+    plot "$resultFile" u (\$1/1e-6):2 w lp t 'Size distribution',\
-         "$resultFile" u (\$3/1e-6):(\$2/\$1) axis x1y2 w l t '{/Symbol k}',\
+         "$resultFile" u (\$1/1e-6):(\$4/\$3) axis x1y2 w l t '{/Symbol k}',\
-         "$resultFile" u (\$3/1e-6):(6./((6./pi*\$1)**(1./3.))) axis x1y2 w l t '{/Symbol k} lower bound'
+         "$resultFile" u (\$1/1e-6):(6./((6./pi*\$3)**(1./3.))) axis x1y2 w l t '{/Symbol k} lower bound'
 EOF
 
 #------------------------------------------------------------------------------

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/system/controlDict

@@ -60,26 +60,18 @@ functions
         fields=(TiCl4.vapor O2.vapor Cl2.vapor alpha.particles)
     )
 
-    numberConcentration
+    #includeFunc sizeDistribution
-    {
+    (
-        type                sizeDistribution;
+        populationBalance=aggregates,
-        functionObjectLibs  ("libmultiphaseEulerFoamFunctionObjects.so");
+        regionType=cellZone,
-
+        name=outlet,
-        writeControl        outputTime;
+        functionType=numberDensity,
-        writeInterval       1;
+        coordinateType=projectedAreaDiameter,
-        setFormat           raw;
+        allCoordinates=yes,
-        log                 yes;
+        normalise=yes,
-
+        logTransform=yes,
-        regionType          cellZone;
+        funcName=numberDensity
-        name                outlet;
+    )
-
-        populationBalance   aggregates;
-        functionType        number;
-        coordinateType      projectedAreaDiameter;
-        densityFunction     yes;
-        normalise           yes;
-        geometric           yes;
-    }
 
     #includeFunc writeObjects
     (

tutorials/multiphase/multiphaseEulerFoam/laminar/titaniaSynthesisSurface/validation/createGraphs

@@ -27,7 +27,7 @@ gnuplot<<EOF
          "$graphFile" u 1:5 axis x1y2 w l t '{/Symbol a}_{particles}'
 EOF
 
-resultFile=../postProcessing/numberConcentration/$time/numberConcentration.xy
+resultFile=../postProcessing/numberDensity/$time/numberDensity.xy
 
 gnuplot<<EOF
     set terminal postscript eps color enhanced font "Helvetica,20"
@@ -43,9 +43,9 @@ gnuplot<<EOF
     set xrange [1e-3:1]
     set yrange [1e-3:3]
 
-    plot "$resultFile" u (\$3/1e-6):4 w lp t 'Size distribution',\
+    plot "$resultFile" u (\$1/1e-6):2 w lp t 'Size distribution',\
-         "$resultFile" u (\$3/1e-6):(\$2/\$1) axis x1y2 w l t '{/Symbol k}',\
+         "$resultFile" u (\$1/1e-6):(\$4/\$3) axis x1y2 w l t '{/Symbol k}',\
-         "$resultFile" u (\$3/1e-6):(6./((6./pi*\$1)**(1./3.))) axis x1y2 w l t '{/Symbol k} lower bound'
+         "$resultFile" u (\$1/1e-6):(6./((6./pi*\$3)**(1./3.))) axis x1y2 w l t '{/Symbol k} lower bound'
 EOF
 
 #------------------------------------------------------------------------------