ENH: selection mechanism for cloudInfo function object (#2390)

- can restrict calculation of D32 and other spray properties to a
  subset of parcels. Uses a predicate selection mechanism similar to
  vtkCloud etc.

ENH: code cleanup in scalar predicates

- pass by value not reference in predicates
- additional assign() method to refactor common code

src/OpenFOAM/primitives/predicates/scalar/scalarPredicates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -182,30 +182,8 @@ Foam::predicates::scalars::scalars
 (
     std::initializer_list<std::pair<word, scalar>> entries
 )
-:
-    List<unary>(entries.size())
 {
-    // Access
+    assign(entries);
-    const auto get0 =
-        [](const std::pair<word,scalar>& entry) { return entry.first; };
-    const auto get1 =
-        [](const std::pair<word,scalar>& entry) { return entry.second; };
-
-    // Check for bad/unknown operations
-    if (hasBadEntries(entries, get0))
-    {
-        printBadEntries(FatalErrorInFunction, entries, get0, get1)
-            << exit(FatalError);
-    }
-
-    // Appears to be good
-    auto& list = *this;
-    label idx = 0;
-    for (const auto& entry : entries)
-    {
-        list[idx] = predicates::scalars::operation(entry);
-        ++idx;
-    }
 }
 
 
@@ -213,14 +191,29 @@ Foam::predicates::scalars::scalars
 (
     const UList<Tuple2<word, scalar>>& entries
 )
-:
-    List<unary>(entries.size())
 {
+    assign(entries);
+}
+
+
+Foam::predicates::scalars::scalars(Istream& is)
+{
+    is >> *this;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::predicates::scalars::assign
+(
+    std::initializer_list<std::pair<word, scalar>> entries
+)
+{
+    typedef std::pair<word, scalar> tuple_type;
+
     // Access
-    const auto get0 =
+    const auto get0 = [](const tuple_type& entry) { return entry.first; };
-        [](const Tuple2<word,scalar>& entry) { return entry.first(); };
+    const auto get1 = [](const tuple_type& entry) { return entry.second; };
-    const auto get1 =
-        [](const Tuple2<word,scalar>& entry) { return entry.second(); };
 
     // Check for bad/unknown operations
     if (hasBadEntries(entries, get0))
@@ -229,27 +222,48 @@ Foam::predicates::scalars::scalars
             << exit(FatalError);
     }
 
-    // Appears to be good
+    // Appears to be good, fill the list
-    auto& list = *this;
+    this->resize_nocopy(entries.size());
-    label idx = 0;
+    auto iter = this->begin();
-    for (const auto& entry : entries)
+
+    for (const tuple_type& entry : entries)
     {
-        list[idx] = predicates::scalars::operation(entry);
+        *iter = predicates::scalars::operation(entry);
-        ++idx;
+        ++iter;
     }
 }
 
 
-Foam::predicates::scalars::scalars(Istream& is)
+void Foam::predicates::scalars::assign
-:
+(
-    List<unary>()
+    const UList<Tuple2<word, scalar>>& entries
+)
 {
-    is >> *this;
+    typedef Tuple2<word, scalar> tuple_type;
+
+    // Access
+    const auto get0 = [](const tuple_type& entry) { return entry.first(); };
+    const auto get1 = [](const tuple_type& entry) { return entry.second(); };
+
+    // Check for bad/unknown operations
+    if (hasBadEntries(entries, get0))
+    {
+        printBadEntries(FatalErrorInFunction, entries, get0, get1)
+            << exit(FatalError);
+    }
+
+    // Appears to be good, fill the list
+    this->resize_nocopy(entries.size());
+    auto iter = this->begin();
+
+    for (const tuple_type& entry : entries)
+    {
+        *iter = predicates::scalars::operation(entry);
+        ++iter;
+    }
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
 Foam::label Foam::predicates::scalars::find
 (
     const scalar value,
@@ -310,29 +324,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Foam::predicates::scalars& list)
     // Read tuples
     List<Tuple2<word, scalar>> entries(is);
 
-    // Access
+    list.assign(entries);
-    const auto get0 =
-        [](const Tuple2<word,scalar>& entry) { return entry.first(); };
-    const auto get1 =
-        [](const Tuple2<word,scalar>& entry) { return entry.second(); };
-
-
-    // Check for bad/unknown operations
-    if (hasBadEntries(entries, get0))
-    {
-        printBadEntries(FatalIOErrorInFunction(is), entries, get0, get1)
-            << exit(FatalIOError);
-    }
-
-    // Appears to be good
-    list.resize(entries.size());
-
-    label idx = 0;
-    for (const Tuple2<word, scalar>& entry : entries)
-    {
-        list[idx] = predicates::scalars::operation(entry);
-        ++idx;
-    }
 
     return is;
 }

src/OpenFOAM/primitives/predicates/scalar/scalarPredicates.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -256,6 +256,14 @@ public:
 
     // Member Functions
 
+        //- Assign new predicates from an initializer list of
+        //- (opName opValue) tuples
+        void assign(std::initializer_list<std::pair<word, scalar>> entries);
+
+        //- Assign new predicates from a list of (opName opValue) tuples
+        void assign(const UList<Tuple2<word, scalar>>& entries);
+
+
     // Search
 
         //- Index of the first match for the value.
@@ -279,23 +287,23 @@ public:
         //- Match any condition in the list.
         //
         //  \return True if the value matches any condition in the list.
-        inline bool match(const scalar& value) const;
+        inline bool match(const scalar value) const;
 
         //- Match any condition in the list.
         //
         //  \return True if the value matches any condition in the list.
-        inline bool matchAny(const scalar& value) const;
+        inline bool matchAny(const scalar value) const;
 
         //- Match all conditions in the list.
         //
         //  \return True if the value matches all conditions in the list.
-        inline bool matchAll(const scalar& value) const;
+        inline bool matchAll(const scalar value) const;
 
         //- Extract list indices for all matches.
         //
         //  \return The locations (indices) in the input list where match()
         //      is true
-        inline labelList matching(const scalar& value) const;
+        inline labelList matching(const scalar value) const;
 
         //- Extract list indices for all matches.
         //
@@ -313,7 +321,7 @@ public:
     // Member Operators
 
         //- Identical to found(), match(), for use as a predicate.
-        inline bool operator()(const scalar& value) const;
+        inline bool operator()(const scalar value) const;
 };
 
 

src/OpenFOAM/primitives/predicates/scalar/scalarPredicatesI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -70,13 +70,13 @@ inline bool Foam::predicates::scalars::found
 }
 
 
-inline bool Foam::predicates::scalars::match(const scalar& value) const
+inline bool Foam::predicates::scalars::match(const scalar value) const
 {
     return this->matchAny(value);
 }
 
 
-inline bool Foam::predicates::scalars::matchAny(const scalar& value) const
+inline bool Foam::predicates::scalars::matchAny(const scalar value) const
 {
     for (const unary& test : *this)
     {
@@ -90,7 +90,7 @@ inline bool Foam::predicates::scalars::matchAny(const scalar& value) const
 }
 
 
-inline bool Foam::predicates::scalars::matchAll(const scalar& value) const
+inline bool Foam::predicates::scalars::matchAll(const scalar value) const
 {
     for (const unary& test : *this)
     {
@@ -106,7 +106,7 @@ inline bool Foam::predicates::scalars::matchAll(const scalar& value) const
 
 inline Foam::labelList Foam::predicates::scalars::matching
 (
-    const scalar& value
+    const scalar value
 ) const
 {
     labelList indices(this->size());
@@ -154,7 +154,7 @@ inline Foam::labelList Foam::predicates::scalars::matching
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline bool Foam::predicates::scalars::operator()(const scalar& value) const
+inline bool Foam::predicates::scalars::operator()(const scalar value) const
 {
     return this->found(value);
 }

src/functionObjects/lagrangian/cloudInfo/cloudInfo.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +27,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cloudInfo.H"
+#include "cloud.H"
 #include "kinematicCloud.H"
 #include "dictionary.H"
+#include "mathematicalConstants.H"
 #include "PstreamReduceOps.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -74,8 +76,10 @@ Foam::functionObjects::cloudInfo::cloudInfo
     const dictionary& dict
 )
 :
-    regionFunctionObject(name, runTime, dict),
+    functionObjects::regionFunctionObject(name, runTime, dict),
-    logFiles(obr_, name, dict)
+    functionObjects::logFiles(obr_, name, dict),
+    verbose_(false),
+    onExecute_(false)
 {
     read(dict);
 }
@@ -85,55 +89,175 @@ Foam::functionObjects::cloudInfo::cloudInfo
 
 bool Foam::functionObjects::cloudInfo::read(const dictionary& dict)
 {
+    parcelSelect_.clear();
+    verbose_ = false;
+    onExecute_ = false;
+
     if (regionFunctionObject::read(dict) && logFiles::read(dict))
     {
         logFiles::resetNames(dict.get<wordList>("clouds"));
 
         Info<< type() << " " << name() << ": ";
-        if (writeToFile() && names().size())
+        if (names().size())
         {
             Info<< "applying to clouds:" << nl;
-            forAll(names(), cloudi)
+            for (const word& cldName : names())
             {
-                Info<< "    " << names()[cloudi] << nl;
+                Info<< "    " << cldName << nl;
-                writeFileHeader(files(cloudi));
             }
             Info<< endl;
+
+            // Actions to define selection
+            parcelSelect_ = dict.subOrEmptyDict("selection");
+
+            verbose_ = dict.getOrDefault("verbose", false);
+            onExecute_ = dict.getOrDefault("sampleOnExecute", false);
         }
         else
         {
             Info<< "no clouds to be processed" << nl << endl;
         }
+
+        if (writeToFile())
+        {
+            forAll(names(), cloudi)
+            {
+                writeFileHeader(files(cloudi));
+            }
+        }
     }
 
     return true;
 }
 
 
-bool Foam::functionObjects::cloudInfo::execute()
+bool Foam::functionObjects::cloudInfo::performAction(unsigned request)
 {
-    return true;
+    if (!request || names().empty())
-}
+    {
+        return true;
+    }
 
-
-bool Foam::functionObjects::cloudInfo::write()
-{
     forAll(names(), cloudi)
     {
+        // The reported quantities
+        label nTotParcels = 0;
+        scalar totMass = 0, Dmax = 0, D10 = 0, D32 = 0;
+        bool applyFilter = false;
+
         const word& cloudName = names()[cloudi];
 
-        const kinematicCloud& cloud =
+        const auto* kinCloudPtr = obr_.cfindObject<kinematicCloud>(cloudName);
-            obr_.lookupObject<kinematicCloud>(cloudName);
 
-        const label nTotParcels =
+        if (!kinCloudPtr)
-            returnReduce(cloud.nParcels(), sumOp<label>());
+        {
+            // Safety
+            continue;
+        }
 
-        const scalar totMass =
+        const auto& kinCloud = *kinCloudPtr;
-            returnReduce(cloud.massInSystem(), sumOp<scalar>());
+        const auto* plainCloudPtr = isA<cloud>(kinCloud);
 
-        const scalar Dmax = cloud.Dmax();
+        if (!parcelSelect_.empty() && plainCloudPtr)
-        const scalar D10 = cloud.Dij(1, 0);
+        {
-        const scalar D32 = cloud.Dij(3, 2);
+            const auto& plainCloud = *plainCloudPtr;
+
+            // Filtering - simply use cloud methods
+
+            objectRegistry obrTmp
+            (
+                IOobject
+                (
+                    "tmp::cloudInfo::" + cloudName,
+                    obr_.time().constant(),
+                    obr_,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                )
+            );
+
+            plainCloud.writeObjects(obrTmp);
+
+            // Apply output filter (for the current cloud)
+            applyFilter = calculateFilter(obrTmp, log);
+
+            // Expected/required fields
+            const auto* diamFldPtr = obrTmp.cfindObject<IOField<scalar>>("d");
+            const auto* rhoFldPtr = obrTmp.cfindObject<IOField<scalar>>("rho");
+            const auto* nParticleFldPtr =
+                obrTmp.cfindObject<IOField<scalar>>("nParticle");
+
+            do
+            {
+                #undef  doLocalCode
+                #define doLocalCode(FldPtr, FldName)                          \
+                if (applyFilter && !FldPtr)                                   \
+                {                                                             \
+                    WarningInFunction                                         \
+                        << "Missing \"" << #FldName                           \
+                        << "\" field - disabling filter" << nl;               \
+                    applyFilter = false;                                      \
+                    break;                                                    \
+                }
+
+                doLocalCode(diamFldPtr, d);
+                doLocalCode(rhoFldPtr, rho);
+                doLocalCode(nParticleFldPtr, nParticle);
+                #undef doLocalCode
+            }
+            while (false);
+
+            if (applyFilter)
+            {
+                // Filtered. Need to do everything by hand!
+
+                const auto& diams = *diamFldPtr;
+                const auto& rhos = *rhoFldPtr;
+                const auto& nParts = *nParticleFldPtr;
+
+                FixedList<scalar, 4> Dsums(Zero);
+
+                for (const label particlei : parcelAddr_)
+                {
+                    ++nTotParcels;
+
+                    const scalar d = diams[particlei];
+                    const scalar rho = rhos[particlei];
+                    const scalar np = nParts[particlei];
+
+                    totMass += np*rho*pow3(d);
+                    Dmax = max(Dmax, d);
+
+                    Dsums[0] += np;
+                    Dsums[1] += np*(d);
+                    Dsums[2] += np*(sqr(d));
+                    Dsums[3] += np*(pow3(d));
+                }
+
+                reduce(nTotParcels, sumOp<label>());
+                reduce(totMass, sumOp<scalar>());
+                reduce(Dmax, maxOp<scalar>());
+                reduce(Dsums, sumOp<scalar>());
+
+                totMass *= (constant::mathematical::pi/6.0);
+
+                Dmax = max(0, Dmax);
+                D10 = Dsums[1]/(max(Dsums[0], VSMALL));
+                D32 = Dsums[3]/(max(Dsums[2], VSMALL));
+            }
+        }
+
+        if (!applyFilter)
+        {
+            // No filter - use regular cloud methods
+            nTotParcels = returnReduce(kinCloud.nParcels(), sumOp<label>());
+            totMass = returnReduce(kinCloud.massInSystem(), sumOp<scalar>());
+
+            Dmax = kinCloud.Dmax();
+            D10 = kinCloud.Dij(1, 0);
+            D32 = kinCloud.Dij(3, 2);
+        }
 
         Log << type() << " " << name() <<  " write:" << nl
             << "    number of parcels : " << nTotParcels << nl
@@ -143,7 +267,7 @@ bool Foam::functionObjects::cloudInfo::write()
             << "    D32 diameter      : " << D32 << nl
             << endl;
 
-        if (writeToFile())
+        if ((request & ACTION_WRITE) && writeToFile())
         {
             auto& os = files(cloudi);
 
@@ -162,4 +286,21 @@ bool Foam::functionObjects::cloudInfo::write()
 }
 
 
+bool Foam::functionObjects::cloudInfo::execute()
+{
+    if (onExecute_)
+    {
+        return performAction(ACTION_ALL & ~ACTION_WRITE);
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::cloudInfo::write()
+{
+    return performAction(ACTION_ALL);
+}
+
+
 // ************************************************************************* //

src/functionObjects/lagrangian/cloudInfo/cloudInfo.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -58,6 +58,8 @@ Usage
         Property     | Description                          | Required | Default
         type         | type name: cloudInfo                 | yes |
         clouds       | list of clouds names to process      | yes |
+        selection    | Parcel selection control         | no  | empty-dict
+        sampleOnExecute| Sample/report (on execute) without writing | no | false
     \endtable
 
     The output data of each cloud is written to a file named \<cloudName\>.dat
@@ -72,11 +74,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef functionObjects_cloudInfo_H
+#ifndef Foam_functionObjects_cloudInfo_H
-#define functionObjects_cloudInfo_H
+#define Foam_functionObjects_cloudInfo_H
 
 #include "regionFunctionObject.H"
 #include "logFiles.H"
+#include "parcelSelectionDetail.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -91,13 +94,31 @@ namespace functionObjects
 
 class cloudInfo
 :
-    public regionFunctionObject,
+    public functionObjects::regionFunctionObject,
-    public logFiles
+    public functionObjects::logFiles,
+    public Foam::Detail::parcelSelection
 {
 protected:
 
-        //- Switch to send output to Info as well
+    // Data Types
-        Switch log_;
+
+        //- Local control for sampling actions
+        enum sampleActionType : unsigned
+        {
+            ACTION_NONE  = 0,
+            ACTION_WRITE = 0x1,
+            ACTION_STORE = 0x2,
+            ACTION_ALL = 0xF
+        };
+
+
+    // Protected Data
+
+        //- Additional verbosity
+        bool verbose_;
+
+        //- Perform sample actions on execute as well
+        bool onExecute_;
 
         //- List of cloud names
         wordList cloudNames_;
@@ -105,11 +126,14 @@ protected:
         //- Output file per cloud
         PtrList<OFstream> filePtrs_;
 
+
     // Protected Member Functions
 
         //- File header information
         virtual void writeFileHeader(Ostream& os) const;
 
+        //- Perform operation report/write
+        bool performAction(unsigned request);
 
         //- No copy construct
         cloudInfo(const cloudInfo&) = delete;